diff --git a/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver.json b/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver.json index 0ca6bde570a..eef503e7274 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver.json +++ b/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver.json @@ -17,4 +17,4 @@ "particle_position_z": 10066.329600000008, "particle_weight": 19969036501.910976 } -} \ No newline at end of file +} diff --git a/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt b/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt index 6e16f19084c..04abc9d3e91 100644 --- a/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt +++ b/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt @@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS) warpx_set_suffix_dims(SD ${D}) target_sources(lib_${SD} PRIVATE + ImplicitSolver.cpp SemiImplicitEM.cpp ThetaImplicitEM.cpp WarpXImplicitOps.cpp diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H index 88ad6a058fd..09165519f64 100644 --- a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H +++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H @@ -1,4 +1,4 @@ -/* Copyright 2024 Justin Angus +/* Copyright 2024 Justin Angus, Debojyoti Ghosh * * This file is part of WarpX. * @@ -7,11 +7,13 @@ #ifndef Implicit_Solver_H_ #define Implicit_Solver_H_ -#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" -#include "NonlinearSolvers/NonlinearSolverLibrary.H" - #include #include +#include + +#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" +#include "NonlinearSolvers/NonlinearSolverLibrary.H" +#include "Utils/WarpXAlgorithmSelection.H" /** * \brief Base class for implicit time solvers. The base functions are those @@ -85,6 +87,16 @@ public: int a_nl_iter, bool a_from_jacobian ) = 0; + [[nodiscard]] virtual amrex::Real theta () const { return 1.0; } + + [[nodiscard]] int numAMRLevels () const { return m_num_amr_levels; } + + [[nodiscard]] const amrex::Geometry& GetGeometry (int) const; + [[nodiscard]] const amrex::Array& GetFieldBoundaryLo () const; + [[nodiscard]] const amrex::Array& GetFieldBoundaryHi () const; + [[nodiscard]] amrex::Array GetLinOpBCLo () const; + [[nodiscard]] amrex::Array GetLinOpBCHi () const; + protected: /** @@ -94,6 +106,11 @@ protected: bool m_is_defined = false; + /** + * \brief Number of AMR levels + */ + int m_num_amr_levels = 1; + /** * \brief Nonlinear solver type and object */ @@ -140,6 +157,11 @@ protected: } + /** + * \brief Convert from WarpX FieldBoundaryType to amrex::LinOpBCType + */ + [[nodiscard]] amrex::Array convertFieldBCToLinOpBC ( const amrex::Array& ) const; + }; #endif diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp new file mode 100644 index 00000000000..76f5509e4a6 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp @@ -0,0 +1,58 @@ +#include "ImplicitSolver.H" +#include "WarpX.H" + +using namespace amrex; + +const Geometry& ImplicitSolver::GetGeometry (const int a_lvl) const +{ + AMREX_ASSERT((a_lvl >= 0) && (a_lvl < m_num_amr_levels)); + return m_WarpX->GetGeometry(a_lvl); +} + +const Array& ImplicitSolver::GetFieldBoundaryLo () const +{ + return m_WarpX->GetFieldBoundaryLo(); +} + +const Array& ImplicitSolver::GetFieldBoundaryHi () const +{ + return m_WarpX->GetFieldBoundaryHi(); +} + +Array ImplicitSolver::GetLinOpBCLo () const +{ + return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo()); +} + +Array ImplicitSolver::GetLinOpBCHi () const +{ + return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryHi()); +} + +Array ImplicitSolver::convertFieldBCToLinOpBC (const Array& a_fbc) const +{ + Array lbc; + for (auto& bc : lbc) { bc = LinOpBCType::interior; } + for (int i = 0; i < AMREX_SPACEDIM; i++) { + if (a_fbc[i] == FieldBoundaryType::PML) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } else if (a_fbc[i] == FieldBoundaryType::Periodic) { + lbc[i] = LinOpBCType::Periodic; + } else if (a_fbc[i] == FieldBoundaryType::PEC) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } else if (a_fbc[i] == FieldBoundaryType::PMC) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } else if (a_fbc[i] == FieldBoundaryType::Damped) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } else if (a_fbc[i] == FieldBoundaryType::Absorbing_SilverMueller) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } else if (a_fbc[i] == FieldBoundaryType::Neumann) { + lbc[i] = LinOpBCType::Neumann; + } else if (a_fbc[i] == FieldBoundaryType::None) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } else if (a_fbc[i] == FieldBoundaryType::Open) { + WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType"); + } + } + return lbc; +} diff --git a/Source/FieldSolver/ImplicitSolvers/Make.package b/Source/FieldSolver/ImplicitSolvers/Make.package index a4543f94dd3..16cd4003490 100644 --- a/Source/FieldSolver/ImplicitSolvers/Make.package +++ b/Source/FieldSolver/ImplicitSolvers/Make.package @@ -1,3 +1,4 @@ +CEXE_sources += ImplicitSolver.cpp CEXE_sources += SemiImplicitEM.cpp CEXE_sources += ThetaImplicitEM.cpp CEXE_sources += WarpXImplicitOps.cpp diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H index aba66782154..a645f66df91 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H @@ -79,7 +79,7 @@ public: int a_nl_iter, bool a_from_jacobian ) override; - [[nodiscard]] amrex::Real theta () const { return m_theta; } + [[nodiscard]] amrex::Real theta () const override { return m_theta; } private: diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp index 4cd5de4f24f..e5b8431a930 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp @@ -19,6 +19,7 @@ void ThetaImplicitEM::Define ( WarpX* const a_WarpX ) // Retain a pointer back to main WarpX class m_WarpX = a_WarpX; + m_num_amr_levels = 1; // Define E and Eold vectors m_E.Define( m_WarpX, "Efield_fp" ); @@ -26,8 +27,7 @@ void ThetaImplicitEM::Define ( WarpX* const a_WarpX ) // Define B_old MultiFabs using ablastr::fields::Direction; - const int num_levels = 1; - for (int lev = 0; lev < num_levels; ++lev) { + for (int lev = 0; lev < m_num_amr_levels; ++lev) { const auto& ba_Bx = m_WarpX->m_fields.get(FieldType::Bfield_fp, Direction{0}, lev)->boxArray(); const auto& ba_By = m_WarpX->m_fields.get(FieldType::Bfield_fp, Direction{1}, lev)->boxArray(); const auto& ba_Bz = m_WarpX->m_fields.get(FieldType::Bfield_fp, Direction{2}, lev)->boxArray(); diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H index 29c808b48cd..d864f239e42 100644 --- a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H +++ b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H @@ -75,6 +75,7 @@ public: void Define ( const WarpXSolverVec& a_solver_vec ) { assertIsDefined( a_solver_vec ); + m_num_amr_levels = a_solver_vec.m_num_amr_levels; Define( WarpXSolverVec::m_WarpX, a_solver_vec.getVectorType(), a_solver_vec.getScalarType() ); @@ -300,7 +301,7 @@ private: std::string m_scalar_type_name = "none"; static constexpr int m_ncomp = 1; - static constexpr int m_num_amr_levels = 1; + int m_num_amr_levels = 1; inline static bool m_warpx_ptr_defined = false; inline static WarpX* m_WarpX = nullptr; diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp index 6a0e6bb8a91..22c3b1d67c1 100644 --- a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp +++ b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp @@ -34,6 +34,8 @@ void WarpXSolverVec::Define ( WarpX* a_WarpX, m_warpx_ptr_defined = true; } + m_num_amr_levels = 1; + m_vector_type_name = a_vector_type_name; m_scalar_type_name = a_scalar_type_name; diff --git a/Source/NonlinearSolvers/CurlCurlMLMGPC.H b/Source/NonlinearSolvers/CurlCurlMLMGPC.H new file mode 100644 index 00000000000..c0a7045d05e --- /dev/null +++ b/Source/NonlinearSolvers/CurlCurlMLMGPC.H @@ -0,0 +1,453 @@ +/* Copyright 2024 Debojyoti Ghosh + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef CURL_CURL_MLMG_PC_H_ +#define CURL_CURL_MLMG_PC_H_ + +#include +#include +#include +#include +#include + +// currently not implemented in 1D +#ifndef WARPX_DIM_1D_Z +#include +#include +#include +#include +#endif + +#include + +#include "Fields.H" +#include "Utils/WarpXConst.H" +#include "Preconditioner.H" + +/** + * \brief Curl-curl Preconditioner + * + * Preconditioner that solves the curl-curl equation for the E-field, given + * a RHS. Uses AMReX's curl-curl linear operator and multigrid solver. + * + * The equation solves for Eg in: + * curl ( alpha * curl ( Eg ) ) + beta * Eg = b + * where + * + alpha is a scalar + * + beta can either be a scalar that is constant in space or a MultiFab + * + Eg is the electric field. + * + b is a specified RHS with the same layout as Eg + * + * This class is templated on a solution-type class T and an operator class Ops. + * + * The Ops class must have the following function: + * + Return number of AMR levels + * + Return the amrex::Geometry object given an AMR level + * + Return hi and lo linear operator boundaries + * + Return the time step factor (theta) for the time integration scheme + * + * The T class must have the following functions: + * + Return underlying vector of amrex::MultiFab arrays + */ + +template +class CurlCurlMLMGPC : public Preconditioner +{ + public: + + using RT = typename T::value_type; + + /** + * \brief Default constructor + */ + CurlCurlMLMGPC () = default; + + /** + * \brief Default destructor + */ + ~CurlCurlMLMGPC () override = default; + + // Default move and copy operations + CurlCurlMLMGPC(const CurlCurlMLMGPC&) = default; + CurlCurlMLMGPC& operator=(const CurlCurlMLMGPC&) = default; + CurlCurlMLMGPC(CurlCurlMLMGPC&&) noexcept = default; + CurlCurlMLMGPC& operator=(CurlCurlMLMGPC&&) noexcept = default; + + /** + * \brief Define the preconditioner + */ + void Define (const T&, Ops*) override; + + /** + * \brief Update the preconditioner + */ + void Update (const T&) override; + + /** + * \brief Apply (solve) the preconditioner given a RHS + * + * Given a right-hand-side b, solve: + * A x = b + * where A is the linear operator, in this case, the curl-curl operator: + * A x = curl (alpha * curl (x) ) + beta * x + */ + void Apply (T&, const T&) override; + + /** + * \brief Print parameters + */ + void printParameters() const override; + + /** + * \brief Check if the nonlinear solver has been defined. + */ + [[nodiscard]] inline bool IsDefined () const override { return m_is_defined; } + + protected: + + using MFArr = amrex::Array; + + bool m_is_defined = false; + + bool m_verbose = true; + bool m_bottom_verbose = false; + bool m_agglomeration = true; + bool m_consolidation = true; + bool m_use_gmres = false; + bool m_use_gmres_pc = true; + + int m_max_iter = 10; + int m_max_coarsening_level = 30; + + RT m_atol = 1.0e-16; + RT m_rtol = 1.0e-4; + + Ops* m_ops; + + int m_num_amr_levels = 0; + amrex::Vector m_geom; + amrex::Vector m_grids; + amrex::Vector m_dmap; + amrex::IntVect m_gv; + +// currently not implemented in 1D +#ifndef WARPX_DIM_1D_Z + amrex::Array m_bc_lo; + amrex::Array m_bc_hi; + + std::unique_ptr m_info; + std::unique_ptr m_curl_curl; + std::unique_ptr> m_solver; + std::unique_ptr> m_gmres_solver; +#endif + + /** + * \brief Read parameters + */ + void readParameters(); + + /** + * \brief Define staggered-grid MultiFabs for AMReX's curl-curl operator given WarpX's fields + */ + void defineMLCCAMFFromWarpXAMF ( amrex::Array&, + const ablastr::fields::VectorField&, + const amrex::IntVect& ); + + /** + * \brief Copy to AMReX's staggered-grid MultiFabs from WarpX's staggered-grid MultiFabs + */ + void copyMLCCAMFFromWarpXAMF ( amrex::Array&, + const ablastr::fields::VectorField&, + const amrex::IntVect& ); + + /** + * \brief Copy to WarpX's staggered-grid MultiFabs from AMReX's staggered-grid MultiFabs + */ + void copyWarpXAMFFromMLCCAMF ( ablastr::fields::VectorField&, + const amrex::Array&, + const amrex::IntVect& ); + + private: + +}; + +template +void CurlCurlMLMGPC::printParameters() const +{ + using namespace amrex; + Print() << PreconditionerTypes::curl_curl_mlmg << " verbose: " << (m_verbose?"true":"false") << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " bottom verbose: " << (m_bottom_verbose?"true":"false") << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " max iter: " << m_max_iter << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " agglomeration: " << m_agglomeration << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " consolidation: " << m_consolidation << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " max_coarsening_level: " << m_max_coarsening_level << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " absolute tolerance: " << m_atol << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " relative tolerance: " << m_rtol << "\n"; + Print() << PreconditionerTypes::curl_curl_mlmg << " use GMRES: " << (m_use_gmres?"true":"false") << "\n"; + if (m_use_gmres) { + Print() << PreconditionerTypes::curl_curl_mlmg + << " use PC for GMRES: " + << (m_use_gmres_pc?"true":"false") << "\n"; + } +} + +template +void CurlCurlMLMGPC::readParameters() +{ + const amrex::ParmParse pp(PreconditionerTypes::curl_curl_mlmg); + pp.query("verbose", m_verbose); + pp.query("bottom_verbose", m_bottom_verbose); + pp.query("max_iter", m_max_iter); + pp.query("agglomeration", m_agglomeration); + pp.query("consolidation", m_consolidation); + pp.query("max_coarsening_level", m_max_coarsening_level); + pp.query("absolute_tolerance", m_atol); + pp.query("relative_tolerance", m_rtol); + pp.query("use_gmres", m_use_gmres); + pp.query("use_gmres_pc", m_use_gmres_pc); +} + +template +void CurlCurlMLMGPC::Define ( const T& a_U, + Ops* const a_ops ) +{ + using namespace amrex; + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !IsDefined(), + "CurlCurlMLMGPC::Define() called on defined object" ); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (a_ops != nullptr), + "CurlCurlMLMGPC::Define(): a_ops is nullptr" ); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + a_U.getArrayVecType()==warpx::fields::FieldType::Efield_fp, + "CurlCurlMLMGPC::Define() must be called with Efield_fp type"); + + m_ops = a_ops; + // read preconditioner parameters + readParameters(); + +// currently not implemented in 1D +#ifdef WARPX_DIM_1D_Z + WARPX_ABORT_WITH_MESSAGE("CurlCurlMLMGPC not yet implemented for 1D"); +#else + // create info object for curl-curl op + m_info = std::make_unique(); + m_info->setAgglomeration(m_agglomeration); + m_info->setConsolidation(m_consolidation); + m_info->setMaxCoarseningLevel(m_max_coarsening_level); + + // Get data vectors from a_U + auto& u_mfarrvec = a_U.getArrayVec(); + + // Set number of AMR levels and create geometry, grids, and + // distribution mapping vectors. + m_num_amr_levels = m_ops->numAMRLevels(); + m_geom.resize(m_num_amr_levels); + m_grids.resize(m_num_amr_levels); + m_dmap.resize(m_num_amr_levels); + for (int n = 0; n < m_num_amr_levels; n++) { + m_geom[n] = m_ops->GetGeometry(n); + m_dmap[n] = u_mfarrvec[n][0]->DistributionMap(); + + BoxArray ba = u_mfarrvec[n][0]->boxArray(); + m_grids[n] = ba.enclosedCells(); + } + + // Construct the curl-curl linear operator and set its BCs + m_curl_curl = std::make_unique(m_geom, m_grids, m_dmap, *m_info); + m_curl_curl->setDomainBC(m_ops->GetLinOpBCLo(), m_ops->GetLinOpBCHi()); + + // Dummy value for alpha and beta to avoid abort due to degenerate matrix by MLMG solver + m_curl_curl->setScalars(1.0, 1.0); + + // Construct the MLMG solver + m_solver = std::make_unique>(*m_curl_curl); + m_solver->setMaxIter(m_max_iter); + m_solver->setVerbose(static_cast(m_verbose)); + m_solver->setBottomVerbose(static_cast(m_bottom_verbose)); + + // If using GMRES solver, construct it + if (m_use_gmres) { + m_gmres_solver = std::make_unique>(*m_solver); + m_gmres_solver->usePrecond(m_use_gmres_pc); + m_gmres_solver->setPrecondNumIters(m_max_iter); + m_gmres_solver->setVerbose(static_cast(m_verbose)); + } +#endif + + m_is_defined = true; +} + +template +void CurlCurlMLMGPC::Update (const T& a_U) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + IsDefined(), + "CurlCurlMLMGPC::Update() called on undefined object" ); + + // a_U is not needed for a linear operator + amrex::ignore_unused(a_U); + + // set the coefficients alpha and beta for curl-curl op + const RT alpha = (m_ops->theta()*this->m_dt*PhysConst::c) * (m_ops->theta()*this->m_dt*PhysConst::c); + const RT beta = RT(1.0); + +// currently not implemented in 1D +#ifndef WARPX_DIM_1D_Z + m_curl_curl->setScalars(alpha, beta); +#endif + + if (m_verbose) { + amrex::Print() << "Updating " << PreconditionerTypes::curl_curl_mlmg + << ": dt = " << this->m_dt << ", " + << " coefficients: " + << "alpha = " << alpha << ", " + << "beta = " << beta << "\n"; + } +} + +template +void CurlCurlMLMGPC::defineMLCCAMFFromWarpXAMF ( amrex::Array& a_mfarr_mlcc, + const ablastr::fields::VectorField& a_mfarr_warpx, + const amrex::IntVect& a_ng ) +{ +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(a_mfarr_mlcc); + amrex::ignore_unused(a_mfarr_warpx); + amrex::ignore_unused(a_ng); + WARPX_ABORT_WITH_MESSAGE("CurlCurlMLMGPC::defineMLCCAMFFromWarpXAMF not yet implemented for 1D"); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + // In WarpX, missing dimension is y + // In AMReX, missing dimension is z + { + a_mfarr_mlcc[0].define( a_mfarr_warpx[0]->boxArray(), + a_mfarr_warpx[0]->DistributionMap(), + a_mfarr_warpx[0]->nComp(), + a_ng ); + } + { + a_mfarr_mlcc[1].define( a_mfarr_warpx[2]->boxArray(), + a_mfarr_warpx[2]->DistributionMap(), + a_mfarr_warpx[2]->nComp(), + a_ng ); + } + { + a_mfarr_mlcc[2].define( a_mfarr_warpx[1]->boxArray(), + a_mfarr_warpx[1]->DistributionMap(), + a_mfarr_warpx[1]->nComp(), + a_ng ); + } +#elif defined(WARPX_DIM_3D) + for (int d = 0; d < 3; d++) { + a_mfarr_mlcc[d].define( a_mfarr_warpx[d]->boxArray(), + a_mfarr_warpx[d]->DistributionMap(), + a_mfarr_warpx[d]->nComp(), + a_ng ); + } +#endif +} + +template +void CurlCurlMLMGPC::copyMLCCAMFFromWarpXAMF ( amrex::Array& a_dst, + const ablastr::fields::VectorField& a_src, + const amrex::IntVect& a_ng ) +{ +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(a_dst); + amrex::ignore_unused(a_src); + amrex::ignore_unused(a_ng); + WARPX_ABORT_WITH_MESSAGE("CurlCurlMLMGPC::copyMLCCAMFFromWarpXAMF not yet implemented for 1D"); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + // In WarpX, missing dimension is y + // In AMReX, missing dimension is z + amrex::MultiFab::Copy( a_dst[0], *(a_src[0]), 0, 0, a_dst[0].nComp(), a_ng); + amrex::MultiFab::Copy( a_dst[1], *(a_src[2]), 0, 0, a_dst[1].nComp(), a_ng); + amrex::MultiFab::Copy( a_dst[2], *(a_src[1]), 0, 0, a_dst[2].nComp(), a_ng); +#elif defined(WARPX_DIM_3D) + for (int d = 0; d < 3; d++) { + amrex::MultiFab::Copy( a_dst[d], *(a_src[d]), 0, 0, a_dst[d].nComp(), a_ng); + } +#endif +} + +template +void CurlCurlMLMGPC::copyWarpXAMFFromMLCCAMF ( ablastr::fields::VectorField& a_dst, + const amrex::Array& a_src, + const amrex::IntVect& a_ng ) +{ +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(a_dst); + amrex::ignore_unused(a_src); + amrex::ignore_unused(a_ng); + WARPX_ABORT_WITH_MESSAGE("CurlCurlMLMGPC::copyMLCCAMFFromWarpXAMF not yet implemented for 1D"); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + // In WarpX, missing dimension is y + // In AMReX, missing dimension is z + amrex::MultiFab::Copy( *(a_dst[0]), a_src[0], 0, 0, a_dst[0]->nComp(), a_ng); + amrex::MultiFab::Copy( *(a_dst[1]), a_src[2], 0, 0, a_dst[1]->nComp(), a_ng); + amrex::MultiFab::Copy( *(a_dst[2]), a_src[1], 0, 0, a_dst[2]->nComp(), a_ng); +#elif defined(WARPX_DIM_3D) + for (int d = 0; d < 3; d++) { + amrex::MultiFab::Copy( *(a_dst[d]), a_src[d], 0, 0, a_dst[d]->nComp(), a_ng); + } +#endif +} + +template +void CurlCurlMLMGPC::Apply (T& a_x, const T& a_b) +{ + // Given a right-hand-side b, solve: + // A x = b + // where A is the linear operator, in this case, the curl-curl + // operator: + // A x = curl (alpha * curl (x) ) + beta * x + + using namespace amrex; + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + IsDefined(), + "CurlCurlMLMGPC::Apply() called on undefined object" ); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + a_x.getArrayVecType()==warpx::fields::FieldType::Efield_fp, + "CurlCurlMLMGPC::Apply() - a_x must be Efield_fp type"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + a_b.getArrayVecType()==warpx::fields::FieldType::Efield_fp, + "CurlCurlMLMGPC::Apply() - a_b must be Efield_fp type"); + + // Get the data vectors + auto& b_mfarrvec = a_b.getArrayVec(); + auto& x_mfarrvec = a_x.getArrayVec(); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + ((b_mfarrvec.size() == m_num_amr_levels) && (x_mfarrvec.size() == m_num_amr_levels)), + "Error in CurlCurlMLMGPC::Apply() - mismatch in number of levels." ); + + for (int n = 0; n < m_num_amr_levels; n++) { + + // Copy initial guess to local object + Array solution, rhs; + defineMLCCAMFFromWarpXAMF(solution, x_mfarrvec[n], IntVect::TheUnitVector()); + defineMLCCAMFFromWarpXAMF(rhs, b_mfarrvec[n], IntVect::TheZeroVector()); + copyMLCCAMFFromWarpXAMF(solution, x_mfarrvec[n], IntVect::TheZeroVector()); + copyMLCCAMFFromWarpXAMF(rhs, b_mfarrvec[n], IntVect::TheZeroVector()); + +// currently not implemented in 1D +#ifndef WARPX_DIM_1D_Z + m_curl_curl->prepareRHS({&rhs}); + if (m_use_gmres) { + m_gmres_solver->solve(solution, rhs, m_rtol, m_atol); + } else { + m_solver->solve({&solution}, {&rhs}, m_rtol, m_atol); + } +#endif + + // Copy solution in local object to a_x + copyWarpXAMFFromMLCCAMF(x_mfarrvec[n], solution, IntVect::TheZeroVector()); + + } +} + +#endif diff --git a/Source/NonlinearSolvers/JacobianFunctionMF.H b/Source/NonlinearSolvers/JacobianFunctionMF.H index d5c2b6cbac9..11ef0828105 100644 --- a/Source/NonlinearSolvers/JacobianFunctionMF.H +++ b/Source/NonlinearSolvers/JacobianFunctionMF.H @@ -7,6 +7,9 @@ #ifndef JacobianFunctionMF_H_ #define JacobianFunctionMF_H_ +#include +#include "CurlCurlMLMGPC.H" + /** * \brief This is a linear function class for computing the action of a * Jacobian on a vector using a matrix-free finite-difference method. @@ -35,14 +38,18 @@ class JacobianFunctionMF inline void precond ( T& a_U, const T& a_X ) { - if (m_usePreCond) { a_U.zero(); } - else { a_U.Copy(a_X); } + if (m_usePreCond) { + a_U.zero(); + m_preCond->Apply(a_U, a_X); + } else { + a_U.Copy(a_X); + } } inline void updatePreCondMat ( const T& a_X ) { - amrex::ignore_unused(a_X); + if (m_usePreCond) { m_preCond->Update(a_X); } } inline @@ -133,15 +140,25 @@ class JacobianFunctionMF void curTime ( RT a_time ) { m_cur_time = a_time; + if (m_usePreCond) { m_preCond->CurTime(a_time); } } inline void curTimeStep ( RT a_dt ) { m_dt = a_dt; + if (m_usePreCond) { m_preCond->CurTimeStep(a_dt); } + } + + inline + void printParams () const + { + if (m_pc_type != "none") { + m_preCond->printParameters(); + } } - void define( const T&, Ops* ); + void define( const T&, Ops*, const std::string& ); private: @@ -151,16 +168,18 @@ class JacobianFunctionMF RT m_epsJFNK = RT(1.0e-6); RT m_normY0; RT m_cur_time, m_dt; - std::string m_pc_type; - T m_Z, m_Y0, m_R0, m_R; - Ops* m_ops; + std::string m_pc_type = "none"; + T m_Z, m_Y0, m_R0, m_R; + Ops* m_ops = nullptr; + std::unique_ptr> m_preCond = nullptr; }; template void JacobianFunctionMF::define ( const T& a_U, - Ops* a_ops ) + Ops* a_ops, + const std::string& a_pc_type ) { m_Z.Define(a_U); m_Y0.Define(a_U); @@ -169,6 +188,20 @@ void JacobianFunctionMF::define ( const T& a_U, m_ops = a_ops; + m_usePreCond = (a_pc_type != "none"); + if (m_usePreCond) { + m_pc_type = a_pc_type; + if (m_pc_type == PreconditionerTypes::curl_curl_mlmg) { + m_preCond = std::make_unique>(); + } else { + std::stringstream convergenceMsg; + convergenceMsg << "JacobianFunctionMF::define(): " << m_pc_type << + " is not a valid preconditioner type."; + WARPX_ABORT_WITH_MESSAGE(convergenceMsg.str()); + } + m_preCond->Define(a_U, a_ops); + } + m_is_defined = true; } diff --git a/Source/NonlinearSolvers/NewtonSolver.H b/Source/NonlinearSolvers/NewtonSolver.H index 742e139a5f5..8ec1cee2170 100644 --- a/Source/NonlinearSolvers/NewtonSolver.H +++ b/Source/NonlinearSolvers/NewtonSolver.H @@ -79,6 +79,9 @@ public: amrex::Print() << "GMRES max iterations: " << m_gmres_maxits << "\n"; amrex::Print() << "GMRES relative tolerance: " << m_gmres_rtol << "\n"; amrex::Print() << "GMRES absolute tolerance: " << m_gmres_atol << "\n"; + amrex::Print() << "Preconditioner type: " << m_pc_type << "\n"; + + m_linear_function->printParams(); } private: @@ -138,9 +141,12 @@ private: */ int m_gmres_restart_length = 30; + /** + * \brief Preconditioner type + */ + std::string m_pc_type = "none"; + mutable amrex::Real m_cur_time, m_dt; - mutable bool m_update_pc = false; - mutable bool m_update_pc_init = false; /** * \brief The linear function used by GMRES to compute A*v. @@ -184,7 +190,7 @@ void NewtonSolver::Define ( const Vec& a_U, m_ops = a_ops; m_linear_function = std::make_unique>(); - m_linear_function->define(m_F, m_ops); + m_linear_function->define(m_F, m_ops, m_pc_type); m_linear_solver = std::make_unique>>(); m_linear_solver->define(*m_linear_function); @@ -212,6 +218,9 @@ void NewtonSolver::ParseParameters () pp_gmres.query("absolute_tolerance", m_gmres_atol); pp_gmres.query("relative_tolerance", m_gmres_rtol); pp_gmres.query("max_iterations", m_gmres_maxits); + + const amrex::ParmParse pp_jac("jacobian"); + pp_jac.query("pc_type", m_pc_type); } template @@ -330,10 +339,7 @@ void NewtonSolver::EvalResidual ( Vec& a_F, m_linear_function->setBaseRHS(m_R); // update preconditioner - if (m_update_pc || m_update_pc_init) { - m_linear_function->updatePreCondMat(a_U); - } - m_update_pc_init = false; + m_linear_function->updatePreCondMat(a_U); // Compute residual: F(U) = U - b - R(U) a_F.Copy(a_U); diff --git a/Source/NonlinearSolvers/Preconditioner.H b/Source/NonlinearSolvers/Preconditioner.H new file mode 100644 index 00000000000..3f9b7d72dc0 --- /dev/null +++ b/Source/NonlinearSolvers/Preconditioner.H @@ -0,0 +1,104 @@ +/* Copyright 2024 Debojyoti Ghosh + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PRECONDITIONER_H_ +#define WARPX_PRECONDITIONER_H_ + +/** + * \brief Types for preconditioners for field solvers + */ +namespace PreconditionerTypes +{ + /** + * MLMG-based solver for the curl-curl operator + */ + const std::string curl_curl_mlmg = "pc_curl_curl_mlmg"; +} + +/** + * \brief Base class for preconditioners + * + * This class is templated on a solution-type class T and an operator class Ops. + * + * The Ops class must have the following function: + * (this will depend on the specific preconditioners inheriting from this class) + * + * The T class must have the following functions: + * (this will depend on the specific preconditioners inheriting from this class) + */ + +template +class Preconditioner +{ + public: + + using RT = typename T::value_type; + + /** + * \brief Default constructor + */ + Preconditioner () = default; + + /** + * \brief Default destructor + */ + virtual ~Preconditioner () = default; + + // Default move and copy operations + Preconditioner(const Preconditioner&) = default; + Preconditioner& operator=(const Preconditioner&) = default; + Preconditioner(Preconditioner&&) noexcept = default; + Preconditioner& operator=(Preconditioner&&) noexcept = default; + + /** + * \brief Define the preconditioner + */ + virtual void Define (const T&, Ops*) = 0; + + /** + * \brief Update the preconditioner + */ + virtual void Update ( const T& ) = 0; + + /** + * \brief Apply (solve) the preconditioner given a RHS + * + * Given a right-hand-side b, solve: + * A x = b + * where A is a linear operator. + */ + virtual void Apply (T&, const T&) = 0; + + /** + * \brief Check if the nonlinear solver has been defined. + */ + [[nodiscard]] virtual bool IsDefined () const = 0; + + /** + * \brief Print parameters + */ + virtual void printParameters() const { } + + /** + * \brief Set the current time. + */ + inline void CurTime (const RT a_time) { m_time = a_time; } + + /** + * \brief Set the current time step size. + */ + inline void CurTimeStep (const RT a_dt) { m_dt = a_dt; } + + protected: + + RT m_time = 0.0; + RT m_dt = 0.0; + + private: + +}; + +#endif diff --git a/Source/WarpX.H b/Source/WarpX.H index 83b1880f2b1..c873d911d12 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -113,6 +113,21 @@ public: [[nodiscard]] int Verbose () const { return verbose; } + [[nodiscard]] const amrex::Geometry& GetGeometry ( const int a_lev ) const + { + return GetInstance().Geom(a_lev); + } + + [[nodiscard]] const amrex::Array& GetFieldBoundaryLo () const + { + return field_boundary_lo; + } + + [[nodiscard]] const amrex::Array& GetFieldBoundaryHi () const + { + return field_boundary_hi; + } + void InitData (); void Evolve (int numsteps = -1);