diff --git a/Docs/sphinx_documentation/source/LinearSolvers.rst b/Docs/sphinx_documentation/source/LinearSolvers.rst index 6d9d6b782a6..a68336a200a 100644 --- a/Docs/sphinx_documentation/source/LinearSolvers.rst +++ b/Docs/sphinx_documentation/source/LinearSolvers.rst @@ -456,11 +456,11 @@ residual correction form of the original problem. To build Hypre, follow the nex To use hypre, one must include ``amrex/Src/Extern/HYPRE`` in the build system. For examples of using hypre, we refer the reader to -`ABecLaplacian`_ or `Nodal Projection EB`_. +`ABecLaplacian`_ or `NodeTensorLap`_. .. _`ABecLaplacian`: https://amrex-codes.github.io/amrex/tutorials_html/LinearSolvers_Tutorial.html -.. _`Nodal Projection EB`: https://amrex-codes.github.io/amrex/tutorials_html/LinearSolvers_Tutorial.html +.. _`NodeTensorLap`: https://amrex-codes.github.io/amrex/tutorials_html/LinearSolvers_Tutorial.html The following parameter should be set to True if the problem to be solved has a singular matrix. In this case, the solution is only defined to within a constant. Setting this parameter to True @@ -525,309 +525,6 @@ To use PETSc, one must include ``amrex/Src/Extern/PETSc`` in the build system. For an example of using PETSc, we refer the reader to the tutorial, `ABecLaplacian`_. -MAC Projection -========================= - -Some codes define a velocity field :math:`U = (u,v,w)` on faces, i.e. -:math:`u` is defined on x-faces, :math:`v` is defined on y-faces, -and :math:`w` is defined on z-faces. We refer to the exact projection -of this velocity field as a MAC projection, in which we solve - -.. math:: - - D( \beta \nabla \phi) = D(U^*) - S - -for :math:`\phi` and then set - -.. math:: - - U = U^* - \beta \nabla \phi - - -where :math:`U^*` is a vector field (typically velocity) that we want to satisfy -:math:`D(U) = S`. For incompressible flow, :math:`S = 0`. - -The MacProjection class can be defined and used to perform the MAC projection without explicitly -calling the solver directly. In addition to solving the variable coefficient Poisson equation, -the MacProjector internally computes the divergence of the vector field, :math:`D(U^*)`, -to compute the right-hand-side, and after the solve, subtracts the weighted gradient term to -make the vector field result satisfy the divergence constraint. - -In the simplest form of the call, :math:`S` is assumed to be zero and does not need to be specified. -Typically, the user does not allocate the solution array, but it is also possible to create and pass -in the solution array and have :math:`\phi` returned as well as :math:`U`. - -Caveat: Currently the MAC projection only works when the base level covers the full domain; it does -not yet have the interface to pass boundary conditions for a fine level that come from coarser data. - -Also note that any Dirichlet or Neumann boundary conditions at domain boundaries -are assumed to be homogeneous. The call to the :cpp:`MLLinOp` member function -:cpp:`setLevelBC` occurs inside the MacProjection class; one does not need to call that -explicitly when using the MacProjection class. - -The code below is taken from the file -``Tutorials/LinearSolvers/MAC_Projection_EB/main.cpp`` in `MAC Projection EB`_ and demonstrates how to set up -the MACProjector object and use it to perform a MAC projection. - -.. _`MAC Projection EB`: https://amrex-codes.github.io/amrex/tutorials_html/LinearSolvers_Tutorial.html - -.. highlight:: c++ - -:: - - EBFArrayBoxFactory factory(eb_level, geom, grids, dmap, ng_ebs, ebs); - - // allocate face-centered velocities and face-centered beta coefficient - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - vel[idim].define (amrex::convert(grids,IntVect::TheDimensionVector(idim)), dmap, 1, 1, - MFInfo(), factory); - beta[idim].define(amrex::convert(grids,IntVect::TheDimensionVector(idim)), dmap, 1, 0, - MFInfo(), factory); - beta[idim].setVal(1.0); // set beta to 1 - } - - // If we want to use phi elsewhere, we must create an array in which to return the solution - // MultiFab phi_inout(grids, dmap, 1, 1, MFInfo(), factory); - - // If we want to supply a non-zero S we must allocate and fill it outside the solver - // MultiFab S(grids, dmap, 1, 0, MFInfo(), factory); - // Set S here ... - - // set initial velocity to U=(1,0,0) - AMREX_D_TERM(vel[0].setVal(1.0);, - vel[1].setVal(0.0);, - vel[2].setVal(0.0);); - - LPInfo lp_info; - - // If we want to use hypre to solve the full problem we do not need to coarsen the GMG stencils - if (use_hypre_as_full_solver) - lp_info.setMaxCoarseningLevel(0); - - MacProjector macproj({amrex::GetArrOfPtrs(vel)}, // face-based velocity - {amrex::GetArrOfConstPtrs(beta)}, // beta - {geom}, // the geometry object - lp_info); // structure for passing info to the operator - - // Here we specify the desired divergence S - // MacProjector macproj({amrex::GetArrOfPtrs(vel)}, // face-based velocity - // {amrex::GetArrOfConstPtrs(beta)}, // beta - // {geom}, // the geometry object - // lp_info, // structure for passing info to the operator - // {&S}); // defines the specified RHS divergence - - // Set bottom-solver to use hypre instead of native BiCGStab - if (use_hypre_as_full_solver || use_hypre_as_bottom_solver) - macproj.setBottomSolver(MLMG::BottomSolver::hypre); - - // Set boundary conditions. - // Here we use Neumann on the low x-face, Dirichlet on the high x-face, - // and periodic in the other two directions - // (the first argument is for the low end, the second is for the high end) - // Note that Dirichlet and Neumann boundary conditions are assumed to be homogeneous. - macproj.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann, - LinOpBCType::Periodic, - LinOpBCType::Periodic)}, - {AMREX_D_DECL(LinOpBCType::Dirichlet, - LinOpBCType::Periodic, - LinOpBCType::Periodic)}); - - macproj.setVerbose(mg_verbose); - macproj.setBottomVerbose(bottom_verbose); - - // Define the relative tolerance - Real reltol = 1.e-8; - - // Define the absolute tolerance; note that this argument is optional - Real abstol = 1.e-15; - - // Solve for phi and subtract from the velocity to make it divergence-free - // Note that when we build with USE_EB = TRUE, we must specify whether the velocities live - // at face centers (MLMG::Location::FaceCenter) or face centroids (MLMG::Location::FaceCentroid) - macproj.project(reltol,abstol,MLMG::Location::FaceCenter); - - // If we want to use phi elsewhere, we can pass in an array in which to return the solution - // macproj.project({&phi_inout},reltol,abstol,MLMG::Location::FaceCenter); - -See the `MAC Projection EB`_ tutorial for the complete working example. - -Nodal Projection -================ - -Some codes define a velocity field :math:`U = (u,v,w)` with all -components co-located on cell centers. The nodal solver in AMReX -can be used to compute an approximate projection of the cell-centered -velocity field, with pressure and velocity divergence defined on nodes. -When we use the nodal solver this way, and subtract only the cell average -of the gradient from the velocity, it is effectively an approximate projection. - -As with the MAC projection, consider that we want to solve - -.. math:: - - D( \beta \nabla \phi) = D(U^*) - S - -for :math:`\phi` and then set - -.. math:: - - U = U^* - \beta \nabla \phi - -where :math:`U^*` is a vector field defined on cell centers and we want to satisfy -:math:`D(U) = S`. For incompressible flow, :math:`S = 0`. - -Currently this nodal approximate projection does not exist in a separate -operator like the MAC projection; instead we demonstrate below the steps needed -to compute the approximate projection. This means we must compute explicitly the -right-hand-side , including the the divergence of the vector field, :math:`D(U^*)`, -solve the variable coefficient Poisson equation, then subtract the weighted -gradient term to make the vector field result satisfy the divergence constraint. - -.. highlight:: c++ - -:: - - // - // Given a cell-centered velocity (vel) field, a cell-centered - // scalar field (sigma) field, and a source term S (either node- - // or cell-centered )solve: - // - // div( sigma * grad(phi) ) = div(vel) - S - // - // and then perform the projection: - // - // vel = vel - sigma * grad(phi) - // - - // - // Create the EB factory - // - EBFArrayBoxFactory factory(eb_level, geom, grids, dmap, ng_ebs, ebs); - - // - // Create the cell-centered velocity field we want to project - // - MultiFab vel(grids, dmap, AMREX_SPACEDIM, 1, MFInfo(), factory); - - // Set velocity field to (1,0,0) including ghost cells for this example - vel.setVal(1.0, 0, 1, 1); - vel.setVal(0.0, 1, AMREX_SPACEDIM-1, 1); - - // - // Setup linear operator, AKA the nodal Laplacian - // - LPInfo lp_info; - - // If we want to use hypre to solve the full problem we do not need to coarsen the GMG stencils - // if (use_hypre_as_full_solver) - // lp_info.setMaxCoarseningLevel(0); - - MLNodeLaplacian matrix({geom}, {grids}, {dmap}, lp_info, - Vector{&factory}); - - // Set boundary conditions. - // Here we use Neumann on the low x-face, Dirichlet on the high x-face, - // and periodic in the other two directions - // (the first argument is for the low end, the second is for the high end) - // Note that Dirichlet boundary conditions are assumed to be homogeneous (i.e. phi = 0) - matrix.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann, - LinOpBCType::Periodic, - LinOpBCType::Periodic)}, - {AMREX_D_DECL(LinOpBCType::Dirichlet, - LinOpBCType::Periodic, - LinOpBCType::Periodic)}); - - // Set matrix attributes to be used by MLMG solver - matrix.setGaussSeidel(true); - matrix.setHarmonicAverage(false); - - // - // Compute RHS - // - // NOTE: it's up to the user to compute the RHS. as opposed - // to the MAC projection case !!! - // - // NOTE: do this operation AFTER setting up the linear operator so - // that compRHS method can be used - // - - // RHS is nodal - const BoxArray & nd_grids = amrex::convert(grids, IntVect{1,1,1}); // nodal grids - - // MultiFab to host RHS - MultiFab rhs(nd_grids, dmap, 1, 1, MFInfo(), factory); - - // Cell-centered contributions to RHS - MultiFab S_cc(grids, dmap, 1, 1, MFInfo(), factory); - S_cc.setVal(0.0); // Set it to zero for this example - - // Node-centered contributions to RHS - MultiFab S_nd(nd_grids, dmap, 1, 1, MFInfo(), factory); - S_nd.setVal(0.0); // Set it to zero for this example - - // Compute RHS -- vel must be cell-centered - matrix.compRHS({&rhs}, {&vel}, {&S_nd}, {&S_cc}); - - // - // Create the cell-centered sigma field and set it to 1 for this example - // - MultiFab sigma(grids, dmap, 1, 1, MFInfo(), factory); - sigma.setVal(1.0); - - // Set sigma - matrix.setSigma(0, sigma); - - // - // Create node-centered phi - // - MultiFab phi(nd_grids, dmap, 1, 1, MFInfo(), factory); - phi.setVal(0.0); - - // - // Setup MLMG solver - // - MLMG nodal_solver(matrix); - - // We can specify the maximum number of iterations - nodal_solver.setMaxIter(mg_maxiter); - nodal_solver.setBottomMaxIter(mg_bottom_maxiter); - - nodal_solver.setVerbose(mg_verbose); - nodal_solver.setBottomVerbose(mg_bottom_verbose); - - // Set bottom-solver to use hypre instead of native BiCGStab - // ( we could also have set this to cg, bicgcg, cgbicg) - // if (use_hypre_as_full_solver || use_hypre_as_bottom_solver) - // nodal_solver.setBottomSolver(MLMG::BottomSolver::hypre); - - // Define the relative tolerance - Real reltol = 1.e-8; - - // Define the absolute tolerance; note that this argument is optional - Real abstol = 1.e-15; - - // - // Solve div( sigma * grad(phi) ) = RHS - // - nodal_solver.solve( {&phi}, {&rhs}, reltol, abstol); - - // - // Create cell-centered MultiFab to hold value of -sigma*grad(phi) at cell-centers - // - // - MultiFab fluxes(grids, dmap, AMREX_SPACEDIM, 1, MFInfo(), factory); - fluxes.setVal(0.0); - - // Get fluxes from solver - nodal_solver.getFluxes( {&fluxes} ); - - // - // Apply projection explicitly -- vel = vel - sigma * grad(phi) - // - MultiFab::Add( *vel, *fluxes, 0, 0, AMREX_SPACEDIM, 0); - -See the `Nodal Projection EB`_ tutorial for the complete working example. - Tensor Solve ============ diff --git a/GNUmakefile.in b/GNUmakefile.in index 098faa19274..8a6ce69df09 100644 --- a/GNUmakefile.in +++ b/GNUmakefile.in @@ -18,7 +18,7 @@ ifeq ($(USE_FORTRAN_INTERFACE),TRUE) Pdirs += F_Interfaces/Base F_Interfaces/Octree F_Interfaces/AmrCore endif ifeq ($(USE_LINEAR_SOLVERS),TRUE) - Pdirs += LinearSolvers/MLMG LinearSolvers/Projections + Pdirs += LinearSolvers/MLMG ifeq ($(USE_FORTRAN_INTERFACE),TRUE) Pdirs += F_Interfaces/LinearSolvers endif diff --git a/Src/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt index 2c084bf77cf..38ed7d363a2 100644 --- a/Src/LinearSolvers/CMakeLists.txt +++ b/Src/LinearSolvers/CMakeLists.txt @@ -2,7 +2,6 @@ # Sources in subdirectory MLMG # target_include_directories(amrex PUBLIC $) -target_include_directories(amrex PUBLIC $) target_sources(amrex PRIVATE @@ -53,10 +52,6 @@ target_sources(amrex MLMG/AMReX_MLTensorOp_grad.cpp MLMG/AMReX_MLTensor_K.H MLMG/AMReX_MLTensor_${AMReX_SPACEDIM}D_K.H - Projections/AMReX_MacProjector.H - Projections/AMReX_MacProjector.cpp - Projections/AMReX_NodalProjector.H - Projections/AMReX_NodalProjector.cpp ) if (AMReX_SPACEDIM EQUAL 3) diff --git a/Src/LinearSolvers/Projections/AMReX_MacProjector.H b/Src/LinearSolvers/Projections/AMReX_MacProjector.H deleted file mode 100644 index 98c084ff6b9..00000000000 --- a/Src/LinearSolvers/Projections/AMReX_MacProjector.H +++ /dev/null @@ -1,195 +0,0 @@ -#ifndef AMREX_MAC_PROJECTOR_H_ -#define AMREX_MAC_PROJECTOR_H_ -#include - -#include -#include -#include - -#ifdef AMREX_USE_EB -#include -#endif - -namespace amrex { - -class MacProjector -{ -public: - MacProjector( - const Vector& a_geom, -#ifndef AMREX_USE_EB - MLMG::Location a_umac_loc = MLMG::Location::FaceCenter, - MLMG::Location a_beta_loc = MLMG::Location::FaceCenter, - MLMG::Location a_phi_loc = MLMG::Location::CellCenter, - MLMG::Location a_divu_loc = MLMG::Location::CellCenter -#else - MLMG::Location a_umac_loc, - MLMG::Location a_beta_loc, - MLMG::Location a_phi_loc, - MLMG::Location a_divu_loc = MLMG::Location::CellCenter -#endif - ); - - // - // Constructors - // - MacProjector (const Vector >& a_umac, - MLMG::Location a_umac_loc, - const Vector >& a_beta, - MLMG::Location a_beta_loc, - MLMG::Location a_phi_loc, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const Vector& a_divu, - MLMG::Location a_divu_loc, - const Vector& a_overset_mask = {}); - - MacProjector (const Vector >& a_umac, - MLMG::Location a_umac_loc, - const Vector >& a_beta, - MLMG::Location a_beta_loc, - MLMG::Location a_phi_loc, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const Vector& a_divu = {}) - : MacProjector(a_umac, a_umac_loc, a_beta, a_beta_loc, a_phi_loc, - a_geom, a_lpinfo, a_divu, MLMG::Location::CellCenter) {} - -#ifndef AMREX_USE_EB - MacProjector (const Vector >& a_umac, - const Vector >& a_beta, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const Vector& a_divu = {}) - : MacProjector(a_umac, MLMG::Location::FaceCenter, - a_beta, MLMG::Location::FaceCenter, MLMG::Location::CellCenter, - a_geom, a_lpinfo, a_divu, MLMG::Location::CellCenter) {} - - MacProjector (const Vector >& a_umac, - const Vector >& a_beta, - const Vector& a_geom, - const Vector& a_divu = {}) - : MacProjector(a_umac, MLMG::Location::FaceCenter, - a_beta, MLMG::Location::FaceCenter, MLMG::Location::CellCenter, - a_geom, LPInfo(), a_divu, MLMG::Location::CellCenter) {} - - MacProjector (const Vector >& a_umac, - const Real a_const_beta, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const Vector& a_overset_mask = {}, - const Vector& a_divu = {}); -#endif - - /** Initialize the underlying linear operator and MLMG instances - */ - void initProjector ( - const LPInfo& a_lpinfo, - const Vector >& a_beta, - const Vector& a_overset_mask = {}); - - //! Update Bcoeffs for the linear operator - void updateBeta (const Vector >&); - -#ifndef AMREX_USE_EB - void initProjector (Vector const& a_grids, - Vector const& a_dmap, - const LPInfo& a_lpinfo, Real const a_const_beta, - const Vector& a_overset_mask = {}); - - void updateBeta (Real a_const_beta); -#endif - - //! Set Umac before calling the projection step - void setUMAC(const Vector >&); - - //! Set div(U) - void setDivU(const Vector&); - - // - // Methods to set BCs and coarse/fine values - // - // These methods are wrappers of the linop methods of the same name - // However, use of these is preferred to make sure operations - // are performed in the correct order - // - void setDomainBC (const Array& lobc, - const Array& hibc); - - void setLevelBC (int amrlev, const MultiFab* levelbcdata); - - void setCoarseFineBC (const MultiFab* crse, int crse_ratio) - { m_linop->setCoarseFineBC(crse, crse_ratio);} - - // - // Methods to perform projection - // - void project (const Vector& phi_in, Real reltol, Real atol); - void project (Real reltol, Real atol); - - // - // Get Fluxes. DO NOT USE LinOp to get fluxes!!! - // - void getFluxes (const Vector >& a_flux, - const Vector& a_sol, MLMG::Location a_loc) const; - - // - // Setters and getters - // - void setVerbose (int v) noexcept - { m_verbose = v; - m_mlmg->setVerbose(m_verbose); } - - // Methods to get underlying objects - // Use these to modify properties of MLMG and linear operator - MLLinOp& getLinOp () noexcept { return *m_linop; } - MLMG& getMLMG () noexcept { return *m_mlmg; } - - bool needInitialization() const noexcept { return m_needs_init; } - -private: - void setOptions (); - - void averageDownVelocity (); - - std::unique_ptr m_poisson; - std::unique_ptr m_abeclap; -#ifdef AMREX_USE_EB - std::unique_ptr m_eb_abeclap; - Vector m_eb_factory; -#endif - MLLinOp* m_linop = nullptr; - - Real m_const_beta = 0.; - - std::unique_ptr m_mlmg; - - Vector > m_umac; - Vector m_rhs; - Vector m_phi; - Vector m_divu; - Vector > m_fluxes; - - Vector m_geom; - - int m_verbose = 0; - - bool m_needs_domain_bcs = true; - Vector m_needs_level_bcs; - - // Location of umac -- face center vs face centroid - MLMG::Location m_umac_loc; - - MLMG::Location m_beta_loc; - - MLMG::Location m_phi_loc; - - // Location of divu (RHS -- optional) -- cell center vs cell centroid - MLMG::Location m_divu_loc; - - bool m_needs_init = true; -}; - -} - -#endif diff --git a/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp b/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp deleted file mode 100644 index fd66aff5f49..00000000000 --- a/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp +++ /dev/null @@ -1,513 +0,0 @@ -#ifdef AMREX_USE_EB -#include -#endif - -#include -#include -#include - -namespace amrex { - -MacProjector::MacProjector( - const Vector& a_geom, - MLMG::Location a_umac_loc, - MLMG::Location a_beta_loc, - MLMG::Location a_phi_loc, - MLMG::Location a_divu_loc) - : m_geom(a_geom), - m_needs_level_bcs(a_geom.size(),true), - m_umac_loc(a_umac_loc), - m_beta_loc(a_beta_loc), - m_phi_loc(a_phi_loc), - m_divu_loc(a_divu_loc) -{ - amrex::ignore_unused(m_divu_loc, m_beta_loc, m_phi_loc, m_umac_loc); -} - -MacProjector::MacProjector (const Vector >& a_umac, - MLMG::Location a_umac_loc, - const Vector >& a_beta, - MLMG::Location a_beta_loc, - MLMG::Location a_phi_loc, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const Vector& a_divu, - MLMG::Location a_divu_loc, - const Vector& a_overset_mask) - : m_umac(a_umac), - m_geom(a_geom), - m_needs_level_bcs(a_geom.size(),true), - m_umac_loc(a_umac_loc), - m_beta_loc(a_beta_loc), - m_phi_loc(a_phi_loc), - m_divu_loc(a_divu_loc) -{ - amrex::ignore_unused(m_divu_loc, m_beta_loc, m_phi_loc, m_umac_loc); - initProjector(a_lpinfo, a_beta, a_overset_mask); - setDivU(a_divu); -} - -void MacProjector::initProjector ( - const LPInfo& a_lpinfo, - const Vector >& a_beta, - const Vector& a_overset_mask) -{ - const int nlevs = a_beta.size(); - Vector ba(nlevs); - Vector dm(nlevs); - for (int ilev = 0; ilev < nlevs; ++ilev) { - ba[ilev] = amrex::convert( - a_beta[ilev][0]->boxArray(), IntVect::TheZeroVector()); - dm[ilev] = a_beta[ilev][0]->DistributionMap(); - } - - m_rhs.resize(nlevs); - m_phi.resize(nlevs); - m_fluxes.resize(nlevs); - m_divu.resize(nlevs); - -#ifdef AMREX_USE_EB - bool has_eb = a_beta[0][0]->hasEBFabFactory(); - if (has_eb) { - m_eb_factory.resize(nlevs, nullptr); - for (int ilev = 0; ilev < nlevs; ++ilev) { - m_eb_factory[ilev] = dynamic_cast( - &(a_beta[ilev][0]->Factory())); - m_rhs[ilev].define( - ba[ilev], dm[ilev], 1, 0, MFInfo(), a_beta[ilev][0]->Factory()); - m_phi[ilev].define( - ba[ilev], dm[ilev], 1, 1, MFInfo(), a_beta[ilev][0]->Factory()); - m_rhs[ilev].setVal(0.0); - m_phi[ilev].setVal(0.0); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - m_fluxes[ilev][idim].define( - amrex::convert(ba[ilev], IntVect::TheDimensionVector(idim)), - dm[ilev], 1, 0, MFInfo(), a_beta[ilev][0]->Factory()); - } - } - - m_eb_abeclap = std::make_unique(m_geom, ba, dm, a_lpinfo, m_eb_factory); - m_linop = m_eb_abeclap.get(); - - if (m_phi_loc == MLMG::Location::CellCentroid) - m_eb_abeclap->setPhiOnCentroid(); - - m_eb_abeclap->setScalars(0.0, 1.0); - for (int ilev = 0; ilev < nlevs; ++ilev) { - m_eb_abeclap->setBCoeffs(ilev, a_beta[ilev], m_beta_loc); - } - } else -#endif - { - for (int ilev = 0; ilev < nlevs; ++ilev) { - m_rhs[ilev].define(ba[ilev], dm[ilev], 1, 0); - m_phi[ilev].define(ba[ilev], dm[ilev], 1, 1); - m_rhs[ilev].setVal(0.0); - m_phi[ilev].setVal(0.0); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - m_fluxes[ilev][idim].define( - amrex::convert(ba[ilev], IntVect::TheDimensionVector(idim)), - dm[ilev], 1, 0); - } - } - - if (a_overset_mask.empty()) { - m_abeclap = std::make_unique(m_geom, ba, dm, a_lpinfo); - } else { - m_abeclap = std::make_unique(m_geom, ba, dm, a_overset_mask, a_lpinfo); - } - - m_linop = m_abeclap.get(); - - m_abeclap->setScalars(0.0, 1.0); - for (int ilev = 0; ilev < nlevs; ++ilev) { - m_abeclap->setBCoeffs(ilev, a_beta[ilev]); - } - } - - m_mlmg = std::make_unique(*m_linop); - - setOptions(); - - m_needs_init = false; -} - -void MacProjector::updateBeta ( - const Vector>& a_beta) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_linop != nullptr, - "MacProjector::updateBeta: initProjector must be called before calling this method"); - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_poisson == nullptr, - "MacProjector::updateBeta: should not be called for constant beta"); - - const int nlevs = a_beta.size(); -#ifdef AMREX_USE_EB - const bool has_eb = a_beta[0][0]->hasEBFabFactory(); - if (has_eb) { - for (int ilev=0; ilev < nlevs; ++ilev) - m_eb_abeclap->setBCoeffs(ilev, a_beta[ilev], m_beta_loc); - } else -#endif - { - for (int ilev=0; ilev < nlevs; ++ilev) - m_abeclap->setBCoeffs(ilev, a_beta[ilev]); - } -} - -void MacProjector::setUMAC( - const Vector>& a_umac) -{ - m_umac = a_umac; -} - -void MacProjector::setDivU(const Vector& a_divu) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_linop != nullptr, - "MacProjector::setDivU: initProjector must be called before calling this method"); - - for (int ilev = 0, N = a_divu.size(); ilev < N; ++ilev) { - if (a_divu[ilev]) { - if (!m_divu[ilev].ok()) { -#ifdef AMREX_USE_EB - m_divu[ilev].define( - a_divu[ilev]->boxArray(), - a_divu[ilev]->DistributionMap(), - 1,0,MFInfo(), a_divu[ilev]->Factory()); -#else - m_divu[ilev].define( - a_divu[ilev]->boxArray(), - a_divu[ilev]->DistributionMap(), 1, 0); -#endif - } - MultiFab::Copy(m_divu[ilev], *a_divu[ilev], 0, 0, 1, 0); - } - } -} - -void -MacProjector::setDomainBC (const Array& lobc, - const Array& hibc) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_linop != nullptr, - "MacProjector::setDomainBC: initProjector must be called before calling this method"); - m_linop->setDomainBC(lobc, hibc); - m_needs_domain_bcs = false; -} - - -void -MacProjector::setLevelBC (int amrlev, const MultiFab* levelbcdata) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_needs_domain_bcs, - "setDomainBC must be called before setLevelBC"); - m_linop->setLevelBC(amrlev, levelbcdata); - m_needs_level_bcs[amrlev] = false; -} - - - -void -MacProjector::project (Real reltol, Real atol) -{ - const int nlevs = m_rhs.size(); - - for (int ilev = 0; ilev < nlevs; ++ilev) { - if (m_needs_level_bcs[ilev]) { - m_linop->setLevelBC(ilev, nullptr); - m_needs_level_bcs[ilev] = false; - } - } - - if ( m_umac[0][0] ) - averageDownVelocity(); - - for (int ilev = 0; ilev < nlevs; ++ilev) - { - if ( m_umac[0][0] ) - { - Array u; - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - u[idim] = m_umac[ilev][idim]; - } -#ifdef AMREX_USE_EB - if (m_umac_loc != MLMG::Location::FaceCentroid) - { - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_umac[ilev][idim]->nGrow() > 0, - "MacProjector: with EB, umac must have at least one ghost cell if not already_on_centroid"); - m_umac[ilev][idim]->FillBoundary(m_geom[ilev].periodicity()); - } - } - - EB_computeDivergence(m_rhs[ilev], u, m_geom[ilev], (m_umac_loc == MLMG::Location::FaceCentroid)); -#else - computeDivergence(m_rhs[ilev], u, m_geom[ilev]); -#endif - - // For mlabeclaplacian, we solve -del dot (beta grad phi) = rhs - // and set up RHS as (m_divu - divu), where m_divu is a user-provided source term - // For mlpoisson, we solve `del dot grad phi = rhs/(-const_beta)` - // and set up RHS as (m_divu - divu)*(-1/const_beta) - AMREX_ASSERT(m_poisson == nullptr || m_const_beta != Real(0.0)); - m_rhs[ilev].mult(m_poisson ? Real(1.0)/m_const_beta : Real(-1.0)); - } - //else m_rhs already initialized to 0 - - if (m_divu[ilev].ok()) - { - MultiFab::Saxpy(m_rhs[ilev], m_poisson ? Real(-1.0)/m_const_beta : Real(1.0), - m_divu[ilev], 0, 0, 1, 0); - } - - // Always reset initial phi to be zero. This is needed to handle the - // situation where the MacProjector is being reused. - m_phi[ilev].setVal(0.0); - } - - m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), reltol, atol); - - if ( m_umac[0][0] ) - { - m_mlmg->getFluxes(amrex::GetVecOfArrOfPtrs(m_fluxes), m_umac_loc); - - for (int ilev = 0; ilev < nlevs; ++ilev) { - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (m_poisson) { - MultiFab::Saxpy(*m_umac[ilev][idim], m_const_beta, m_fluxes[ilev][idim], 0,0,1,0); - } else { - MultiFab::Add(*m_umac[ilev][idim], m_fluxes[ilev][idim], 0, 0, 1, 0); - } -#ifdef AMREX_USE_EB - EB_set_covered_faces(m_umac[ilev], 0.0); -#endif - } - } - - averageDownVelocity(); - } - -} - -void -MacProjector::project (const Vector& phi_inout, Real reltol, Real atol) -{ - const int nlevs = m_rhs.size(); - for (int ilev = 0; ilev < nlevs; ++ilev) { - MultiFab::Copy(m_phi[ilev], *phi_inout[ilev], 0, 0, 1, 0); - } - - project(reltol, atol); - - for (int ilev = 0; ilev < nlevs; ++ilev) { - MultiFab::Copy(*phi_inout[ilev], m_phi[ilev], 0, 0, 1, 0); - } -} - -void -MacProjector::getFluxes (const Vector >& a_flux, - const Vector& a_sol, MLMG::Location a_loc) const -{ - int ilev = 0; - if (m_needs_level_bcs[ilev]) - m_linop->setLevelBC(ilev, nullptr); - - m_linop->getFluxes(a_flux, a_sol, a_loc); - if (m_poisson) { - for (auto const& mfarr : a_flux) { - for (auto const& mfp : mfarr) { - mfp->mult(m_const_beta); - } - } - } -} - -// -// Set options by using default values and values read in input file -// -void -MacProjector::setOptions () -{ - // Default values - int maxorder(3); - int bottom_verbose(0); - int maxiter(200); - int bottom_maxiter(200); - Real bottom_rtol(1.0e-4_rt); - Real bottom_atol(-1.0_rt); - std::string bottom_solver("bicg"); - - int num_pre_smooth(2); - int num_post_smooth(2); - - // Read from input file - ParmParse pp("mac_proj"); - pp.query( "verbose" , m_verbose ); - pp.query( "maxorder" , maxorder ); - pp.query( "bottom_verbose", bottom_verbose ); - pp.query( "maxiter" , maxiter ); - pp.query( "bottom_maxiter", bottom_maxiter ); - pp.query( "bottom_rtol" , bottom_rtol ); - pp.query( "bottom_atol" , bottom_atol ); - pp.query( "bottom_solver" , bottom_solver ); - - pp.query( "num_pre_smooth" , num_pre_smooth ); - pp.query( "num_post_smooth" , num_post_smooth ); - - // Set default/input values - m_linop->setMaxOrder(maxorder); - m_mlmg->setVerbose(m_verbose); - m_mlmg->setBottomVerbose(bottom_verbose); - m_mlmg->setMaxIter(maxiter); - m_mlmg->setBottomMaxIter(bottom_maxiter); - m_mlmg->setBottomTolerance(bottom_rtol); - m_mlmg->setBottomToleranceAbs(bottom_atol); - - m_mlmg->setPreSmooth(num_pre_smooth); - m_mlmg->setPostSmooth(num_post_smooth); - - if (bottom_solver == "smoother") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::smoother); - } - else if (bottom_solver == "bicg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::bicgstab); - } - else if (bottom_solver == "cg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::cg); - } - else if (bottom_solver == "bicgcg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::bicgcg); - } - else if (bottom_solver == "cgbicg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::cgbicg); - } - else if (bottom_solver == "hypre") - { -#ifdef AMREX_USE_HYPRE - m_mlmg->setBottomSolver(MLMG::BottomSolver::hypre); -#else - amrex::Abort("AMReX was not built with HYPRE support"); -#endif - } -} - -void -MacProjector::averageDownVelocity () -{ - int finest_level = m_umac.size() - 1; - - - for (int lev = finest_level; lev > 0; --lev) - { - - IntVect rr = m_geom[lev].Domain().size() / m_geom[lev-1].Domain().size(); - -#ifdef AMREX_USE_EB - EB_average_down_faces(GetArrOfConstPtrs(m_umac[lev]), - m_umac[lev-1], - rr, m_geom[lev-1]); -#else - average_down_faces(GetArrOfConstPtrs(m_umac[lev]), - m_umac[lev-1], - rr, m_geom[lev-1]); -#endif - } -} - -#ifndef AMREX_USE_EB -void MacProjector::initProjector (Vector const& a_grids, - Vector const& a_dmap, - const LPInfo& a_lpinfo, Real const a_const_beta, - const Vector& a_overset_mask) -{ - m_const_beta = a_const_beta; - - const int nlevs = a_grids.size(); - Vector ba(nlevs); - for (int ilev = 0; ilev < nlevs; ++ilev) { - ba[ilev] = amrex::convert(a_grids[ilev], IntVect::TheZeroVector()); - } - auto const& dm = a_dmap; - - m_rhs.resize(nlevs); - m_phi.resize(nlevs); - m_fluxes.resize(nlevs); - m_divu.resize(nlevs); - - for (int ilev = 0; ilev < nlevs; ++ilev) { - m_rhs[ilev].define(ba[ilev], dm[ilev], 1, 0); - m_phi[ilev].define(ba[ilev], dm[ilev], 1, 1); - m_rhs[ilev].setVal(0.0); - m_phi[ilev].setVal(0.0); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - m_fluxes[ilev][idim].define( - amrex::convert(ba[ilev], IntVect::TheDimensionVector(idim)), - dm[ilev], 1, 0); - } - } - - if (a_overset_mask.empty()) { - m_poisson = std::make_unique(m_geom, ba, dm, a_lpinfo); - } else { - m_poisson = std::make_unique(m_geom, ba, dm, a_overset_mask, a_lpinfo); - } - - m_linop = m_poisson.get(); - - m_mlmg = std::make_unique(*m_linop); - - setOptions(); - - m_needs_init = false; -} - -MacProjector::MacProjector (const Vector >& a_umac, - const Real a_const_beta, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const Vector& a_overset_mask, - const Vector& a_divu) - : m_const_beta(a_const_beta), - m_umac(a_umac), - m_geom(a_geom), - m_umac_loc(MLMG::Location::FaceCenter), - m_beta_loc(MLMG::Location::FaceCenter), - m_phi_loc(MLMG::Location::CellCenter), - m_divu_loc(MLMG::Location::CellCenter) -{ - const int nlevs = a_umac.size(); - Vector ba(nlevs); - Vector dm(nlevs); - for (int ilev = 0; ilev < nlevs; ++ilev) { - ba[ilev] = a_umac[ilev][0]->boxArray(); - dm[ilev] = a_umac[ilev][0]->DistributionMap(); - } - initProjector(ba, dm, a_lpinfo, a_const_beta, a_overset_mask); - setDivU(a_divu); -} - -void MacProjector::updateBeta (Real a_const_beta) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_linop != nullptr, - "MacProjector::updateBeta: initProjector must be called before calling this method"); - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_poisson != nullptr, - "MacProjector::updateBeta: should not be called for variable beta"); - - m_const_beta = a_const_beta; -} - -#endif - -} diff --git a/Src/LinearSolvers/Projections/AMReX_NodalProjector.H b/Src/LinearSolvers/Projections/AMReX_NodalProjector.H deleted file mode 100644 index 6fb4f3ff8c3..00000000000 --- a/Src/LinearSolvers/Projections/AMReX_NodalProjector.H +++ /dev/null @@ -1,170 +0,0 @@ -#ifndef NODAL_PROJECTION_H -#define NODAL_PROJECTION_H -#include - -#include -#include -#include -#include -#include - -// -// -// *************************** DEFAULT MODE *************************** -// -// Solves -// -// div(sigma*grad(phi)) = div(vel) + S_nd + S_cc -// -// and performs the projection -// -// vel = vel - sigma * grad(phi) -// -// vel, sigma, and S_cc are cell-centered variables, while -// phi and S_nd are nodal-centered variables. -// -// *************************** CUSTOM MODE *************************** -// -// Solves -// -// div(sigma*grad(phi)) = rhs -// -// and performs the projection -// -// vel = vel - (sigma/alpha) * grad(phi) -// -// alpha is a cell-centered variable, while rhs is nodal-centered. -// -// In this mode, the user provides rhs and alpha -// -// By default alpha is assumed to be 1. Use setAlpha to change the default. -// -// Use setCustomRHS to provide a custom RHS, else "div(vel) + S_nd + S_cc" -// is used. -// -// Example: rhs = div(alpha*vel) -// -namespace amrex { - -class NodalProjector -{ - -public: - - NodalProjector ( const amrex::Vector& a_vel, - const amrex::Vector& a_sigma, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const amrex::Vector& a_S_cc = {}, - const amrex::Vector& a_S_nd = {} ); - - NodalProjector ( const amrex::Vector& a_vel, - const amrex::Vector& a_sigma, - const Vector& a_geom, - const amrex::Vector& a_S_cc = {}, - const amrex::Vector& a_S_nd = {} ) - : NodalProjector(a_vel, a_sigma, a_geom, LPInfo(), a_S_cc, a_S_nd ) {} - - - NodalProjector ( const amrex::Vector& a_vel, - const amrex::Real a_const_sigma, - const Vector& a_geom, - const LPInfo& a_lpinfo, - const amrex::Vector& a_S_cc = {}, - const amrex::Vector& a_S_nd = {} ); - - void project ( amrex::Real a_rtol = Real(1.0e-11), amrex::Real a_atol = Real(1.0e-14) ); - void project ( const Vector& a_phi, amrex::Real a_rtol = Real(1.0e-11), - amrex::Real a_atol = Real(1.0e-14) ); - - Vector< MultiFab* > getGradPhi () {return GetVecOfPtrs(m_fluxes);} - Vector< const MultiFab* > getGradPhiConst () const {return GetVecOfConstPtrs(m_fluxes);} - Vector< MultiFab* > getPhi () {return GetVecOfPtrs(m_phi);} - Vector< const MultiFab* > getPhiConst () const {return GetVecOfConstPtrs(m_phi);} - - void computeRHS ( const amrex::Vector& a_rhs, - const amrex::Vector& a_vel, - const amrex::Vector& a_S_cc = {}, - const amrex::Vector& a_S_nd = {} ); - - void setAlpha (const amrex::Vector a_alpha) - {m_alpha=a_alpha;m_has_alpha=true;} - void setCustomRHS (const amrex::Vector a_rhs); - - - // Methods to set verbosity - void setVerbose (int v) noexcept { m_verbose = v; } - - - // Set domain BC - void setDomainBC ( std::array a_bc_lo, - std::array a_bc_hi ); - - // Methods to get underlying objects - // Use these to modify properties of MLMG and linear operator - MLNodeLaplacian& getLinOp () noexcept { return *m_linop; } - MLMG& getMLMG () noexcept { return *m_mlmg; } - - // Methods to set MF for sync - void setSyncResidualFine (MultiFab* a_sync_resid_fine) {m_sync_resid_fine=a_sync_resid_fine;} - void setSyncResidualCrse (MultiFab* a_sync_resid_crse, IntVect a_ref_ratio, BoxArray a_fine_grids ) - {m_sync_resid_crse=a_sync_resid_crse; m_ref_ratio=a_ref_ratio; m_fine_grids=a_fine_grids;} - -private: - - void setOptions (); - void setCoarseBoundaryVelocityForSync (); - void computeSyncResidual (); - void averageDown (const amrex::Vector a_var); - void define (LPInfo const& a_lpinfo); - - bool m_has_rhs = false; - bool m_has_alpha = false; - bool m_need_bcs = true; - - // Verbosity - int m_verbose = 0; - - // Geometry - Vector m_geom; - - // EB factory if any -#ifdef AMREX_USE_EB - Vector m_ebfactory; -#endif - - // Cell-centered data - Vector m_vel; - Vector m_fluxes; - Vector m_alpha; - Vector m_S_cc; - Vector m_sigma; - amrex::Real m_const_sigma = 0.0; - - // Node-centered data - Vector m_phi; - Vector m_rhs; - Vector m_S_nd; - - // Linear operator - std::unique_ptr< MLNodeLaplacian > m_linop; - - // Solver - std::unique_ptr< MLMG > m_mlmg; - - // Boundary conditions - std::array m_bc_lo; - std::array m_bc_hi; - - // Members for synchronization - IntVect m_ref_ratio; - BoxArray m_fine_grids; // Grid at crse level + 1 - MultiFab* m_sync_resid_crse = nullptr; - MultiFab* m_sync_resid_fine = nullptr; - - void printInfo (); -}; - -} - -#endif diff --git a/Src/LinearSolvers/Projections/AMReX_NodalProjector.cpp b/Src/LinearSolvers/Projections/AMReX_NodalProjector.cpp deleted file mode 100644 index f8fb938b5e1..00000000000 --- a/Src/LinearSolvers/Projections/AMReX_NodalProjector.cpp +++ /dev/null @@ -1,524 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -namespace amrex { - - -// -// Add defaults for boundary conditions??? -// -NodalProjector::NodalProjector ( const amrex::Vector& a_vel, - const amrex::Vector& a_sigma, - const amrex::Vector& a_geom, - const LPInfo& a_lpinfo, - const amrex::Vector& a_S_cc, - const amrex::Vector& a_S_nd ) - : m_geom(a_geom), - m_vel(a_vel), - m_S_cc(a_S_cc), - m_sigma(a_sigma), - m_S_nd(a_S_nd) -{ - define(a_lpinfo); -} - -NodalProjector::NodalProjector ( const amrex::Vector& a_vel, - const amrex::Real a_const_sigma, - const amrex::Vector& a_geom, - const LPInfo& a_lpinfo, - const amrex::Vector& a_S_cc, - const amrex::Vector& a_S_nd ) - : m_geom(a_geom), - m_vel(a_vel), - m_S_cc(a_S_cc), - m_const_sigma(a_const_sigma), - m_S_nd(a_S_nd) -{ - define(a_lpinfo); -} - -void NodalProjector::define (LPInfo const& a_lpinfo) -{ - int nlevs = m_vel.size(); - - Vector ba(nlevs); - Vector dm(nlevs); - for (int lev = 0; lev < nlevs; ++lev) - { - ba[lev] = m_vel[lev]->boxArray(); - dm[lev] = m_vel[lev]->DistributionMap(); - } - - // Resize member data - m_phi.resize(nlevs); - m_fluxes.resize(nlevs); - m_rhs.resize(nlevs); - - -#ifdef AMREX_USE_EB - bool has_eb = m_vel[0] -> hasEBFabFactory(); - if (has_eb) - { - m_ebfactory.resize(nlevs,nullptr); - for (int lev = 0; lev < nlevs; ++lev ) - { - m_ebfactory[lev] = dynamic_cast(&(m_vel[lev]->Factory())); - - // Cell-centered data - m_fluxes[lev].define(ba[lev], dm[lev], AMREX_SPACEDIM, 0, MFInfo(), m_vel[lev]->Factory()); - - // Node-centered data - auto tmp = ba[lev]; - const auto& ba_nd = tmp.surroundingNodes(); - m_phi[lev].define(ba_nd, dm[lev], 1, 1, MFInfo(), m_vel[lev]->Factory()); - m_rhs[lev].define(ba_nd, dm[lev], 1, 0, MFInfo(), m_vel[lev]->Factory()); - } - } - else -#endif - { - for (int lev = 0; lev < nlevs; ++lev ) - { - // Cell-centered data - m_fluxes[lev].define(ba[lev], dm[lev], AMREX_SPACEDIM, 0); - - // Node-centered data - BoxArray tmp = ba[lev]; - const auto& ba_nd = tmp.surroundingNodes(); - m_phi[lev].define(ba_nd, dm[lev], 1, 1); - m_rhs[lev].define(ba_nd, dm[lev], 1, 0); - } - } - - // Initialize all variables - for (int lev(0); lev < m_phi.size(); ++lev) - { - m_phi[lev].setVal(0.0); - m_fluxes[lev].setVal(0.0); - m_rhs[lev].setVal(0.0); - } - - // - // Setup linear operator - // - if (m_sigma.empty()) { -#ifdef AMREX_USE_EB - m_linop = std::make_unique(m_geom, ba, dm, a_lpinfo, m_ebfactory, - m_const_sigma); -#else - m_linop = std::make_unique(m_geom, ba, dm, a_lpinfo, - Vector const*>{}, - m_const_sigma); -#endif - } else { -#ifdef AMREX_USE_EB - m_linop = std::make_unique(m_geom, ba, dm, a_lpinfo, m_ebfactory); -#else - m_linop = std::make_unique(m_geom, ba, dm, a_lpinfo); -#endif - } - - m_linop->setGaussSeidel(true); - m_linop->setHarmonicAverage(false); - - // - // Setup solver - // - m_mlmg = std::make_unique(*m_linop); - - setOptions(); -} - - -// -// Set options by using default values and values read in input file -// -void -NodalProjector::setOptions () -{ - // Default values - int bottom_verbose(0); - int maxiter(100); - int bottom_maxiter(100); - Real bottom_rtol(1.0e-4_rt); - Real bottom_atol(-1.0_rt); - std::string bottom_solver("bicgcg"); - - int num_pre_smooth (2); - int num_post_smooth(2); - - Real normalization_threshold(-1.); - - // Read from input file - ParmParse pp("nodal_proj"); - pp.query( "verbose" , m_verbose ); - pp.query( "bottom_verbose", bottom_verbose ); - pp.query( "maxiter" , maxiter ); - pp.query( "bottom_maxiter", bottom_maxiter ); - pp.query( "bottom_rtol" , bottom_rtol ); - pp.query( "bottom_atol" , bottom_atol ); - pp.query( "bottom_solver" , bottom_solver ); - - pp.query( "normalization_threshold" , normalization_threshold); - - pp.query( "num_pre_smooth" , num_pre_smooth ); - pp.query( "num_post_smooth" , num_post_smooth ); - - // This is only used by the Krylov solvers but we pass it through the nodal operator - // if it is set here. Otherwise we use the default set in AMReX_NodeLaplacian.H - if (normalization_threshold > 0.) - m_linop->setNormalizationThreshold(normalization_threshold); - - // Set default/input values - m_mlmg->setVerbose(m_verbose); - m_mlmg->setBottomVerbose(bottom_verbose); - m_mlmg->setMaxIter(maxiter); - m_mlmg->setBottomMaxIter(bottom_maxiter); - m_mlmg->setBottomTolerance(bottom_rtol); - m_mlmg->setBottomToleranceAbs(bottom_atol); - - m_mlmg->setPreSmooth(num_pre_smooth); - m_mlmg->setPostSmooth(num_post_smooth); - - if (bottom_solver == "smoother") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::smoother); - } - else if (bottom_solver == "bicg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::bicgstab); - } - else if (bottom_solver == "cg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::cg); - } - else if (bottom_solver == "bicgcg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::bicgcg); - } - else if (bottom_solver == "cgbicg") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::cgbicg); - } -#ifdef AMREX_USE_HYPRE - else if (bottom_solver == "hypre") - { - m_mlmg->setBottomSolver(MLMG::BottomSolver::hypre); - } -#endif -} - -void -NodalProjector::setDomainBC ( std::array a_bc_lo, - std::array a_bc_hi ) -{ - m_bc_lo=a_bc_lo; - m_bc_hi=a_bc_hi; - m_linop->setDomainBC(m_bc_lo,m_bc_hi); - m_need_bcs = false; -} - -void -NodalProjector::setCustomRHS (const amrex::Vector a_rhs) -{ - AMREX_ALWAYS_ASSERT(m_rhs.size()==a_rhs.size()); - - for (int lev=0; lev < m_rhs.size(); ++lev) - { - MultiFab::Copy(m_rhs[lev], *a_rhs[lev], 0, 0, 1, 0); - } - - m_has_rhs = true; -} - - -void -NodalProjector::project ( Real a_rtol, Real a_atol ) -{ - BL_PROFILE("NodalProjector::project"); - AMREX_ALWAYS_ASSERT(!m_need_bcs); - - if (m_verbose > 0) - amrex::Print() << "Nodal Projection:" << std::endl; - - // - // Average fine grid velocity values down to the coarse grid - // By doing this operation at this time we: - // - // 1) fill regions covered by a finer grid with some "valid" data, i.e. not NaNs - // for example. This is required internally by MLMG. - // 2) make the velocity field ready for projection on the entirety of each level. - // - averageDown(m_vel); - - // Set matrix coefficients - for (int lev = 0; lev < m_sigma.size(); ++lev) - { - m_linop -> setSigma(lev, *m_sigma[lev]); - } - - // Compute RHS if necessary - if (!m_has_rhs) - { - computeRHS( GetVecOfPtrs(m_rhs), m_vel, m_S_cc, m_S_nd ); - } - - // Print diagnostics - if (m_verbose > 0) - { - amrex::Print() << " >> Before projection:" << std::endl; - printInfo(); - amrex::Print() << std::endl; - } - - // Solve - // phi comes out already averaged-down and ready to be used by caller if needed - m_mlmg -> solve( GetVecOfPtrs(m_phi), GetVecOfConstPtrs(m_rhs), a_rtol, a_atol ); - - // Get fluxes -- fluxes = - sigma * grad(phi) - m_mlmg -> getFluxes( GetVecOfPtrs(m_fluxes) ); - - // Compute sync residual BEFORE performing projection - computeSyncResidual(); - - // Perform projection - for (int lev(0); lev < m_phi.size(); ++lev) - { - if (m_has_alpha) - { - // fluxes -> fluxes/alpha = - ( sigma / alpha ) * grad(phi) - for (int n = 0; n < AMREX_SPACEDIM; ++n) - { - MultiFab::Divide( m_fluxes[lev], *m_alpha[lev], 0, n, 1, 0 ); - } - } - - // - // vel = vel + fluxes = vel - ( sigma / alpha ) * grad(phi) - // - MultiFab::Add( *m_vel[lev], m_fluxes[lev], 0, 0, AMREX_SPACEDIM, 0); - - // set m_fluxes = grad(phi) - m_linop->compGrad(lev,m_fluxes[lev],m_phi[lev]); - } - - // - // At this time, results are "correct" only on regions not covered by finer grids. - // We average them down so that they are "correct" everywhere in each level. - // - averageDown(GetVecOfPtrs(m_fluxes)); - averageDown(m_vel); - - - // Print diagnostics - if ( (m_verbose > 0) && (!m_has_rhs)) - { - computeRHS( GetVecOfPtrs(m_rhs), m_vel, m_S_cc, m_S_nd ); - amrex::Print() << " >> After projection:" << std::endl; - printInfo(); - amrex::Print() << std::endl; - } - -} - -void -NodalProjector::project ( const Vector& a_phi, Real a_rtol, Real a_atol ) -{ - AMREX_ALWAYS_ASSERT(a_phi.size()==m_phi.size()); - - for (int lev=0; lev < m_phi.size(); ++lev ) - { - MultiFab::Copy(m_phi[lev],*a_phi[lev],0,0,1,m_phi[lev].nGrow()); - } - - project(a_rtol, a_atol); - - for (int lev=0; lev < m_phi.size(); ++lev ) - { - MultiFab::Copy(*a_phi[lev],m_phi[lev],0,0,1,m_phi[lev].nGrow()); - } -} - - -// -// Compute RHS: div(u) + S_nd + S_cc -// -void -NodalProjector::computeRHS ( const amrex::Vector& a_rhs, - const amrex::Vector& a_vel, - const amrex::Vector& a_S_cc, - const amrex::Vector& a_S_nd ) -{ - AMREX_ALWAYS_ASSERT(!m_need_bcs); // This is needed to use linop - AMREX_ALWAYS_ASSERT(a_rhs.size()==m_phi.size()); - AMREX_ALWAYS_ASSERT(a_vel.size()==m_phi.size()); - AMREX_ALWAYS_ASSERT((a_S_cc.size()==0) || (a_S_cc.size()==m_phi.size()) ); - AMREX_ALWAYS_ASSERT((a_S_nd.size()==0) || (a_S_nd.size()==m_phi.size()) ); - - BL_PROFILE("NodalProjector::computeRHS"); - - bool has_S_cc(!a_S_cc.empty()); - bool has_S_nd(!a_S_nd.empty()); - - // Check the type of grids used - for (int lev(0); lev < m_phi.size(); ++lev) - { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_rhs[lev]->ixType().nodeCentered(), "a_rhs is not node centered"); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_vel[lev]->ixType().cellCentered(), "a_vel is not cell centered"); - - if (has_S_cc && (a_S_cc[lev]!=nullptr) ) - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_S_cc[lev]->ixType().cellCentered(),"a_S_cc is not cell centered"); - - if (has_S_nd && (a_S_nd[lev]!=nullptr) ) - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_S_nd[lev]->ixType().nodeCentered(),"a_S_nd is not node centered"); - } - - m_linop -> compRHS( a_rhs, a_vel, a_S_nd, a_S_cc ); -} - - -void -NodalProjector::printInfo () -{ - for (int lev(0); lev < m_rhs.size(); ++lev) - { - amrex::Print() << " * On lev " << lev - << " max(abs(rhs)) = " - << m_rhs[lev].norm0(0,0,false,true) - << std::endl; - } -} - - -void -NodalProjector::computeSyncResidual () -{ - BL_PROFILE("NodalProjector::computeSyncResidual"); - - if ( (m_sync_resid_fine != nullptr) || (m_sync_resid_crse != nullptr) ) - { - int c_lev = 0; - - setCoarseBoundaryVelocityForSync(); - - if (m_sync_resid_fine != nullptr) - { - MultiFab* rhptr = nullptr; - if (!m_S_cc.empty()) - rhptr = m_S_cc[c_lev]; - m_linop->compSyncResidualFine(*m_sync_resid_fine, m_phi[c_lev], *m_vel[c_lev], rhptr); - } - - if (m_sync_resid_crse != nullptr) - { - MultiFab* rhptr = nullptr; - if (!m_S_cc.empty()) - rhptr = m_S_cc[c_lev]; - - // this requires sigma to have 2 ghost cells (valid at -2) - m_linop->compSyncResidualCoarse(*m_sync_resid_crse, m_phi[c_lev], *m_vel[c_lev], - rhptr, m_fine_grids, m_ref_ratio); - } - - - } - -} - - -// Set velocity in ghost cells to zero except for inflow -// 1) At non-inflow faces, the normal component of velocity will be completely zero'd -// 2) The normal velocity at corners -- even periodic corners -- just outside inflow faces will be zero'd -void -NodalProjector::setCoarseBoundaryVelocityForSync () -{ - - const BoxArray& grids = m_vel[0]->boxArray(); - const Box& domainBox = m_geom[0].Domain(); - - for (int idir=0; idir < AMREX_SPACEDIM; idir++) - { - if (m_bc_lo[idir] != LinOpBCType::inflow && m_bc_hi[idir] != LinOpBCType::inflow) - { - m_vel[0]->setBndry(0.0, idir, 1); - } - else - { -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(*m_vel[0]); mfi.isValid(); ++mfi) - { - int i = mfi.index(); - FArrayBox& v_fab = (*m_vel[0])[mfi]; - - const Box& reg = grids[i]; - const Box& bxg1 = amrex::grow(reg,1); - BoxList bxlist(reg); - - if (m_bc_lo[idir] == LinOpBCType::inflow && reg.smallEnd(idir) == domainBox.smallEnd(idir)) - { - Box bx; // bx is the region we *protect* from zero'ing - bx = amrex::adjCellLo(reg, idir); - bxlist.push_back(bx); - } - - if (m_bc_hi[idir] == LinOpBCType::inflow && reg.bigEnd(idir) == domainBox.bigEnd(idir)) { - Box bx; // bx is the region we *protect* from zero'ing - bx = amrex::adjCellHi(reg, idir); - bxlist.push_back(bx); - } - - BoxList bxlist2 = amrex::complementIn(bxg1, bxlist); - - for (BoxList::iterator it=bxlist2.begin(); it != bxlist2.end(); ++it) - { - Box ovlp = *it & v_fab.box(); - if (ovlp.ok()) - { - v_fab.setVal(0.0, ovlp, idir, 1); - } - } - } - } - } - -} - - -void -NodalProjector::averageDown (const amrex::Vector a_var) -{ - - int f_lev = a_var.size()-1; - int c_lev = 0; - - for (int lev = f_lev-1; lev >= c_lev; --lev) - { - IntVect rr = m_geom[lev+1].Domain().size() / m_geom[lev].Domain().size(); - -#ifdef AMREX_USE_EB - const auto ebf = dynamic_cast(a_var[lev+1]->Factory()); - - amrex::MultiFab volume(a_var[lev+1]->boxArray(),a_var[lev+1]->DistributionMap(),1,0); - m_geom[lev+1].GetVolume(volume); - - EB_average_down(*a_var[lev+1], *a_var[lev], volume, ebf.getVolFrac(), - 0, a_var[lev]->nComp(), rr); -#else - average_down(*a_var[lev+1], *a_var[lev], m_geom[lev+1], m_geom[lev], - 0, a_var[lev]->nComp(), rr); -#endif - - } - - -} - - -} diff --git a/Src/LinearSolvers/Projections/Make.package b/Src/LinearSolvers/Projections/Make.package deleted file mode 100644 index 43b595ac3b0..00000000000 --- a/Src/LinearSolvers/Projections/Make.package +++ /dev/null @@ -1,8 +0,0 @@ -CEXE_headers += AMReX_MacProjector.H -CEXE_headers += AMReX_NodalProjector.H - -CEXE_sources += AMReX_MacProjector.cpp -CEXE_sources += AMReX_NodalProjector.cpp - -VPATH_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/Projections -INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/Projections diff --git a/Tests/LinearSolvers/MAC_Projection_EB/CMakeLists.txt b/Tests/LinearSolvers/MAC_Projection_EB/CMakeLists.txt deleted file mode 100644 index 429354d5b31..00000000000 --- a/Tests/LinearSolvers/MAC_Projection_EB/CMakeLists.txt +++ /dev/null @@ -1,11 +0,0 @@ -if ( (NOT AMReX_EB) OR NOT (AMReX_SPACEDIM EQUAL 3) ) - return() -endif () - -set(_sources main.cpp) -set(_input_files inputs_3d) - -setup_test(_sources _input_files) - -unset(_sources) -unset(_input_files) diff --git a/Tests/LinearSolvers/MAC_Projection_EB/GNUmakefile b/Tests/LinearSolvers/MAC_Projection_EB/GNUmakefile deleted file mode 100644 index 93bf0ab7224..00000000000 --- a/Tests/LinearSolvers/MAC_Projection_EB/GNUmakefile +++ /dev/null @@ -1,31 +0,0 @@ -USE_MPI = TRUE -USE_OMP = FALSE - -COMP = gnu - -DIM = 3 - -DEBUG = FALSE - -AMREX_HOME = ../../.. - -USE_EB = TRUE - -USE_HYPRE = TRUE -USE_HYPRE = FALSE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package - -Pdirs := AmrCore -Pdirs += Base -Pdirs += Boundary -Pdirs += EB -Pdirs += LinearSolvers/MLMG LinearSolvers/Projections - -Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) - -include $(Ppack) - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Tests/LinearSolvers/MAC_Projection_EB/Make.package b/Tests/LinearSolvers/MAC_Projection_EB/Make.package deleted file mode 100644 index 6b4b865e8fc..00000000000 --- a/Tests/LinearSolvers/MAC_Projection_EB/Make.package +++ /dev/null @@ -1 +0,0 @@ -CEXE_sources += main.cpp diff --git a/Tests/LinearSolvers/MAC_Projection_EB/README b/Tests/LinearSolvers/MAC_Projection_EB/README deleted file mode 100644 index 88bfe43f3f2..00000000000 --- a/Tests/LinearSolvers/MAC_Projection_EB/README +++ /dev/null @@ -1,66 +0,0 @@ - -This tutorial demonstrates the solution of a Poisson equation -to compute a potential flow field around nine obstacles. - -**************************************************************************************************** - -To run it in serial, - -./main3d.ex inputs - -To run it in parallel, for example on 4 ranks: - -mpirun -n 4 ./main3d.ex inputs - -The following parameters can be set at run-time -- these are currently set in the inputs -file but you can also set them on the command line. - -n_cell = 128 # number of cells in x-direction; we double this in the y-direction -max_grid_size = 64 # the maximum number of cells in any direction in a single grid - -**************************************************************************************************** - -The output from your run should look something like this: - -******************************************************************** - Let's project the initial velocity to find - the flow field around the obstacles ... - The domain has 256 cells in the x-direction - The maximum grid size is 64 -******************************************************************** - -MLMG: # of AMR levels: 1 - # of MG levels on the coarsest AMR level: 3 -MLMG: Initial rhs = 128 -MLMG: Initial residual (resid0) = 128 -MLCGSolver_BiCGStab: Initial error (error0) = 0.00585083279 -MLCGSolver_BiCGStab: Final: Iteration 131 rel. err. 8.316094303e-05 -MLMG: Iteration 1 Fine resid/bnorm = 0.5844966018 -MLCGSolver_BiCGStab: Initial error (error0) = 0.0001712929542 -MLCGSolver_BiCGStab: Final: Iteration 117 rel. err. 9.471637432e-05 -MLMG: Iteration 2 Fine resid/bnorm = 0.008993801295 -MLCGSolver_BiCGStab: Initial error (error0) = 7.062486538e-06 -MLCGSolver_BiCGStab: Final: Iteration 126 rel. err. 8.921885015e-05 -MLMG: Iteration 3 Fine resid/bnorm = 0.0003448363733 -MLCGSolver_BiCGStab: Initial error (error0) = 3.522653229e-07 -MLCGSolver_BiCGStab: Final: Iteration 131 rel. err. 6.102633163e-05 -MLMG: Iteration 4 Fine resid/bnorm = 1.818381726e-05 -MLCGSolver_BiCGStab: Initial error (error0) = 1.869194089e-08 -MLCGSolver_BiCGStab: Final: Iteration 131 rel. err. 6.35413132e-05 -MLMG: Iteration 5 Fine resid/bnorm = 1.023441528e-06 -MLCGSolver_BiCGStab: Initial error (error0) = 1.017130677e-09 -MLCGSolver_BiCGStab: Final: Iteration 93 rel. err. 9.168336176e-05 -MLMG: Iteration 6 Fine resid/bnorm = 5.871105663e-08 -MLCGSolver_BiCGStab: Initial error (error0) = 5.645789832e-11 -MLCGSolver_BiCGStab: Final: Iteration 132 rel. err. 8.280742281e-05 -MLMG: Iteration 7 Fine resid/bnorm = 3.394263815e-09 -MLMG: Final Iter. 7 resid, resid/bnorm = 4.344657683e-07, 3.394263815e-09 -MLMG: Timers: Solve = 1.15662553 Iter = 1.132159153 Bottom = 0.499076789 - -******************************************************************** - Done! -******************************************************************** - -Writing plt00000 -Total run time 1.435699991 - diff --git a/Tests/LinearSolvers/MAC_Projection_EB/inputs_3d b/Tests/LinearSolvers/MAC_Projection_EB/inputs_3d deleted file mode 100644 index b2ddd7d7fc2..00000000000 --- a/Tests/LinearSolvers/MAC_Projection_EB/inputs_3d +++ /dev/null @@ -1,13 +0,0 @@ -n_cell = 128 # number of cells in y-direction; we double this in the x-direction and divide by 8 in the z-direction -max_grid_size = 64 # the maximum number of cells in any direction in a single grid (default: 32) - -obstacles = 0 1 2 3 4 5 6 7 8 # this is how we choose which obstacles to include - -#The parameters below specify solver choices - -use_hypre = 0 # if 1 then use hypre instead of geometric multigrid (default: 0) - -mg_verbose = 2 # specify verbosity of geometric multigrid solver (default: 0) -bottom_verbose = 2 # specify verbosity of the bottom solver (default: 0) - -regtest = 0 # If regtest == 1 then we don't plot zvel (used for regression testing with GPUs) diff --git a/Tests/LinearSolvers/MAC_Projection_EB/main.cpp b/Tests/LinearSolvers/MAC_Projection_EB/main.cpp deleted file mode 100644 index aba89ddf96f..00000000000 --- a/Tests/LinearSolvers/MAC_Projection_EB/main.cpp +++ /dev/null @@ -1,264 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace amrex; - -void write_plotfile(const Geometry& geom, const MultiFab& plotmf, int regtest) -{ - std::stringstream sstream; - sstream << "plt00000"; - std::string plotfile_name = sstream.str(); - - amrex::Print() << "Writing " << plotfile_name << std::endl; - -#if (AMREX_SPACEDIM == 2) - EB_WriteSingleLevelPlotfile(plotfile_name, plotmf, - { "proc" ,"xvel", "yvel" }, - geom, 0.0, 0); -#elif (AMREX_SPACEDIM == 3) - - if (regtest == 1) - { - MultiFab small_mf(plotmf.boxArray(), plotmf.DistributionMap(), AMREX_SPACEDIM, 0, MFInfo(), plotmf.Factory()); - MultiFab::Copy(small_mf, plotmf, 0, 0, small_mf.nComp(),0); - - EB_WriteSingleLevelPlotfile(plotfile_name, small_mf, - { "proc", "xvel", "yvel" }, - geom, 0.0, 0); - } else { - EB_WriteSingleLevelPlotfile(plotfile_name, plotmf, - { "proc", "xvel", "yvel", "zvel" }, - geom, 0.0, 0); - } -#endif -} - -int main (int argc, char* argv[]) -{ - // Turn off amrex-related output - amrex::SetVerbose(0); - - amrex::Initialize(argc, argv); - - auto strt_time = amrex::second(); - - { - int mg_verbose = 0; - int bottom_verbose = 0; - int n_cell = 128; - int max_grid_size = 32; - int use_hypre = 0; - int regtest = 0; - - Real obstacle_radius = 0.10; - - // read parameters - { - ParmParse pp; - pp.query("mg_verbose", mg_verbose); - pp.query("bottom_verbose", bottom_verbose); - pp.query("n_cell", n_cell); - pp.query("max_grid_size", max_grid_size); - pp.query("use_hypre", use_hypre); - pp.query("regtest", regtest); - } - -#ifndef AMREX_USE_HYPRE - if (use_hypre == 1) - amrex::Abort("Cant use hypre if we dont build with USE_HYPRE=TRUE"); -#endif - - if (n_cell%8 != 0) - amrex::Abort("n_cell must be a multiple of 8"); - - int n_cell_y = n_cell; - int n_cell_x = 2*n_cell; - int n_cell_z = n_cell/8; - - Real ylen = 1.0; - Real xlen = 2.0 * ylen; - Real zlen = ylen / 8.0; - - Geometry geom; - BoxArray grids; - DistributionMapping dmap; - { - RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(xlen,ylen,zlen)}); - - Array isp{AMREX_D_DECL(0,1,1)}; - Geometry::Setup(&rb, 0, isp.data()); - Box domain(IntVect{AMREX_D_DECL(0,0,0)}, - IntVect{AMREX_D_DECL(n_cell_x-1,n_cell_y-1,n_cell_z-1)}); - geom.define(domain); - - grids.define(domain); - grids.maxSize(max_grid_size); - - dmap.define(grids); - } - - Array vel; - Array beta; - MultiFab plotfile_mf; - - int required_coarsening_level = 0; // typically the same as the max AMR level index - int max_coarsening_level = 100; // typically a huge number so MG coarsens as much as possible - - amrex::Vector obstacle_center = { - {AMREX_D_DECL(0.3,0.2,0.5)}, - {AMREX_D_DECL(0.3,0.5,0.5)}, - {AMREX_D_DECL(0.3,0.8,0.5)}, - {AMREX_D_DECL(0.7,0.25,0.5)}, - {AMREX_D_DECL(0.7,0.60,0.5)}, - {AMREX_D_DECL(0.7,0.85,0.5)}, - {AMREX_D_DECL(1.1,0.2,0.5)}, - {AMREX_D_DECL(1.1,0.5,0.5)}, - {AMREX_D_DECL(1.1,0.8,0.5)}}; - - int direction = 2; - Real height = -1.0; // Putting a negative number for height means it extends beyond the domain - - // The "false" below is the boolean that determines if the fluid is inside ("true") or - // outside ("false") the object(s) - - Array obstacles{ - EB2::CylinderIF( obstacle_radius, height, direction, obstacle_center[ 0], false), - EB2::CylinderIF( obstacle_radius, height, direction, obstacle_center[ 1], false), - EB2::CylinderIF( obstacle_radius, height, direction, obstacle_center[ 2], false), - EB2::CylinderIF(0.9*obstacle_radius, height, direction, obstacle_center[ 3], false), - EB2::CylinderIF(0.9*obstacle_radius, height, direction, obstacle_center[ 4], false), - EB2::CylinderIF(0.9*obstacle_radius, height, direction, obstacle_center[ 5], false), - EB2::CylinderIF( obstacle_radius, height, direction, obstacle_center[ 6], false), - EB2::CylinderIF( obstacle_radius, height, direction, obstacle_center[ 7], false), - EB2::CylinderIF( obstacle_radius, height, direction, obstacle_center[ 8], false)}; - - auto group_1 = EB2::makeUnion(obstacles[0],obstacles[1],obstacles[2]); - auto group_2 = EB2::makeUnion(obstacles[3],obstacles[4],obstacles[5]); - auto group_3 = EB2::makeUnion(obstacles[6],obstacles[7],obstacles[8]); - auto all = EB2::makeUnion(group_1,group_2,group_3); - auto gshop9 = EB2::makeShop(all); - EB2::Build(gshop9, geom, required_coarsening_level, max_coarsening_level); - - const EB2::IndexSpace& eb_is = EB2::IndexSpace::top(); - const EB2::Level& eb_level = eb_is.getLevel(geom); - - // options are basic, volume, or full - EBSupport ebs = EBSupport::full; - - // number of ghost cells for each of the 3 EBSupport types - Vector ng_ebs = {2,2,2}; - - // This object provides access to the EB database in the format of basic AMReX objects - // such as BaseFab, FArrayBox, FabArray, and MultiFab - EBFArrayBoxFactory factory(eb_level, geom, grids, dmap, ng_ebs, ebs); - - // allocate face-centered velocities and face-centered beta coefficient - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - vel[idim].define (amrex::convert(grids,IntVect::TheDimensionVector(idim)), dmap, 1, 1, - MFInfo(), factory); - beta[idim].define(amrex::convert(grids,IntVect::TheDimensionVector(idim)), dmap, 1, 0, - MFInfo(), factory); - beta[idim].setVal(1.0); // set beta to 1 - } - - // If we want to supply a non-zero S we must allocate and fill it outside the solver - // MultiFab S(grids, dmap, 1, 0, MFInfo(), factory); - // Set S here ... - - // store plotfile variables; velocity and processor id - plotfile_mf.define(grids, dmap, AMREX_SPACEDIM+1, 0, MFInfo(), factory); - - // set initial velocity to u=(1,0,0) - AMREX_D_TERM(vel[0].setVal(1.0);, - vel[1].setVal(0.0);, - vel[2].setVal(0.0);); - - LPInfo lp_info; - - // If we want to use hypre to solve the full problem we need to not coarsen inside AMReX - if (use_hypre) - lp_info.setMaxCoarseningLevel(0); - - MacProjector macproj({amrex::GetArrOfPtrs(vel)}, // mac velocity - MLMG::Location::FaceCenter, // Location of vel - {amrex::GetArrOfConstPtrs(beta)}, // beta - MLMG::Location::FaceCenter, // Location of beta - MLMG::Location::CellCenter, // Location of solution variable phi - {geom}, // the geometry object - lp_info); // structure for passing info to the operator - - // Here we specify the desired divergence S - // MacProjector macproj({amrex::GetArrOfPtrs(vel)}, // face-based velocity - // {amrex::GetArrOfConstPtrs(beta)}, // beta - // {geom}, // the geometry object - // lp_info, // structure for passing info to the operator - // {&S}); // defines the specified RHS divergence - - // Set bottom-solver to use hypre instead of native BiCGStab - if (use_hypre) - macproj.getMLMG().setBottomSolver(MLMG::BottomSolver::hypre); - - // Hard-wire the boundary conditions to be Neumann on the low x-face, Dirichlet - // on the high x-face, and periodic in the other two directions - // (the first argument is for the low end, the second is for the high end) - macproj.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann, - LinOpBCType::Periodic, - LinOpBCType::Periodic)}, - {AMREX_D_DECL(LinOpBCType::Dirichlet, - LinOpBCType::Periodic, - LinOpBCType::Periodic)}); - - macproj.setVerbose(mg_verbose); - macproj.getMLMG().setBottomVerbose(bottom_verbose); - - // Define the relative tolerance - Real reltol = 1.e-8; - - // Define the absolute tolerance; note that this argument is optional - Real abstol = 1.e-15; - - amrex::Print() << " \n********************************************************************" << std::endl; - amrex::Print() << " Let's project the initial velocity to find " << std::endl; - amrex::Print() << " the flow field around the obstacles ... " << std::endl; - amrex::Print() << " The domain has " << n_cell_x << " cells in the x-direction " << std::endl; - amrex::Print() << " The maximum grid size is " << max_grid_size << std::endl; - amrex::Print() << "******************************************************************** \n" << std::endl; - - // Solve for phi and subtract from the velocity to make it divergence-free - // Note that the normal velocities are at face centers (not centroids) - macproj.project(reltol,abstol); - - // If we want to use phi elsewhere, we can pass in an array in which to return the solution - // MultiFab phi_inout(grids, dmap, 1, 1, MFInfo(), factory); - // macproj.project_center_vels({&phi_inout},reltol,abstol,MLMG::Location::FaceCenter); - - amrex::Print() << " \n********************************************************************" << std::endl; - amrex::Print() << " Done!" << std::endl; - amrex::Print() << "******************************************************************** \n" << std::endl; - - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - vel[idim].FillBoundary(geom.periodicity()); - } - - // Copy processor id into plotfile_mf - plotfile_mf.setVal(ParallelDescriptor::MyProc(), 0, 1); - - // Copy velocity into plotfile (starting at component 1) - average_face_to_cellcenter(plotfile_mf,1,amrex::GetArrOfConstPtrs(vel)); - - write_plotfile(geom, plotfile_mf, regtest); - } - - auto stop_time = amrex::second() - strt_time; - amrex::Print() << "Total run time " << stop_time << std::endl; - - amrex::Finalize(); -} diff --git a/Tools/RegressionTesting/AMReX-cuda-tests.ini b/Tools/RegressionTesting/AMReX-cuda-tests.ini index 949ad42ee10..534b4b31364 100644 --- a/Tools/RegressionTesting/AMReX-cuda-tests.ini +++ b/Tools/RegressionTesting/AMReX-cuda-tests.ini @@ -39,7 +39,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/AMReX/ for more #MPIcommand = mpiexec -host @host@ -n @nprocs@ @command@ MPIcommand = mpiexec -n @nprocs@ @command@ -MPIhost = +MPIhost = [AMReX] dir = /home/regtester/git/amrex/ @@ -77,7 +77,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[MLMG_PoisLev] +[MLMG_PoisLev] addToCompileString = USE_CUDA=TRUE buildDir = Tests/LinearSolvers/ABecLaplacian_C inputFile = inputs-rt-poisson-lev @@ -192,23 +192,6 @@ doVis = 0 outputFile = plot testSrcTree = C_Src -[EB_MacProj_3D] -addToCompileString = USE_CUDA=TRUE -tolerance = 1.e-9 -buildDir = Tests/LinearSolvers/MAC_Projection_EB -inputFile = inputs_3d -dim = 3 -restartTest = 0 -useMPI = 1 -numprocs = 4 -useOMP = 0 -numthreads = 1 -compileTest = 0 -doVis = 0 -runtime_params = regtest=1 -outputFile = plt00000 -testSrcTree = C_Src - [EB_Node_2D] addToCompileString = USE_CUDA=TRUE tolerance = 5.e-14 diff --git a/Tools/RegressionTesting/AMReX-hip-tests.ini b/Tools/RegressionTesting/AMReX-hip-tests.ini index c8f63401bb3..3d6499910ed 100644 --- a/Tools/RegressionTesting/AMReX-hip-tests.ini +++ b/Tools/RegressionTesting/AMReX-hip-tests.ini @@ -37,7 +37,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/AMReX/ for more #MPIcommand = mpiexec -host @host@ -n @nprocs@ @command@ MPIcommand = mpiexec -n @nprocs@ @command@ -MPIhost = +MPIhost = [AMReX] dir = ../amrex/ @@ -75,7 +75,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[MLMG_PoisLev] +[MLMG_PoisLev] addToCompileString = USE_CUDA=FALSE USE_HIP=TRUE buildDir = Tests/LinearSolvers/ABecLaplacian_C inputFile = inputs-rt-poisson-lev @@ -190,23 +190,6 @@ doVis = 0 outputFile = plot testSrcTree = C_Src -[EB_MacProj_3D] -addToCompileString = USE_CUDA=FALSE USE_HIP=TRUE -tolerance = 1.e-9 -buildDir = Tests/LinearSolvers/MAC_Projection_EB -inputFile = inputs_3d -dim = 3 -restartTest = 0 -useMPI = 0 -numprocs = 4 -useOMP = 0 -numthreads = 1 -compileTest = 0 -doVis = 0 -runtime_params = regtest=1 -outputFile = plt00000 -testSrcTree = C_Src - [EB_Node_2D] addToCompileString = USE_CUDA=FALSE USE_HIP=TRUE tolerance = 5.e-14 diff --git a/Tools/RegressionTesting/AMReX-tests.ini b/Tools/RegressionTesting/AMReX-tests.ini index c6d2132ce28..c80ebd35e4a 100644 --- a/Tools/RegressionTesting/AMReX-tests.ini +++ b/Tools/RegressionTesting/AMReX-tests.ini @@ -37,7 +37,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/AMReX/ for more det #MPIcommand = mpiexec -host @host@ -n @nprocs@ @command@ MPIcommand = mpiexec -n @nprocs@ @command@ -MPIhost = +MPIhost = [AMReX] dir = /home/regtester/AMReX_RegTesting/amrex/ @@ -57,7 +57,7 @@ compileTest = 1 doVis = 0 testSrcTree = C_Src -[MLMG_PoisLev] +[MLMG_PoisLev] buildDir = Tests/LinearSolvers/ABecLaplacian_C inputFile = inputs-rt-poisson-lev dim = 3 @@ -85,7 +85,7 @@ doVis = 0 outputFile = plot testSrcTree = C_Src -[MLMG_FI_PoisCom] +[MLMG_FI_PoisCom] buildDir = Tests/LinearSolvers/ABecLaplacian_F inputFile = inputs-rt-poisson-com dim = 3 @@ -113,7 +113,7 @@ doVis = 0 outputFile = plot testSrcTree = C_Src -[AMR_Adv_C_2D] +[AMR_Adv_C_2D] buildDir = Tests/Amr/Advection_AmrLevel/Exec/UniformVelocity inputFile = inputs.regt dim = 2 @@ -126,7 +126,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[AMR_Adv_C_2D_Tracers] +[AMR_Adv_C_2D_Tracers] buildDir = Tests/Amr/Advection_AmrLevel/Exec/SingleVortex inputFile = inputs.tracers dim = 2 @@ -141,7 +141,7 @@ compareParticles = 1 particleTypes = Tracer testSrcTree = C_Src -[AMR_Adv_C_3D] +[AMR_Adv_C_3D] buildDir = Tests/Amr/Advection_AmrLevel/Exec/SingleVortex inputFile = inputs dim = 3 @@ -154,7 +154,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[AMR_Adv_C_v2_2D] +[AMR_Adv_C_v2_2D] buildDir = Tests/Amr/Advection_AmrCore/Exec/SingleVortex inputFile = inputs dim = 2 @@ -167,7 +167,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[AMR_Adv_C_v2_3D] +[AMR_Adv_C_v2_3D] buildDir = Tests/Amr/Advection_AmrCore/Exec/SingleVortex inputFile = inputs dim = 3 @@ -180,7 +180,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[AMR_Adv_CF_2D] +[AMR_Adv_CF_2D] buildDir = Tests/FortranInterface/Advection_F/Exec/SingleVortex inputFile = inputs.rt dim = 2 @@ -193,7 +193,7 @@ compileTest = 0 doVis = 0 testSrcTree = C_Src -[AMR_Adv_CF_o_2D] +[AMR_Adv_CF_o_2D] buildDir = Tests/FortranInterface/Advection_octree_F/Exec/SingleVortex inputFile = inputs.rt dim = 2 @@ -359,20 +359,6 @@ doVis = 0 outputFile = plot testSrcTree = C_Src -[EB_MacProj_3D] -buildDir = Tests/LinearSolvers/MAC_Projection_EB -inputFile = inputs_3d -dim = 3 -restartTest = 0 -useMPI = 1 -numprocs = 4 -useOMP = 0 -numthreads = 1 -compileTest = 0 -doVis = 0 -outputFile = plt00000 -testSrcTree = C_Src - [Redistribute] buildDir = Tests/Particles/Redistribute inputFile = inputs.rt