From f800188e73f20ded2264d4e880d50795cb79d31c Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 19 Dec 2023 12:46:12 -0800 Subject: [PATCH] WIP: curl of curl --- Src/LinearSolvers/CMakeLists.txt | 7 + Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H | 412 ++++++++++++++++++++++ Src/LinearSolvers/MLMG/AMReX_MLLinOp.H | 19 +- Src/LinearSolvers/MLMG/Make.package | 4 + 4 files changed, 430 insertions(+), 12 deletions(-) create mode 100644 Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H diff --git a/Src/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt index c2851d49959..cbc76d5d247 100644 --- a/Src/LinearSolvers/CMakeLists.txt +++ b/Src/LinearSolvers/CMakeLists.txt @@ -64,6 +64,13 @@ foreach(D IN LISTS AMReX_SPACEDIM) ) endif () + if (D EQUAL 3) + target_sources(amrex_${D}d + PRIVATE + MLMG/AMReX_MLCurlCurl.H + ) + endif () + if (AMReX_EB AND NOT D EQUAL 1) target_sources(amrex_${D}d PRIVATE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H new file mode 100644 index 00000000000..62faf1d47c6 --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H @@ -0,0 +1,412 @@ +#ifndef AMREX_ML_CURL_CURL_H_ +#define AMREX_ML_CURL_CURL_H_ +#include + +#include + +namespace amrex { + +/** + * \brief curl (alpha curl E) + beta E = rhs + * + * Here E is an Array of MultiFabs on staggered grid, alpha is a positive + * scalar, and beta is a non-negative scalar. + * + * TODO: If beta is zero, the system could be singular. + * TODO: Make sure isCellCentered in MLMG makes sense. + */ +template +class MLCurlCurlT + : public MLLinOpT +{ +public: + using RT = typename MLLinOpT::RT; + using BCType = typename MLLinOpT::BCType; + using BCMode = typename MLLinOpT::BCMode; + using StateMode = typename MLLinOpT::StateMode; + using Location = typename MLLinOpT::Location; + + MLCurlCurlT () = default; + MLCurlCurlT (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info = LPInfo()); + + MLCurlCurlT (const MLCurlCurlT&) = delete; + MLCurlCurlT (MLCurlCurlT&&) = delete; + MLCurlCurlT& operator= (const MLCurlCurlT&) = delete; + MLCurlCurlT& operator= (MLCurlCurlT&&) = delete; + + void define (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info = LPInfo()); + + void setScalars (RT a_alpha, RT a_beta) noexcept; + + [[nodiscard]] std::string name () const override { + return std::string("curl of curl"); + } + + void setLevelBC (int amrlev, const MF* levelbcdata, + const MF* robinbc_a = nullptr, + const MF* robinbc_b = nullptr, + const MF* robinbc_f = nullptr) override; + + void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const override; + + void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override; + + void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const override; + + void interpolationAmr (int famrlev, MF& fine, const MF& crse, + IntVect const& nghost) const override; + + void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs, + const MF& fine_sol, const MF& fine_rhs) override; + + void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode, + StateMode s_mode, const MLMGBndryT* bndry=nullptr) const override; + + void smooth (int amrlev, int mglev, MF& sol, const MF& rhs, + bool skip_fillboundary=false) const override; + + void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, + const MF* crse_bcdata=nullptr) override; + + void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, + const MF& b, BCMode bc_mode, + const MF* crse_bcdata=nullptr) override; + + void reflux (int crse_amrlev, + MF& res, const MF& crse_sol, const MF& crse_rhs, + MF& fine_res, MF& fine_sol, const MF& fine_rhs) const override; + + void compFlux (int amrlev, const Array& fluxes, + MF& sol, Location loc) const override; + + void compGrad (int amrlev, const Array& grad, + MF& sol, Location loc) const override; + + void applyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const override {} + + void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const override {} + + void prepareForSolve () override; + + [[nodiscard]] bool isSingular (int /*amrlev*/) const override { return false; } + [[nodiscard]] bool isBottomSingular () const override { return false; } + + RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const override; + + std::unique_ptr> makeNLinOp (int /*grid_size*/) const override { return nullptr; } + + [[nodiscard]] RT normInf (int amrlev, MF const& mf, bool local) const override; + + void averageDownAndSync (Vector& sol) const override; + + void avgDownResAmr (int clev, MF& cres, MF const& fres) const override; + + void avgDownResMG (int clev, MF& cres, MF const& fres) const override; + + void make (Vector >& mf, IntVect const& ng) const override; + + [[nodiscard]] MF make (int amrlev, int mglev, IntVect const& ng) const override; + + [[nodiscard]] MF makeAlias (MF const& mf) const override; + + [[nodiscard]] MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const override; + + [[nodiscard]] MF makeCoarseAmr (int famrlev, IntVect const& ng) const override; + +private: + + void applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode s_mode) const; + + RT m_alpha = std::numeric_limits::lowest(); + RT m_beta = std::numeric_limits::lowest(); + Array m_etype{IntVect(0,1,1), + IntVect(1,0,1), + IntVect(1,1,0)}; + static constexpr int m_ncomp = 1; + using mf_type = typename MF::value_type; +}; + +template +MLCurlCurlT::MLCurlCurlT (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info) +{ + define(a_geom, a_grids, a_dmap, a_info); +} + +template +void MLCurlCurlT::define (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info) +{ + MLLinOpT::define(a_geom, a_grids, a_dmap, a_info, {}); +} + +template +void MLCurlCurlT::setScalars (RT a_alpha, RT a_beta) noexcept +{ + m_alpha = a_alpha; + m_beta = a_beta; +} + +template +void MLCurlCurlT::setLevelBC (int amrlev, const MF* levelbcdata, // TODO + const MF* robinbc_a, + const MF* robinbc_b, + const MF* robinbc_f) +{ + amrex::ignore_unused(amrlev, levelbcdata, robinbc_a, robinbc_b, robinbc_f); +} + +template +void +MLCurlCurlT::restriction (int amrlev, int cmglev, MF& crse, MF& fine) const +{ + amrex::ignore_unused(amrlev, cmglev, crse, fine); + amrex::Abort("MLCurlCurlT::restriction: TODO"); +} + +template +void +MLCurlCurlT::interpolation (int amrlev, int fmglev, MF& fine, + const MF& crse) const +{ + amrex::ignore_unused(amrlev, fmglev, fine, crse); + amrex::Abort("MLCurlCurlT::interpolation: TODO"); +} + +template +void +MLCurlCurlT::interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const +{ + amrex::ignore_unused(amrlev, fmglev, fine, crse); + amrex::Abort("MLCurlCurlT::interpAssign: TODO"); +} + +template +void +MLCurlCurlT::interpolationAmr (int famrlev, MF& fine, const MF& crse, + IntVect const& nghost) const +{ + amrex::ignore_unused(famrlev, fine, crse, nghost); + amrex::Abort("MLCurlCurlT::interpolationAmr: TODO"); +} + +template +void +MLCurlCurlT::averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs, + const MF& fine_sol, const MF& fine_rhs) +{ + amrex::ignore_unused(camrlev, crse_sol, crse_rhs, fine_sol, fine_rhs); + amrex::Abort("MLCurlCurlT::averageDownSolutionRHS: TODO"); +} + +template +void +MLCurlCurlT::apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode, + StateMode s_mode, const MLMGBndryT* /*bndry*/) const +{ + applyBC(amrlev, mglev, in, bc_mode, s_mode); + + amrex::Warning("MLCurlCurlT::apply: TODO"); +} + +template +void +MLCurlCurlT::smooth (int amrlev, int mglev, MF& sol, const MF& rhs, + bool skip_fillboundary) const +{ + amrex::ignore_unused(amrlev, mglev, sol, rhs, skip_fillboundary); + amrex::Abort("MLCurlCurlT::smooth: TODO"); +} + +template +void +MLCurlCurlT::solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, + const MF* /*crse_bcdata*/) +{ + BL_PROFILE("MLCurlCurl::solutionResidual()"); + const int mglev = 0; + apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + amrex::Xpay(resid, RT(-1.0), b, 0, 0, m_ncomp, IntVect(0)); + } +} + +template +void +MLCurlCurlT::correctionResidual (int amrlev, int mglev, MF& resid, MF& x, + const MF& b, BCMode bc_mode, + const MF* crse_bcdata) +{ + amrex::ignore_unused(amrlev, mglev, resid, x, b, bc_mode, crse_bcdata); + amrex::Abort("MLCurlCurlT::correctionResidual: TODO"); +} + +template +void +MLCurlCurlT::reflux (int crse_amrlev, + MF& res, const MF& crse_sol, const MF& crse_rhs, + MF& fine_res, MF& fine_sol, const MF& fine_rhs) const +{ + amrex::ignore_unused(crse_amrlev, res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs); + amrex::Abort("MLCurlCurlT::reflux: TODO"); +} + +template +void +MLCurlCurlT::compFlux (int amrlev, const Array& fluxes, + MF& sol, Location loc) const +{ + amrex::ignore_unused(amrlev, fluxes, sol, loc); + amrex::Abort("MLCurlCurlT::compFlux: TODO"); +} + +template +void +MLCurlCurlT::compGrad (int amrlev, const Array& grad, + MF& sol, Location loc) const +{ + amrex::ignore_unused(amrlev, grad, sol, loc); + amrex::Abort("MLCurlCurlT::compGrad: TODO"); +} + +template +void +MLCurlCurlT::prepareForSolve () +{ + amrex::Warning("MLCurlCurlT::prepareForSolve: TODO"); +} + +template +typename MLCurlCurlT::RT +MLCurlCurlT::xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const +{ + amrex::ignore_unused(amrlev, mglev, x, y, local); + amrex::Abort("MLCurlCurlT::xdoty: TODO"); + return 0; +} + +template +typename MLCurlCurlT::RT +MLCurlCurlT::normInf (int /*amrlev*/, MF const& mf, bool local) const +{ + return amrex::norminf(mf, 0, m_ncomp, IntVect(0), local); +} + +template +void +MLCurlCurlT::averageDownAndSync (Vector& sol) const +{ + amrex::ignore_unused(sol); + amrex::Abort("MLCurlCurlT::averageDownAndSync: TODO"); +} + +template +void +MLCurlCurlT::avgDownResAmr (int clev, MF& cres, MF const& fres) const +{ + amrex::ignore_unused(clev, cres, fres); + amrex::Abort("MLCurlCurlT::avgDownResAmr: TODO"); +} + +template +void +MLCurlCurlT::avgDownResMG (int clev, MF& cres, MF const& fres) const +{ + amrex::ignore_unused(clev, cres, fres); + amrex::Abort("MLCurlCurlT::avgDownResMG: TODO"); +} + +template +void +MLCurlCurlT::make (Vector >& mf, IntVect const& ng) const +{ + MLLinOpT::make(mf, ng); +} + +template +MF +MLCurlCurlT::make (int amrlev, int mglev, IntVect const& ng) const +{ + MF r; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + r[idim].define(amrex::convert(this->m_grids[amrlev][mglev], m_etype[idim]), + this->m_dmap[amrlev][mglev], m_ncomp, ng, MFInfo(), + *(this->m_factory)[amrlev][mglev]); + } + return r; +} + +template +MF +MLCurlCurlT::makeAlias (MF const& mf) const +{ + MF r; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + r[idim] = mf_type(mf[idim], amrex::make_alias, 0, mf[idim].nComp()); + } + return r; +} + +template +MF +MLCurlCurlT::makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const +{ + BoxArray cba = this->m_grids[amrlev][mglev]; + IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[mglev]; + cba.coarsen(ratio); + + MF r; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + r[idim].define(amrex::convert(cba, m_etype[idim]), + this->m_dmap[amrlev][mglev], m_ncomp, ng); + } + return r; +} + +template +MF +MLCurlCurlT::makeCoarseAmr (int famrlev, IntVect const& ng) const +{ + BoxArray cba = this->m_grids[famrlev][0]; + IntVect ratio(this->AMRRefRatio(famrlev-1)); + cba.coarsen(ratio); + + MF r; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + r[idim].define(amrex::convert(cba, m_etype[idim]), + this->m_dmap[famrlev][0], m_ncomp, ng); + } + return r; +} + +template +void +MLCurlCurlT::applyBC (int amrlev, int mglev, MF& in, BCMode /*bc_mode*/, + StateMode /*s_mode*/) const +{ + Vector mfs{AMREX_D_DECL(&(in[0]),&(in[1]),&(in[2]))}; + FillBoundary(mfs, this->m_geom[amrlev][mglev].periodicity()); +} + +template <> +struct LinOpData , void> +{ + using fab_type = typename MultiFab::fab_type; + using value_type = typename MultiFab::value_type; +}; + +using MLCurlCurl = MLCurlCurlT >; + +} + +#endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index f0dca07f3ab..90c9d1b8f3f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H @@ -1384,18 +1384,13 @@ template void MLLinOpT::make (Vector >& mf, IntVect const& ng) const { - if constexpr (IsMultiFabLike_v) { - mf.clear(); - mf.resize(m_num_amr_levels); - for (int alev = 0; alev < m_num_amr_levels; ++alev) { - mf[alev].resize(m_num_mg_levels[alev]); - for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { - mf[alev][mlev] = make(alev, mlev, ng); - } + mf.clear(); + mf.resize(m_num_amr_levels); + for (int alev = 0; alev < m_num_amr_levels; ++alev) { + mf[alev].resize(m_num_mg_levels[alev]); + for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { + mf[alev][mlev] = make(alev, mlev, ng); } - } else { - amrex::ignore_unused(mf, ng); - amrex::Abort("MLLinOpT::make: how did we get here?"); } } @@ -1483,8 +1478,8 @@ template void MLLinOpT::avgDownResMG (int clev, MF& cres, MF const& fres) const { - const int ncomp = this->getNComp(); if constexpr (amrex::IsFabArray::value) { + const int ncomp = this->getNComp(); #ifdef AMREX_USE_EB if (!fres.isAllRegular()) { if constexpr (std::is_same()) { diff --git a/Src/LinearSolvers/MLMG/Make.package b/Src/LinearSolvers/MLMG/Make.package index d66d64ec0eb..4a354a751af 100644 --- a/Src/LinearSolvers/MLMG/Make.package +++ b/Src/LinearSolvers/MLMG/Make.package @@ -85,5 +85,9 @@ ifneq ($(BL_NO_FORT),TRUE) F90EXE_sources += AMReX_MLLinOp_nd.F90 endif +ifeq ($(DIM),3) + CEXE_headers += AMReX_MLCurlCurl.H +endif + VPATH_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/MLMG INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/MLMG