From c1cd7ab012a1fbe599c3643121b5ef83a3b42b88 Mon Sep 17 00:00:00 2001 From: Debojyoti Ghosh Date: Mon, 4 Nov 2024 12:20:47 -0800 Subject: [PATCH] Implicit Field Solve Preconditioner based on Curl-Curl Operator (#5286) Implemented a preconditioner for the implicit E-field solve using the AMReX curl-curl operator and the MLMG solver. + Introduced a `Preconditioner` base class that defines the action of a preconditioner for the JFNK algorithm. + Implemented the `CurlCurlMLMGPC` that uses the multigrid solution for the curl-curl operator (implemented in `AMReX`) to precondition the E-field JFNK solve. Other changes needed for this: + Partially implemented a mapping between WarpX field boundary types and AMReX's linear operator boundary types. + Added some functionalities to `ImplicitSolver` class that allows preconditioners to access `WarpX` info (like `Geometry`, boundaries, etc). Some premilinary wall times for: ``` Test: inputs_vandb_2d Grid: 160 X 160 dt: 0.125/wpe = 2.22e-18 (dt_CFL = 7.84e-19 s, CFL = 2.83) Time iterations: 20 Solver parameters: newton.max_iterations = 10 newton.relative_tolerance = 1.0e-12 newton.absolute_tolerance = 0.0 gmres.max_iterations = 1000 gmres.relative_tolerance = 1.0e-8 gmres.absolute_tolerance = 0.0 Avg GMRES iterations: ~3 (wPC), ~27 (noPC) ``` with `32^2` particles per cell: ``` Lassen (MPI + CUDA) ------------------- Box GPU Walltime (s) wPC noPC 1 1 2324.7 15004.1 4 1 2306.8 14356.8 4 4 758.9 3647.3 Dane (MPI + OMP) ---------------- Box CPU Threads Walltime (s) wPC noPC 1 1 1 6709.3 43200.0* 1 1 2 3279.1 22296.1 1 1 4 1696.3 11613.2 1 1 8 1085.0 6911.4 1 1 16 724.3 4729.0 4 1 1 5525.9 33288.8 16 1 1 4419.4 28467.8 4 4 1 1324.4 9121.1 16 16 1 524.9 3658.8 * 43200.0 seconds is 12 hours (max job duration on Dane); the simulation was almost done (started the 20th step). ``` with `10^2` particles per cell: ``` Lassen (MPI + CUDA) ------------------- Box GPU Walltime (s) wPC noPC 1 1 365.0 1443.5 4 1 254.1 927.8 4 4 133.1 301.5 Dane (MPI + OMP) ---------------- Box CPU Threads Walltime (s) wPC noPC 1 1 1 440.8 2360.5 1 1 2 241.7 1175.8 1 1 4 129.3 727.0 1 1 8 94.2 407.5 1 1 16 74.3 245.6 4 1 1 393.3 1932.5 16 1 1 337.6 1618.7 4 4 1 92.2 479.1 16 16 1 58.1 192.6 ``` --------- Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Co-authored-by: Remi Lehe Co-authored-by: Justin Angus Co-authored-by: Weiqun Zhang --- .../ImplicitSolvers/CMakeLists.txt | 1 + .../ImplicitSolvers/ImplicitSolver.H | 24 +- .../ImplicitSolvers/ImplicitSolver.cpp | 60 +++ .../FieldSolver/ImplicitSolvers/Make.package | 1 + .../ImplicitSolvers/ThetaImplicitEM.H | 5 +- .../ImplicitSolvers/ThetaImplicitEM.cpp | 4 +- .../ImplicitSolvers/WarpXSolverVec.H | 3 +- .../ImplicitSolvers/WarpXSolverVec.cpp | 2 + Source/NonlinearSolvers/CurlCurlMLMGPC.H | 355 ++++++++++++++++++ Source/NonlinearSolvers/JacobianFunctionMF.H | 51 ++- Source/NonlinearSolvers/NewtonSolver.H | 23 +- Source/NonlinearSolvers/Preconditioner.H | 100 +++++ Source/WarpX.H | 10 + 13 files changed, 615 insertions(+), 24 deletions(-) create mode 100644 Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp create mode 100644 Source/NonlinearSolvers/CurlCurlMLMGPC.H create mode 100644 Source/NonlinearSolvers/Preconditioner.H 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..ea9af6e2298 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. * @@ -9,9 +9,11 @@ #include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" #include "NonlinearSolvers/NonlinearSolverLibrary.H" +#include "Utils/WarpXAlgorithmSelection.H" #include #include +#include /** * \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..a6cbdfd307d --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp @@ -0,0 +1,60 @@ +#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->Geom(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"); + } else { + WARPX_ABORT_WITH_MESSAGE("Invalid value for 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..69d56c6ddc5 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H @@ -8,13 +8,12 @@ #define THETA_IMPLICIT_EM_H_ #include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" +#include "ImplicitSolver.H" #include #include #include -#include "ImplicitSolver.H" - /** @file * Theta-implicit electromagnetic time solver class. This is a fully implicit * algorithm where both the fields and particles are treated implicitly. @@ -79,7 +78,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..47d7310995c --- /dev/null +++ b/Source/NonlinearSolvers/CurlCurlMLMGPC.H @@ -0,0 +1,355 @@ +/* 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 "Fields.H" +#include "Utils/WarpXConst.H" +#include "Preconditioner.H" + +#include + +#include +#include +#include +#include +#include +#ifndef WARPX_DIM_1D_Z // currently not implemented in 1D +#include +#include +#include +#include +#endif + +/** + * \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; + + // Prohibit move and copy operations + CurlCurlMLMGPC(const CurlCurlMLMGPC&) = delete; + CurlCurlMLMGPC& operator=(const CurlCurlMLMGPC&) = delete; + CurlCurlMLMGPC(CurlCurlMLMGPC&&) noexcept = delete; + CurlCurlMLMGPC& operator=(CurlCurlMLMGPC&&) noexcept = delete; + + /** + * \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 = nullptr; + + 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(); + + private: + +}; + +template +void CurlCurlMLMGPC::printParameters() const +{ + using namespace amrex; + auto pc_name = getEnumNameString(PreconditionerType::pc_curl_curl_mlmg); + Print() << pc_name << " verbose: " << (m_verbose?"true":"false") << "\n"; + Print() << pc_name << " bottom verbose: " << (m_bottom_verbose?"true":"false") << "\n"; + Print() << pc_name << " max iter: " << m_max_iter << "\n"; + Print() << pc_name << " agglomeration: " << m_agglomeration << "\n"; + Print() << pc_name << " consolidation: " << m_consolidation << "\n"; + Print() << pc_name << " max_coarsening_level: " << m_max_coarsening_level << "\n"; + Print() << pc_name << " absolute tolerance: " << m_atol << "\n"; + Print() << pc_name << " relative tolerance: " << m_rtol << "\n"; + Print() << pc_name << " use GMRES: " << (m_use_gmres?"true":"false") << "\n"; + if (m_use_gmres) { + Print() << pc_name + << " use PC for GMRES: " + << (m_use_gmres_pc?"true":"false") << "\n"; + } +} + +template +void CurlCurlMLMGPC::readParameters() +{ + const amrex::ParmParse pp(amrex::getEnumNameString(PreconditionerType::pc_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->setFixedIter(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 " << amrex::getEnumNameString(PreconditionerType::pc_curl_curl_mlmg) + << ": dt = " << this->m_dt << ", " + << " coefficients: " + << "alpha = " << alpha << ", " + << "beta = " << beta << "\n"; + } +} + +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 +#if defined(WARPX_DIM_1D_Z) + // Missing dimensions is x,y in WarpX and y,z in AMReX + WARPX_ABORT_WITH_MESSAGE("CurlCurlMLMGPC not yet implemented for 1D"); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + // Missing dimension is y in WarpX and z in AMReX + Array solution { MultiFab(*x_mfarrvec[n][0], make_alias, 0, 1), + MultiFab(*x_mfarrvec[n][2], make_alias, 0, 1), + MultiFab(*x_mfarrvec[n][1], make_alias, 0, 1) }; + Array rhs { MultiFab(*b_mfarrvec[n][0], make_alias, 0, 1), + MultiFab(*b_mfarrvec[n][2], make_alias, 0, 1), + MultiFab(*b_mfarrvec[n][1], make_alias, 0, 1) }; +#elif defined(WARPX_DIM_3D) + Array solution { MultiFab(*x_mfarrvec[n][0], make_alias, 0, 1), + MultiFab(*x_mfarrvec[n][1], make_alias, 0, 1), + MultiFab(*x_mfarrvec[n][2], make_alias, 0, 1) }; + Array rhs { MultiFab(*b_mfarrvec[n][0], make_alias, 0, 1), + MultiFab(*b_mfarrvec[n][1], make_alias, 0, 1), + MultiFab(*b_mfarrvec[n][2], make_alias, 0, 1) }; +#endif + +// 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 + } +} + +#endif diff --git a/Source/NonlinearSolvers/JacobianFunctionMF.H b/Source/NonlinearSolvers/JacobianFunctionMF.H index d5c2b6cbac9..a3222214381 100644 --- a/Source/NonlinearSolvers/JacobianFunctionMF.H +++ b/Source/NonlinearSolvers/JacobianFunctionMF.H @@ -7,6 +7,9 @@ #ifndef JacobianFunctionMF_H_ #define JacobianFunctionMF_H_ +#include "CurlCurlMLMGPC.H" +#include + /** * \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 != PreconditionerType::none) { + m_preCond->printParameters(); + } } - void define( const T&, Ops* ); + void define( const T&, Ops*, const PreconditionerType& ); 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; + PreconditionerType m_pc_type = PreconditionerType::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 ) +void JacobianFunctionMF::define ( const T& a_U, + Ops* a_ops, + const PreconditionerType& 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 != PreconditionerType::none); + if (m_usePreCond) { + m_pc_type = a_pc_type; + if (m_pc_type == PreconditionerType::pc_curl_curl_mlmg) { + m_preCond = std::make_unique>(); + } else { + std::stringstream convergenceMsg; + convergenceMsg << "JacobianFunctionMF::define(): " << amrex::getEnumNameString(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..9c73c44e69e 100644 --- a/Source/NonlinearSolvers/NewtonSolver.H +++ b/Source/NonlinearSolvers/NewtonSolver.H @@ -9,10 +9,11 @@ #include "NonlinearSolver.H" #include "JacobianFunctionMF.H" +#include "Preconditioner.H" +#include "Utils/TextMsg.H" #include #include -#include "Utils/TextMsg.H" #include @@ -79,6 +80,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: " << amrex::getEnumNameString(m_pc_type) << "\n"; + + m_linear_function->printParams(); } private: @@ -138,9 +142,12 @@ private: */ int m_gmres_restart_length = 30; + /** + * \brief Preconditioner type + */ + PreconditionerType m_pc_type = PreconditionerType::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 +191,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 +219,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 +340,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..191a48d00bc --- /dev/null +++ b/Source/NonlinearSolvers/Preconditioner.H @@ -0,0 +1,100 @@ +/* Copyright 2024 Debojyoti Ghosh + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PRECONDITIONER_H_ +#define WARPX_PRECONDITIONER_H_ + +#include + +/** + * \brief Types for preconditioners for field solvers + */ +AMREX_ENUM(PreconditionerType, pc_curl_curl_mlmg, none); + +/** + * \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& a_x, const T& a_b) = 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 bad63cd44d9..a635196d044 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -113,6 +113,16 @@ public: [[nodiscard]] int Verbose () const { return verbose; } + [[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);