diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index f3a74f969a0..fb53e81012f 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -86,37 +86,122 @@ Overall simulation parameters * ``explicit``: Use an explicit solver, such as the standard FDTD or PSATD - * ``implicit_picard``: Use an implicit solver with exact energy conservation that uses a Picard iteration to solve the system. - Note that this method is for demonstration only. It is inefficient and does not work well when - :math:`\omega_{pe} \Delta t` is close to or greater than one. - The method is described in `Angus et al., On numerical energy conservation for an implicit particle-in-cell method coupled with a binary Monte-Carlo algorithm for Coulomb collisions `__. - The version implemented is an updated version that is relativistically correct, including the relativistic gamma factor for the particles. - For exact energy conservation, ``algo.current_deposition = direct`` must be used with ``interpolation.galerkin_scheme = 0``, - and ``algo.current_deposition = Esirkepov`` must be used with ``interpolation.galerkin_scheme = 1`` (which is the default, in - which case charge will also be conserved). - - * ``semi_implicit_picard``: Use an energy conserving semi-implicit solver that uses a Picard iteration to solve the system. - Note that this method has the CFL limitation :math:`\Delta t < c/\sqrt( \sum_i 1/\Delta x_i^2 )`. It is inefficient and does not work well or at all when :math:`\omega_{pe} \Delta t` is close to or greater than one. + * ``theta_implicit_em``: Use a fully implicit electromagnetic solver with a time-biasing parameter theta bound between 0.5 and 1.0. Exact energy conservation is achieved using theta = 0.5. Maximal damping of high-k modes is obtained using theta = 1.0. Choices for the nonlinear solver include a Picard iteration scheme and particle-suppressed (PS) JNFK. + The algorithm itself is numerical stable for large time steps. That is, it does not require time steps that resolve the plasma period or the CFL condition for light waves. However, the practicality of using a large time step depends on the nonlinear solver. Note that the Picard solver is for demonstration only. It is inefficient and will most like not converge when + :math:`\omega_{pe} \Delta t` is close to or greater than one or when the CFL condition for light waves is violated. The PS-JFNK method must be used in order to use large time steps. However, the current implementation of PS-JFNK is still inefficient because the JFNK solver is not preconditioned and there is no use of the mass matrices to minimize the cost of a linear iteration. The time step is limited by how many cells a particle can cross in a time step (MPI-related) and by the need to resolve the relavent physics. + The Picard method is described in `Angus et al., On numerical energy conservation for an implicit particle-in-cell method coupled with a binary Monte-Carlo algorithm for Coulomb collisions `__. + The PS-JFNK method is described in `Angus et al., An implicit particle code with exact energy and charge conservation for electromagnetic studies of dense plasmas `__ . (The version implemented in WarpX is an updated version that includes the relativistic gamma factor for the particles.) Also see `Chen et al., An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm. `__ . + Exact energy conservation requires that the interpolation stencil used for the field gather match that used for the current deposition. ``algo.current_deposition = direct`` must be used with ``interpolation.galerkin_scheme = 0``, and ``algo.current_deposition = Esirkepov`` must be used with ``interpolation.galerkin_scheme = 1``. If using ``algo.current_deposition = villasenor``, the corresponding field gather routine will automatically be selected and the ``interpolation.galerkin_scheme`` flag does not need to be specified. The Esirkepov and villasenor deposition schemes are charge-conserving. + + * ``semi_implicit_em``: Use an approximately energy conserving semi-implicit electromagnetic solver. Choices for the nonlinear solver include a Picard iteration scheme and particle-suppressed JFNK. + Note that this method has the CFL limitation :math:`\Delta t < c/\sqrt( \sum_i 1/\Delta x_i^2 )`. The Picard solver for this method can only be expected to work well when :math:`\omega_{pe} \Delta t` is less than one. The method is described in `Chen et al., A semi-implicit, energy- and charge-conserving particle-in-cell algorithm for the relativistic Vlasov-Maxwell equations `__. - For energy conservation, ``algo.current_deposition = direct`` must be used with ``interpolation.galerkin_scheme = 0``, - and ``algo.current_deposition = Esirkepov`` must be used with ``interpolation.galerkin_scheme = 1`` (which is the default, in - which case charge will also be conserved). - -* ``algo.max_picard_iterations`` (`integer`, default: 10) - When `algo.evolve_scheme` is either `implicit_picard` or `semi_implicit_picard`, this sets the maximum number of Picard - itearations that are done each time step. - -* ``algo.picard_iteration_tolerance`` (`float`, default: 1.e-7) - When `algo.evolve_scheme` is either `implicit_picard` or `semi_implicit_picard`, this sets the convergence tolerance of - the iterations, the maximum of the relative change of the L2 norm of the field from one iteration to the next. - If this is set to zero, the maximum number of iterations will always be done with the change only calculated on the last - iteration (for a slight optimization). - -* ``algo.require_picard_convergence`` (`bool`, default: 1) - When `algo.evolve_scheme` is either `implicit_picard` or `semi_implicit_picard`, this sets whether the iteration each step - is required to converge. - If it is required, an abort is raised if it does not converge and the code then exits. - If not, then a warning is issued and the calculation continues. + Exact energy conservation requires that the interpolation stencil used for the field gather match that used for the current deposition. ``algo.current_deposition = direct`` must be used with ``interpolation.galerkin_scheme = 0``, and ``algo.current_deposition = Esirkepov`` must be used with ``interpolation.galerkin_scheme = 1``. If using ``algo.current_deposition = villasenor``, the corresponding field gather routine will automatically be selected and the ``interpolation.galerkin_scheme`` flag does not need to be specified. The Esirkepov and villasenor deposition schemes are charge-conserving. + +* ``implicit_evolve.theta`` (`float`, default: 0.5) + When `algo.evolve_scheme = theta_implicit_em`, the fields used on the RHS of the equations for the implicit advance + are computed as (1-theta)*E_{n} + theta*E_{n+1}. theta is bound between 0.5 and 1. The default value of theta = 0.5 + is needed for exact energy conservation. For theta > 0.5, high-k modes will be damped and the method will not be + exactly energy conserving, but the solver may perform better. + +* ``implicit_evolve.nonlinear_solver`` (`string`, default: None) + When `algo.evolve_scheme` is either `theta_implicit_em` or `semi_implicit_em`, this sets the nonlinear solver used + to advance the field-particle system in time. Options are `picard` or `newton`. + +* ``implicit_evolve.max_particle_iterations`` (`integer`, default: 21) + When `algo.evolve_scheme` is either `theta_implicit_em` or `semi_implicit_em` and `implicit_evolve.nonlinear_solver = newton` + , this sets the maximum number of iterations for the method used to obtain a self-consistent update of the particles at + each iteration in the JFNK process. + +* ``implicit_evolve.particle_tolerance`` (`float`, default: 1.e-10) + When `algo.evolve_scheme` is either `theta_implicit_em` or `semi_implicit_em` and `implicit_evolve.nonlinear_solver = newton` + , this sets the relative tolerance for the iterative method used to obtain a self-consistent update of the particles at + each iteration in the JFNK process. + +* ``picard.verbose`` (`bool`, default: 1) + When `implicit_evolve.nonlinear_solver = picard`, this sets the verbosity of the Picard solver. If true, then information + on the nonlinear error are printed to screen at each nonlinear iteration. + +* ``picard.require_convergence`` (`bool`, default: 1) + When `implicit_evolve.nonlinear_solver = picard`, this sets whether the Picard method is required to converge at each + time step. If it is required, an abort is raised if it does not converge and the code then exits. If not, then a warning + is issued and the calculation continues. + +* ``picard.maximum_iterations`` (`int`, default: 100) + When `implicit_evolve.nonlinear_solver = picard`, this sets the maximum iterations used by the Picard method. If + `picard.require_convergence = false`, the solution is considered converged if the iteration count reaches this value, + but a warning is issued. If `picard.require_convergence = true`, then an abort is raised if the iteration count reaches + this value. + +* ``picard.relative_tolerance`` (`float`, default: 1.0e-6) + When `implicit_evolve.nonlinear_solver = picard`, this sets the relative tolerance used by the Picard method for determining + convergence. The absolute error for the Picard method is the L2 norm of the difference of the solution vector between + two successive iterations. The relative error is the absolute error after iteration k > 1 divided by the absolute error + after the first iteration. The Picard method is considered converged when the relative error is below the relative tolerance. + This is the preferred means of determining convergence. + +* ``picard.absolute_tolerance`` (`float`, default: 0.0) + When `implicit_evolve.nonlinear_solver = picard`, this sets the absolute tolerance used by the Picard method for determining + convergence. The default value is 0.0, which means that the absolute tolerance is not used to determine converence. The + solution vector in the nonlinear solvers are in physical units rather than normalized ones. Thus, the absolute scale + of the problem can vary over many orders and magnitude depending on the problem. The relative tolerance is the preferred + means of determining convergence. + +* ``newton.verbose`` (`bool`, default: 1) + When `implicit_evolve.nonlinear_solver = newton`, this sets the verbosity of the Newton solver. If true, then information + on the nonlinear error are printed to screen at each nonlinear iteration. + +* ``newton.require_convergence`` (`bool`, default: 1) + When `implicit_evolve.nonlinear_solver = newton`, this sets whether the Newton method is required to converge at each + time step. If it is required, an abort is raised if it does not converge and the code then exits. If not, then a warning + is issued and the calculation continues. + +* ``newton.maximum_iterations`` (`int`, default: 100) + When `implicit_evolve.nonlinear_solver = newton`, this sets the maximum iterations used by the Newton method. If + `newton.require_convergence = false`, the solution is considered converged if the iteration count reaches this value, + but a warning is issued. If `newton.require_convergence = true`, then an abort is raised if the iteration count reaches + this value. + +* ``newton.relative_tolerance`` (`float`, default: 1.0e-6) + When `implicit_evolve.nonlinear_solver = newton`, this sets the relative tolerance used by the Newton method for determining + convergence. The absolute error for the Newton method is the L2 norm of the residual vector. The relative error is the + absolute error divided by the L2 norm of the initial residual associated with the initial guess. The Newton method is + considered converged when the relative error is below the relative tolerance. This is the preferred means of determining + convergence. + +* ``newton.absolute_tolerance`` (`float`, default: 0.0) + When `implicit_evolve.nonlinear_solver = newton`, this sets the absolute tolerance used by the Newton method for determining + convergence. The default value is 0.0, which means that the absolute tolerance is not used to determine converence. The + residual vector in the nonlinear solvers are in physical units rather than normalized ones. Thus, the absolute scale + of the problem can vary over many orders and magnitude depending on the problem. The relative tolerance is the preferred + means of determining convergence. + +* ``gmres.verbose_int`` (`int`, default: 2) + When `implicit_evolve.nonlinear_solver = newton`, this sets the verbosity of the AMReX::GMRES linear solver. The default + value of 2 gives maximumal verbosity and information about the residual are printed to the screen at each GMRES iteration. + +* ``gmres.restart_length`` (`int`, default: 30) + When `implicit_evolve.nonlinear_solver = newton`, this sets the iteration number at which to do a restart in AMReX::GMRES. + This parameter is used to save memory on building the Krylov subspace basis vectors for linear systems that are ill-conditioned + and require many iterations to converge. + +* ``gmres.relative_tolerance`` (`float`, default: 1.0e-4) + When `implicit_evolve.nonlinear_solver = newton`, this sets the relative tolerance used to determine convergence of the + AMReX::GMRES linear solver used to compute the Newton step in the JNFK process. The absolute error is the L2 norm of the + residual vector. The relative error is the absolute error divided by the L2 norm of the initial residual (typically equal + to the norm of the nonlinear residual from the end of the previous Newton iteration). The linear solver is considered + converged when the relative error is below the relative tolerance. This is the preferred means of determining convergence. + +* ``gmres.absolute_tolerance`` (`float`, default: 0.0) + When `implicit_evolve.nonlinear_solver = newton`, this sets the absolute tolerance used to determine converence of the + GMRES linear solver used to compute the Newton step in the JNFK process. The default value is 0.0, which means that the + absolute tolerance is not used to determine converence. The residual vector in the nonlinear/linear solvers are in physical + units rather than normalized ones. Thus, the absolute scale of the problem can vary over many orders and magnitude depending + on the problem. The relative tolerance is the preferred means of determining convergence. + +* ``gmres.maximum_iterations`` (`int`, default: 1000) + When `implicit_evolve.nonlinear_solver = newton`, this sets the maximum iterations used by the GMRES linear solver. The + solution to the linear system is considered converged if the iteration count reaches this value. * ``warpx.do_electrostatic`` (`string`) optional (default `none`) Specifies the electrostatic mode. When turned on, instead of updating diff --git a/Examples/Tests/Implicit/analysis_1d.py b/Examples/Tests/Implicit/analysis_1d.py index 0f00010a505..0e20b925df5 100755 --- a/Examples/Tests/Implicit/analysis_1d.py +++ b/Examples/Tests/Implicit/analysis_1d.py @@ -31,7 +31,7 @@ if re.match('SemiImplicitPicard_1d', fn): tolerance_rel = 2.5e-5 -elif re.match('ImplicitPicard_1d', fn): +elif re.match('ThetaImplicitPicard_1d', fn): # This case should have near machine precision conservation of energy tolerance_rel = 1.e-14 diff --git a/Examples/Tests/Implicit/analysis_vandb_2d.py b/Examples/Tests/Implicit/analysis_vandb_jfnk_2d.py similarity index 100% rename from Examples/Tests/Implicit/analysis_vandb_2d.py rename to Examples/Tests/Implicit/analysis_vandb_jfnk_2d.py diff --git a/Examples/Tests/Implicit/inputs_1d b/Examples/Tests/Implicit/inputs_1d index 465d9dd6965..3e57689b723 100644 --- a/Examples/Tests/Implicit/inputs_1d +++ b/Examples/Tests/Implicit/inputs_1d @@ -31,10 +31,18 @@ boundary.particle_hi = periodic ############ NUMERICS ########### ################################# +warpx.verbose = 1 warpx.const_dt = dt -algo.evolve_scheme = implicit_picard -algo.max_picard_iterations = 31 -algo.picard_iteration_tolerance = 0. +algo.evolve_scheme = theta_implicit_em + +implicit_evolve.nonlinear_solver = "picard" + +picard.verbose = true +picard.max_iterations = 31 +picard.relative_tolerance = 0.0 +picard.absolute_tolerance = 0.0 +picard.require_convergence = false + algo.current_deposition = esirkepov algo.field_gathering = energy-conserving algo.particle_shape = 2 diff --git a/Examples/Tests/Implicit/inputs_1d_semiimplicit b/Examples/Tests/Implicit/inputs_1d_semiimplicit index 4008a559588..07460e08be8 100644 --- a/Examples/Tests/Implicit/inputs_1d_semiimplicit +++ b/Examples/Tests/Implicit/inputs_1d_semiimplicit @@ -31,10 +31,18 @@ boundary.particle_hi = periodic ############ NUMERICS ########### ################################# +warpx.verbose = 1 warpx.const_dt = dt -algo.evolve_scheme = semi_implicit_picard -algo.max_picard_iterations = 5 -algo.picard_iteration_tolerance = 0. +algo.evolve_scheme = semi_implicit_em + +implicit_evolve.nonlinear_solver = "picard" + +picard.verbose = true +picard.max_iterations = 5 +picard.relative_tolerance = 0.0 +picard.absolute_tolerance = 0.0 +picard.require_convergence = false + algo.current_deposition = esirkepov algo.field_gathering = energy-conserving algo.particle_shape = 2 diff --git a/Examples/Tests/Implicit/inputs_vandb_2d b/Examples/Tests/Implicit/inputs_vandb_jfnk_2d similarity index 78% rename from Examples/Tests/Implicit/inputs_vandb_2d rename to Examples/Tests/Implicit/inputs_vandb_jfnk_2d index 33ce964710a..393a9d90330 100644 --- a/Examples/Tests/Implicit/inputs_vandb_2d +++ b/Examples/Tests/Implicit/inputs_vandb_jfnk_2d @@ -38,10 +38,32 @@ warpx.const_dt = dt warpx.use_filter = 0 algo.maxwell_solver = Yee -algo.evolve_scheme = "implicit_picard" -algo.require_picard_convergence = 0 -algo.max_picard_iterations = 25 -algo.picard_iteration_tolerance = 0.0 #1.0e-12 +algo.evolve_scheme = "theta_implicit_em" +#algo.evolve_scheme = "semi_implicit_em" + +implicit_evolve.theta = 0.5 +implicit_evolve.max_particle_iterations = 21 +implicit_evolve.particle_tolerance = 1.0e-12 + +#implicit_evolve.nonlinear_solver = "picard" +#picard.verbose = true +#picard.max_iterations = 25 +#picard.relative_tolerance = 0.0 #1.0e-12 +#picard.absolute_tolerance = 0.0 #1.0e-24 +#picard.require_convergence = false + +implicit_evolve.nonlinear_solver = "newton" +newton.verbose = true +newton.max_iterations = 20 +newton.relative_tolerance = 1.0e-12 +newton.absolute_tolerance = 0.0 +newton.require_convergence = false + +gmres.verbose_int = 2 +gmres.max_iterations = 1000 +gmres.relative_tolerance = 1.0e-8 +gmres.absolute_tolerance = 0.0 + algo.particle_pusher = "boris" #algo.particle_pusher = "higuera" diff --git a/Regression/Checksum/benchmarks_json/ImplicitPicard_VandB_2d.json b/Regression/Checksum/benchmarks_json/ThetaImplicitJFNK_VandB_2d.json similarity index 100% rename from Regression/Checksum/benchmarks_json/ImplicitPicard_VandB_2d.json rename to Regression/Checksum/benchmarks_json/ThetaImplicitJFNK_VandB_2d.json diff --git a/Regression/Checksum/benchmarks_json/ImplicitPicard_1d.json b/Regression/Checksum/benchmarks_json/ThetaImplicitPicard_1d.json similarity index 100% rename from Regression/Checksum/benchmarks_json/ImplicitPicard_1d.json rename to Regression/Checksum/benchmarks_json/ThetaImplicitPicard_1d.json diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 6ee7ad3929d..e2d69be9c60 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -4666,7 +4666,7 @@ particleTypes = electrons outputFile = Point_of_contact_EB_rz_plt analysisRoutine = Examples/Tests/point_of_contact_EB/analysis.py -[ImplicitPicard_1d] +[ThetaImplicitPicard_1d] buildDir = . inputFile = Examples/Tests/Implicit/inputs_1d runtime_params = warpx.abort_on_warning_threshold=high @@ -4683,9 +4683,9 @@ doVis = 0 compareParticles = 1 analysisRoutine = Examples/Tests/Implicit/analysis_1d.py -[ImplicitPicard_VandB_2d] +[ThetaImplicitJFNK_VandB_2d] buildDir = . -inputFile = Examples/Tests/Implicit/inputs_vandb_2d +inputFile = Examples/Tests/Implicit/inputs_vandb_jfnk_2d runtime_params = warpx.abort_on_warning_threshold=high dim = 2 addToCompileString = @@ -4698,7 +4698,7 @@ numthreads = 1 compileTest = 0 doVis = 0 compareParticles = 1 -analysisRoutine = Examples/Tests/Implicit/analysis_vandb_2d.py +analysisRoutine = Examples/Tests/Implicit/analysis_vandb_jfnk_2d.py [SemiImplicitPicard_1d] buildDir = . diff --git a/Source/Evolve/CMakeLists.txt b/Source/Evolve/CMakeLists.txt index a84d6e1e42b..3ef0dae5e8e 100644 --- a/Source/Evolve/CMakeLists.txt +++ b/Source/Evolve/CMakeLists.txt @@ -4,6 +4,5 @@ foreach(D IN LISTS WarpX_DIMS) PRIVATE WarpXEvolve.cpp WarpXComputeDt.cpp - WarpXOneStepImplicitPicard.cpp ) endforeach() diff --git a/Source/Evolve/Make.package b/Source/Evolve/Make.package index a9a877475b9..275b8bfde5a 100644 --- a/Source/Evolve/Make.package +++ b/Source/Evolve/Make.package @@ -1,5 +1,4 @@ CEXE_sources += WarpXEvolve.cpp CEXE_sources += WarpXComputeDt.cpp -CEXE_sources += WarpXOneStepImplicitPicard.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/Evolve diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 881967a0479..4ac50483ad3 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -125,10 +125,8 @@ WarpX::Evolve (int numsteps) ExecutePythonCallback("particleinjection"); - // TODO: move out - if (evolve_scheme == EvolveScheme::ImplicitPicard || - evolve_scheme == EvolveScheme::SemiImplicitPicard) { - OneStep_ImplicitPicard(cur_time); + if (m_implicit_solver) { + m_implicit_solver->OneStep(cur_time, dt[0], step); } else if ( electromagnetic_solver_id == ElectromagneticSolverAlgo::None || electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC ) diff --git a/Source/Evolve/WarpXOneStepImplicitPicard.cpp b/Source/Evolve/WarpXOneStepImplicitPicard.cpp deleted file mode 100644 index 6df3cd0d4e6..00000000000 --- a/Source/Evolve/WarpXOneStepImplicitPicard.cpp +++ /dev/null @@ -1,423 +0,0 @@ -/* Copyright 2022 David Grote - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#include "WarpX.H" - -#include "BoundaryConditions/PML.H" -#include "Diagnostics/MultiDiagnostics.H" -#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" -#include "Evolve/WarpXDtType.H" -#include "Evolve/WarpXPushType.H" -#ifdef WARPX_USE_FFT -# ifdef WARPX_DIM_RZ -# include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" -# else -# include "FieldSolver/SpectralSolver/SpectralSolver.H" -# endif -#endif -#include "Parallelization/GuardCellManager.H" -#include "Particles/MultiParticleContainer.H" -#include "Particles/ParticleBoundaryBuffer.H" -#include "Python/callbacks.H" -#include "Utils/TextMsg.H" -#include "Utils/WarpXAlgorithmSelection.H" -#include "Utils/WarpXUtil.H" -#include "Utils/WarpXConst.H" -#include "Utils/WarpXProfilerWrapper.H" - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -void -WarpX::EvolveImplicitPicardInit (int lev) -{ - - if (lev == 0) { - // Add space to save the positions and velocities at the start of the time steps - for (auto const& pc : *mypc) { -#if (AMREX_SPACEDIM >= 2) - pc->AddRealComp("x_n"); -#endif -#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - pc->AddRealComp("y_n"); -#endif - pc->AddRealComp("z_n"); - pc->AddRealComp("ux_n"); - pc->AddRealComp("uy_n"); - pc->AddRealComp("uz_n"); - } - } - - // Initialize MultiFabs to hold the E and B fields at the start of the time steps - // Only one refinement level is supported - const int nlevs_max = maxLevel() + 1; - Efield_n.resize(nlevs_max); - Efield_save.resize(nlevs_max); - if (evolve_scheme == EvolveScheme::ImplicitPicard) { - Bfield_n.resize(nlevs_max); - Bfield_save.resize(nlevs_max); - } - - // The Efield_n and Bfield_n will hold the fields at the start of the time step. - // This is needed since in each iteration the fields are advanced from the values - // at the start of the step. - // The Efield_save and Bfield_save will hold the fields from the previous iteration, - // to check the change in the fields after the iterations to check for convergence. - // The Efiel_fp and Bfield_fp will hole the n+theta during the iterations and then - // advance to the n+1 time level after the iterations complete. - AllocInitMultiFabFromModel(Efield_n[lev][0], *Efield_fp[0][0], lev, "Efield_n[0]"); - AllocInitMultiFabFromModel(Efield_n[lev][1], *Efield_fp[0][1], lev, "Efield_n[1]"); - AllocInitMultiFabFromModel(Efield_n[lev][2], *Efield_fp[0][2], lev, "Efield_n[2]"); - AllocInitMultiFabFromModel(Efield_save[lev][0], *Efield_fp[0][0], lev, "Efield_save[0]"); - AllocInitMultiFabFromModel(Efield_save[lev][1], *Efield_fp[0][1], lev, "Efield_save[1]"); - AllocInitMultiFabFromModel(Efield_save[lev][2], *Efield_fp[0][2], lev, "Efield_save[2]"); - - if (evolve_scheme == EvolveScheme::ImplicitPicard) { - AllocInitMultiFabFromModel(Bfield_n[lev][0], *Bfield_fp[0][0], lev, "Bfield_n[0]"); - AllocInitMultiFabFromModel(Bfield_n[lev][1], *Bfield_fp[0][1], lev, "Bfield_n[1]"); - AllocInitMultiFabFromModel(Bfield_n[lev][2], *Bfield_fp[0][2], lev, "Bfield_n[2]"); - AllocInitMultiFabFromModel(Bfield_save[lev][0], *Bfield_fp[0][0], lev, "Bfield_save[0]"); - AllocInitMultiFabFromModel(Bfield_save[lev][1], *Bfield_fp[0][1], lev, "Bfield_save[1]"); - AllocInitMultiFabFromModel(Bfield_save[lev][2], *Bfield_fp[0][2], lev, "Bfield_save[2]"); - } - -} - -void -WarpX::OneStep_ImplicitPicard(amrex::Real cur_time) -{ - using namespace amrex::literals; - - // We have E^{n}. - // Particles have p^{n} and x^{n}. - // With full implicit, B^{n} - // With semi-implicit, B^{n-1/2} - - // Save the values at the start of the time step, - // copying particle data to x_n etc. - for (auto const& pc : *mypc) { - SaveParticlesAtImplicitStepStart (*pc, 0); - } - - // Save the fields at the start of the step - amrex::MultiFab::Copy(*Efield_n[0][0], *Efield_fp[0][0], 0, 0, ncomps, Efield_fp[0][0]->nGrowVect()); - amrex::MultiFab::Copy(*Efield_n[0][1], *Efield_fp[0][1], 0, 0, ncomps, Efield_fp[0][1]->nGrowVect()); - amrex::MultiFab::Copy(*Efield_n[0][2], *Efield_fp[0][2], 0, 0, ncomps, Efield_fp[0][2]->nGrowVect()); - - if (evolve_scheme == EvolveScheme::ImplicitPicard) { - amrex::MultiFab::Copy(*Bfield_n[0][0], *Bfield_fp[0][0], 0, 0, ncomps, Bfield_fp[0][0]->nGrowVect()); - amrex::MultiFab::Copy(*Bfield_n[0][1], *Bfield_fp[0][1], 0, 0, ncomps, Bfield_fp[0][1]->nGrowVect()); - amrex::MultiFab::Copy(*Bfield_n[0][2], *Bfield_fp[0][2], 0, 0, ncomps, Bfield_fp[0][2]->nGrowVect()); - } else if (evolve_scheme == EvolveScheme::SemiImplicitPicard) { - // This updates Bfield_fp so it holds the new B at n+1/2 - EvolveB(dt[0], DtType::Full); - // WarpX::sync_nodal_points is used to avoid instability - FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); - ApplyBfieldBoundary(0, PatchType::fine, DtType::Full); - } - - // Start the iterations - amrex::Real deltaE = 1._rt; - amrex::Real deltaB = 1._rt; - int iteration_count = 0; - while (iteration_count < max_picard_iterations && - (deltaE > picard_iteration_tolerance || deltaB > picard_iteration_tolerance)) { - iteration_count++; - - // Advance the particle positions by 1/2 dt, - // particle velocities by dt, then take average of old and new v, - // deposit currents, giving J at n+1/2 - // This uses Efield_fp and Bfield_fp, the field at n+1/2 from the previous iteration. - const bool skip_current = false; - const PushType push_type = PushType::Implicit; - PushParticlesandDeposit(cur_time, skip_current, push_type); - - SyncCurrentAndRho(); - - if (picard_iteration_tolerance > 0. || iteration_count == max_picard_iterations) { - // Save the E at n+1/2 from the previous iteration so that the change - // in this iteration can be calculated - amrex::MultiFab::Copy(*Efield_save[0][0], *Efield_fp[0][0], 0, 0, ncomps, 0); - amrex::MultiFab::Copy(*Efield_save[0][1], *Efield_fp[0][1], 0, 0, ncomps, 0); - amrex::MultiFab::Copy(*Efield_save[0][2], *Efield_fp[0][2], 0, 0, ncomps, 0); - } - - // Copy Efield_n into Efield_fp since EvolveE updates Efield_fp in place - amrex::MultiFab::Copy(*Efield_fp[0][0], *Efield_n[0][0], 0, 0, ncomps, Efield_n[0][0]->nGrowVect()); - amrex::MultiFab::Copy(*Efield_fp[0][1], *Efield_n[0][1], 0, 0, ncomps, Efield_n[0][1]->nGrowVect()); - amrex::MultiFab::Copy(*Efield_fp[0][2], *Efield_n[0][2], 0, 0, ncomps, Efield_n[0][2]->nGrowVect()); - - // Updates Efield_fp so it holds the new E at n+1/2 - EvolveE(0.5_rt*dt[0]); - // WarpX::sync_nodal_points is used to avoid instability - FillBoundaryE(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); - ApplyEfieldBoundary(0, PatchType::fine); - - if (evolve_scheme == EvolveScheme::ImplicitPicard) { - if (picard_iteration_tolerance > 0. || iteration_count == max_picard_iterations) { - // Save the B at n+1/2 from the previous iteration so that the change - // in this iteration can be calculated - amrex::MultiFab::Copy(*Bfield_save[0][0], *Bfield_fp[0][0], 0, 0, ncomps, 0); - amrex::MultiFab::Copy(*Bfield_save[0][1], *Bfield_fp[0][1], 0, 0, ncomps, 0); - amrex::MultiFab::Copy(*Bfield_save[0][2], *Bfield_fp[0][2], 0, 0, ncomps, 0); - } - - // Copy Bfield_n into Bfield_fp since EvolveB updates Bfield_fp in place - amrex::MultiFab::Copy(*Bfield_fp[0][0], *Bfield_n[0][0], 0, 0, ncomps, Bfield_n[0][0]->nGrowVect()); - amrex::MultiFab::Copy(*Bfield_fp[0][1], *Bfield_n[0][1], 0, 0, ncomps, Bfield_n[0][1]->nGrowVect()); - amrex::MultiFab::Copy(*Bfield_fp[0][2], *Bfield_n[0][2], 0, 0, ncomps, Bfield_n[0][2]->nGrowVect()); - - // This updates Bfield_fp so it holds the new B at n+1/2 - EvolveB(0.5_rt*dt[0], DtType::Full); - // WarpX::sync_nodal_points is used to avoid instability - FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); - ApplyBfieldBoundary(0, PatchType::fine, DtType::Full); - } - - // The B field update needs - if (num_mirrors>0){ - applyMirrors(cur_time); - // E : guard cells are NOT up-to-date from the mirrors - // B : guard cells are NOT up-to-date from the mirrors - } - - if (picard_iteration_tolerance > 0. || iteration_count == max_picard_iterations) { - // Calculate the change in E and B from this iteration - // deltaE = abs(Enew - Eold)/max(abs(Enew)) - Efield_save[0][0]->minus(*Efield_fp[0][0], 0, ncomps, 0); - Efield_save[0][1]->minus(*Efield_fp[0][1], 0, ncomps, 0); - Efield_save[0][2]->minus(*Efield_fp[0][2], 0, ncomps, 0); - const amrex::Real maxE0 = std::max(1._rt, Efield_fp[0][0]->norm0(0, 0)); - const amrex::Real maxE1 = std::max(1._rt, Efield_fp[0][1]->norm0(0, 0)); - const amrex::Real maxE2 = std::max(1._rt, Efield_fp[0][2]->norm0(0, 0)); - const amrex::Real deltaE0 = Efield_save[0][0]->norm0(0, 0)/maxE0; - const amrex::Real deltaE1 = Efield_save[0][1]->norm0(0, 0)/maxE1; - const amrex::Real deltaE2 = Efield_save[0][2]->norm0(0, 0)/maxE2; - deltaE = std::max(std::max(deltaE0, deltaE1), deltaE2); - if (evolve_scheme == EvolveScheme::ImplicitPicard) { - Bfield_save[0][0]->minus(*Bfield_fp[0][0], 0, ncomps, 0); - Bfield_save[0][1]->minus(*Bfield_fp[0][1], 0, ncomps, 0); - Bfield_save[0][2]->minus(*Bfield_fp[0][2], 0, ncomps, 0); - const amrex::Real maxB0 = std::max(1._rt, Bfield_fp[0][0]->norm0(0, 0)); - const amrex::Real maxB1 = std::max(1._rt, Bfield_fp[0][1]->norm0(0, 0)); - const amrex::Real maxB2 = std::max(1._rt, Bfield_fp[0][2]->norm0(0, 0)); - const amrex::Real deltaB0 = Bfield_save[0][0]->norm0(0, 0)/maxB0; - const amrex::Real deltaB1 = Bfield_save[0][1]->norm0(0, 0)/maxB1; - const amrex::Real deltaB2 = Bfield_save[0][2]->norm0(0, 0)/maxB2; - deltaB = std::max(std::max(deltaB0, deltaB1), deltaB2); - } else { - deltaB = 0.; - } - amrex::Print() << "Max delta " << iteration_count << " " << deltaE << " " << deltaB << "\n"; - } - - // Now, the particle positions and velocities and the Efield_fp and Bfield_fp hold - // the new values at n+1/2 - } - - amrex::Print() << "Picard iterations = " << iteration_count << ", Eerror = " << deltaE << ", Berror = " << deltaB << "\n"; - if (picard_iteration_tolerance > 0. && iteration_count == max_picard_iterations) { - std::stringstream convergenceMsg; - convergenceMsg << "The Picard implicit solver failed to converge after " << iteration_count << " iterations, with Eerror = " << deltaE << ", Berror = " << deltaB << " with a tolerance of " << picard_iteration_tolerance; - if (require_picard_convergence) { - WARPX_ABORT_WITH_MESSAGE(convergenceMsg.str()); - } else { - ablastr::warn_manager::WMRecordWarning("PicardSolver", convergenceMsg.str()); - } - } - - // Advance particles to step n+1 - for (auto const& pc : *mypc) { - FinishImplicitParticleUpdate(*pc, 0); - } - - // Advance fields to step n+1 - // WarpX::sync_nodal_points is used to avoid instability - FinishImplicitFieldUpdate(Efield_fp, Efield_n); - FillBoundaryE(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); - if (evolve_scheme == EvolveScheme::ImplicitPicard) { - FinishImplicitFieldUpdate(Bfield_fp, Bfield_n); - FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); - } - -} - -void -WarpX::SaveParticlesAtImplicitStepStart (WarpXParticleContainer& pc, int lev) -{ - -#ifdef AMREX_USE_OMP -#pragma omp parallel -#endif - { - - auto particle_comps = pc.getParticleComps(); - - for (WarpXParIter pti(pc, lev); pti.isValid(); ++pti) { - - const auto getPosition = GetParticlePosition(pti); - - auto& attribs = pti.GetAttribs(); - amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - -#if (AMREX_SPACEDIM >= 2) - amrex::ParticleReal* x_n = pti.GetAttribs(particle_comps["x_n"]).dataPtr(); -#endif -#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - amrex::ParticleReal* y_n = pti.GetAttribs(particle_comps["y_n"]).dataPtr(); -#endif - amrex::ParticleReal* z_n = pti.GetAttribs(particle_comps["z_n"]).dataPtr(); - amrex::ParticleReal* ux_n = pti.GetAttribs(particle_comps["ux_n"]).dataPtr(); - amrex::ParticleReal* uy_n = pti.GetAttribs(particle_comps["uy_n"]).dataPtr(); - amrex::ParticleReal* uz_n = pti.GetAttribs(particle_comps["uz_n"]).dataPtr(); - - const long np = pti.numParticles(); - - amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) - { - amrex::ParticleReal xp, yp, zp; - getPosition(ip, xp, yp, zp); - -#if (AMREX_SPACEDIM >= 2) - x_n[ip] = xp; -#endif -#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - y_n[ip] = yp; -#endif - z_n[ip] = zp; - - ux_n[ip] = ux[ip]; - uy_n[ip] = uy[ip]; - uz_n[ip] = uz[ip]; - - }); - - } - } -} - -void -WarpX::FinishImplicitParticleUpdate (WarpXParticleContainer& pc, int lev) -{ - using namespace amrex::literals; - -#ifdef AMREX_USE_OMP -#pragma omp parallel -#endif - { - - auto particle_comps = pc.getParticleComps(); - - for (WarpXParIter pti(pc, lev); pti.isValid(); ++pti) { - - const auto getPosition = GetParticlePosition(pti); - const auto setPosition = SetParticlePosition(pti); - - auto& attribs = pti.GetAttribs(); - amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - -#if (AMREX_SPACEDIM >= 2) - amrex::ParticleReal* x_n = pti.GetAttribs(particle_comps["x_n"]).dataPtr(); -#endif -#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - amrex::ParticleReal* y_n = pti.GetAttribs(particle_comps["y_n"]).dataPtr(); -#endif - amrex::ParticleReal* z_n = pti.GetAttribs(particle_comps["z_n"]).dataPtr(); - amrex::ParticleReal* ux_n = pti.GetAttribs(particle_comps["ux_n"]).dataPtr(); - amrex::ParticleReal* uy_n = pti.GetAttribs(particle_comps["uy_n"]).dataPtr(); - amrex::ParticleReal* uz_n = pti.GetAttribs(particle_comps["uz_n"]).dataPtr(); - - const long np = pti.numParticles(); - - amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) - { - amrex::ParticleReal xp, yp, zp; - getPosition(ip, xp, yp, zp); - -#if (AMREX_SPACEDIM >= 2) - xp = 2._rt*xp - x_n[ip]; -#endif -#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - yp = 2._rt*yp - y_n[ip]; -#endif - zp = 2._rt*zp - z_n[ip]; - - ux[ip] = 2._rt*ux[ip] - ux_n[ip]; - uy[ip] = 2._rt*uy[ip] - uy_n[ip]; - uz[ip] = 2._rt*uz[ip] - uz_n[ip]; - - setPosition(ip, xp, yp, zp); - }); - - } - } -} - -void -WarpX::FinishImplicitFieldUpdate(amrex::Vector, 3 > >& Field_fp, - amrex::Vector, 3 > >& Field_n) -{ - using namespace amrex::literals; - - for (int lev = 0; lev <= finest_level; ++lev) { - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( amrex::MFIter mfi(*Field_fp[lev][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - - amrex::Array4 const& Fx = Field_fp[lev][0]->array(mfi); - amrex::Array4 const& Fy = Field_fp[lev][1]->array(mfi); - amrex::Array4 const& Fz = Field_fp[lev][2]->array(mfi); - - amrex::Array4 const& Fx_n = Field_n[lev][0]->array(mfi); - amrex::Array4 const& Fy_n = Field_n[lev][1]->array(mfi); - amrex::Array4 const& Fz_n = Field_n[lev][2]->array(mfi); - - amrex::Box const tbx = mfi.tilebox(Field_fp[lev][0]->ixType().toIntVect()); - amrex::Box const tby = mfi.tilebox(Field_fp[lev][1]->ixType().toIntVect()); - amrex::Box const tbz = mfi.tilebox(Field_fp[lev][2]->ixType().toIntVect()); - - amrex::ParallelFor( - tbx, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) - { - Fx(i,j,k,n) = 2._rt*Fx(i,j,k,n) - Fx_n(i,j,k,n); - }, - tby, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) - { - Fy(i,j,k,n) = 2._rt*Fy(i,j,k,n) - Fy_n(i,j,k,n); - }, - tbz, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) - { - Fz(i,j,k,n) = 2._rt*Fz(i,j,k,n) - Fz_n(i,j,k,n); - }); - } - } -} diff --git a/Source/FieldSolver/CMakeLists.txt b/Source/FieldSolver/CMakeLists.txt index cf741aa6e01..859117eb214 100644 --- a/Source/FieldSolver/CMakeLists.txt +++ b/Source/FieldSolver/CMakeLists.txt @@ -11,6 +11,7 @@ endforeach() add_subdirectory(FiniteDifferenceSolver) add_subdirectory(MagnetostaticSolver) +add_subdirectory(ImplicitSolvers) if(WarpX_FFT) add_subdirectory(SpectralSolver) endif() diff --git a/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt b/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt new file mode 100644 index 00000000000..6e16f19084c --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/CMakeLists.txt @@ -0,0 +1,10 @@ +foreach(D IN LISTS WarpX_DIMS) + warpx_set_suffix_dims(SD ${D}) + target_sources(lib_${SD} + PRIVATE + SemiImplicitEM.cpp + ThetaImplicitEM.cpp + WarpXImplicitOps.cpp + WarpXSolverVec.cpp + ) +endforeach() diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H new file mode 100644 index 00000000000..88ad6a058fd --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H @@ -0,0 +1,145 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef Implicit_Solver_H_ +#define Implicit_Solver_H_ + +#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" +#include "NonlinearSolvers/NonlinearSolverLibrary.H" + +#include +#include + +/** + * \brief Base class for implicit time solvers. The base functions are those + * needed by an implicit solver to be used through WarpX and those needed + * to use the nonlinear solvers, such as Picard or Newton (i.e., JFNK). + */ + +class WarpX; +class ImplicitSolver +{ +public: + + ImplicitSolver() = default; + + virtual ~ImplicitSolver() = default; + + // Prohibit Move and Copy operations + ImplicitSolver(const ImplicitSolver&) = delete; + ImplicitSolver& operator=(const ImplicitSolver&) = delete; + ImplicitSolver(ImplicitSolver&&) = delete; + ImplicitSolver& operator=(ImplicitSolver&&) = delete; + + // + // the following routines are called by WarpX + // + + /** + * \brief Read user-provided parameters that control the implicit solver. + * Allocate internal arrays for intermediate field values needed by the solver. + */ + virtual void Define ( WarpX* a_WarpX ) = 0; + + [[nodiscard]] bool IsDefined () const { return m_is_defined; } + + virtual void PrintParameters () const = 0; + + void GetParticleSolverParams (int& a_max_particle_iter, + amrex::ParticleReal& a_particle_tol ) const + { + a_max_particle_iter = m_max_particle_iterations; + a_particle_tol = m_particle_tolerance; + } + + /** + * \brief Advance fields and particles by one time step using the specified implicit algorithm + */ + virtual void OneStep ( amrex::Real a_time, + amrex::Real a_dt, + int a_step ) = 0; + + // + // the following routines are called by the linear and nonlinear solvers + // + + /** + * \brief Computes the RHS of the equation corresponding to the specified implicit algorithm. + * The discrete equations corresponding to numerical integration of ODEs are often + * written in the form U = b + RHS(U), where U is the variable to be solved for (e.g., + * the solution at the next time step), b is a constant (i.e., the solution from the + * previous time step), and RHS(U) is the right-hand-side of the equation. Iterative + * solvers, such as Picard and Newton, and higher-order Runge-Kutta methods, need to + * compute RHS(U) multiple times for different values of U. Thus, a routine that + * returns this value is needed. + * e.g., Ebar - E^n = cvac^2*0.5*dt*(curl(Bbar) - mu0*Jbar(Ebar,Bbar)) + * Here, U = Ebar, b = E^n, and the expression on the right is RHS(U). + */ + virtual void ComputeRHS ( WarpXSolverVec& a_RHS, + const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt, + int a_nl_iter, + bool a_from_jacobian ) = 0; + +protected: + + /** + * \brief Pointer back to main WarpX class + */ + WarpX* m_WarpX; + + bool m_is_defined = false; + + /** + * \brief Nonlinear solver type and object + */ + NonlinearSolverType m_nlsolver_type; + std::unique_ptr> m_nlsolver; + + /** + * \brief tolerance used by the iterative method used to obtain a self-consistent + * update of the particle positions and velocities for given E and B on the grid + */ + amrex::ParticleReal m_particle_tolerance = 1.0e-10; + + /** + * \brief maximum iterations for the iterative method used to obtain a self-consistent + * update of the particle positions and velocities for given E and B on the grid + */ + int m_max_particle_iterations = 21; + + /** + * \brief parse nonlinear solver parameters (if one is used) + */ + void parseNonlinearSolverParams( const amrex::ParmParse& pp ) + { + + std::string nlsolver_type_str; + pp.get("nonlinear_solver", nlsolver_type_str); + + if (nlsolver_type_str=="picard") { + m_nlsolver_type = NonlinearSolverType::Picard; + m_nlsolver = std::make_unique>(); + m_max_particle_iterations = 1; + m_particle_tolerance = 0.0; + } + else if (nlsolver_type_str=="newton") { + m_nlsolver_type = NonlinearSolverType::Newton; + m_nlsolver = std::make_unique>(); + pp.query("max_particle_iterations", m_max_particle_iterations); + pp.query("particle_tolerance", m_particle_tolerance); + } + else { + WARPX_ABORT_WITH_MESSAGE( + "invalid nonlinear_solver specified. Valid options are picard and newton."); + } + + } + +}; + +#endif diff --git a/Source/FieldSolver/ImplicitSolvers/ImplicitSolverLibrary.H b/Source/FieldSolver/ImplicitSolvers/ImplicitSolverLibrary.H new file mode 100644 index 00000000000..423957ef061 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/ImplicitSolverLibrary.H @@ -0,0 +1,13 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef IMPLICIT_SOLVER_LIBRARY_H_ +#define IMPLICIT_SOLVER_LIBRARY_H_ + +#include "SemiImplicitEM.H" // IWYU pragma: export +#include "ThetaImplicitEM.H" // IWYU pragma: export + +#endif diff --git a/Source/FieldSolver/ImplicitSolvers/Make.package b/Source/FieldSolver/ImplicitSolvers/Make.package new file mode 100644 index 00000000000..a4543f94dd3 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/Make.package @@ -0,0 +1,6 @@ +CEXE_sources += SemiImplicitEM.cpp +CEXE_sources += ThetaImplicitEM.cpp +CEXE_sources += WarpXImplicitOps.cpp +CEXE_sources += WarpXSolverVec.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/ImplicitSolvers diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H new file mode 100644 index 00000000000..6e3e5db2c74 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.H @@ -0,0 +1,80 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef SEMI_IMPLICIT_EM_H_ +#define SEMI_IMPLICIT_EM_H_ + +#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" + +#include +#include +#include + +#include "ImplicitSolver.H" + +/** @file + * Semi-implicit electromagnetic time solver class. The electric field and the + * particles are implicitly coupled in this algorithm, but the magnetic field + * is advanced in the standard explicit leap-frog manner (hence semi-implicit). + * + * The time stencil is + * Eg^{n+1} = Eg^n + c^2*dt*( curlBg^{n+1/2} - mu0*Jg^{n+1/2} ) + * Bg^{n+3/2} = Bg^{n+1/2} - dt*curlEg^{n+1} + * xp^{n+1} = xp^n + dt*up^{n+1/2}/(0.5(gammap^n + gammap^{n+1})) + * up^{n+1} = up^n + dt*qp/mp*(Ep^{n+1/2} + up^{n+1/2}/gammap^{n+1/2} x Bp^{n+1/2}) + * where f^{n+1/2} = (f^{n+1} + f^n)/2.0, for all but Bg, which lives at half steps + * + * This algorithm is approximately energy conserving. It is exactly energy conserving + * using a non-standard definition for the magnetic field energy. The advantage of + * this method over the exactly energy-conserving theta-implicit EM method is that + * light wave dispersion is captured much better. However, the CFL condition for light + * waves has to be satisifed for numerical stability (and for the modified definition + * of the magnetic field energy to be well-posed). + * + * See G. Chen, L. Chacon, L. Yin, B.J. Albright, D.J. Stark, R.F. Bird, + * "A semi-implicit energy- and charge-conserving particle-in-cell algorithm for the + * relativistic Vlasov-Maxwell equations.", JCP 407 (2020). + */ + +class SemiImplicitEM : public ImplicitSolver +{ +public: + + SemiImplicitEM() = default; + + ~SemiImplicitEM() override = default; + + // Prohibit Move and Copy operations + SemiImplicitEM(const SemiImplicitEM&) = delete; + SemiImplicitEM& operator=(const SemiImplicitEM&) = delete; + SemiImplicitEM(SemiImplicitEM&&) = delete; + SemiImplicitEM& operator=(SemiImplicitEM&&) = delete; + + void Define ( WarpX* a_WarpX ) override; + + void PrintParameters () const override; + + void OneStep ( amrex::Real a_time, + amrex::Real a_dt, + int a_step ) override; + + void ComputeRHS ( WarpXSolverVec& a_RHS, + const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt, + int a_nl_iter, + bool a_from_jacobian ) override; + +private: + + /** + * \brief Solver vectors for E and Eold + */ + WarpXSolverVec m_E, m_Eold; + +}; + +#endif diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp new file mode 100644 index 00000000000..897bd5c07f3 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp @@ -0,0 +1,117 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "SemiImplicitEM.H" +#include "WarpX.H" + +using namespace warpx::fields; +using namespace amrex::literals; + +void SemiImplicitEM::Define ( WarpX* a_WarpX ) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !m_is_defined, + "SemiImplicitEM object is already defined!"); + + // Retain a pointer back to main WarpX class + m_WarpX = a_WarpX; + + // Define E and Eold vectors + m_E.Define( m_WarpX->getMultiLevelField(FieldType::Efield_fp) ); + m_Eold.Define( m_WarpX->getMultiLevelField(FieldType::Efield_fp) ); + + // Need to define the WarpXSolverVec owned dot_mask to do dot + // product correctly for linear and nonlinear solvers + const amrex::Vector& Geom = m_WarpX->Geom(); + m_E.SetDotMask(Geom); + + // Parse implicit solver parameters + const amrex::ParmParse pp("implicit_evolve"); + parseNonlinearSolverParams( pp ); + + // Define the nonlinear solver + m_nlsolver->Define(m_E, this); + m_is_defined = true; + +} + +void SemiImplicitEM::PrintParameters () const +{ + if (!m_WarpX->Verbose()) { return; } + amrex::Print() << std::endl; + amrex::Print() << "-----------------------------------------------------------" << std::endl; + amrex::Print() << "----------- SEMI IMPLICIT EM SOLVER PARAMETERS ------------" << std::endl; + amrex::Print() << "-----------------------------------------------------------" << std::endl; + amrex::Print() << "max particle iterations: " << m_max_particle_iterations << std::endl; + amrex::Print() << "particle tolerance: " << m_particle_tolerance << std::endl; + if (m_nlsolver_type==NonlinearSolverType::Picard) { + amrex::Print() << "Nonlinear solver type: Picard" << std::endl; + } + else if (m_nlsolver_type==NonlinearSolverType::Newton) { + amrex::Print() << "Nonlinear solver type: Newton" << std::endl; + } + m_nlsolver->PrintParams(); + amrex::Print() << "-----------------------------------------------------------" << std::endl; + amrex::Print() << std::endl; +} + +void SemiImplicitEM::OneStep ( amrex::Real a_time, + amrex::Real a_dt, + int a_step ) +{ + amrex::ignore_unused(a_step); + + // Fields have E^{n}, B^{n-1/2} + // Particles have p^{n} and x^{n}. + + // Save the values at the start of the time step, + m_WarpX->SaveParticlesAtImplicitStepStart ( ); + + // Save the fields at the start of the step + m_Eold.Copy( m_WarpX->getMultiLevelField(FieldType::Efield_fp) ); + m_E.Copy(m_Eold); // initial guess for E + + // Compute Bfield at time n+1/2 + m_WarpX->EvolveB(a_dt, DtType::Full); + m_WarpX->ApplyMagneticFieldBCs(); + + const amrex::Real half_time = a_time + 0.5_rt*a_dt; + + // Solve nonlinear system for E at t_{n+1/2} + // Particles will be advanced to t_{n+1/2} + m_nlsolver->Solve( m_E, m_Eold, half_time, a_dt ); + + // Update WarpX owned Efield_fp to t_{n+1/2} + m_WarpX->SetElectricFieldAndApplyBCs( m_E ); + + // Advance particles from time n+1/2 to time n+1 + m_WarpX->FinishImplicitParticleUpdate(); + + // Advance E fields from time n+1/2 to time n+1 + // Eg^{n+1} = 2.0*E_g^{n+1/2} - E_g^n + m_E.linComb( 2._rt, m_E, -1._rt, m_Eold ); + m_WarpX->SetElectricFieldAndApplyBCs( m_E ); + +} + +void SemiImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS, + const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt, + int a_nl_iter, + bool a_from_jacobian ) +{ + // update WarpX-owned Efield_fp using current state of E from + // the nonlinear solver at time n+theta + m_WarpX->SetElectricFieldAndApplyBCs( a_E ); + + // Self consistently update particle positions and velocities using the + // current state of the fields E and B. Deposit current density at time n+1/2. + m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian ); + + // RHS = cvac^2*0.5*dt*( curl(B^{n+1/2}) - mu0*J^{n+1/2} ) + m_WarpX->ImplicitComputeRHSE(0.5_rt*a_dt, a_RHS); +} diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H new file mode 100644 index 00000000000..009c2c7e546 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H @@ -0,0 +1,124 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef THETA_IMPLICIT_EM_H_ +#define THETA_IMPLICIT_EM_H_ + +#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.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. + * + * The time stencil is + * Eg^{n+1} = Eg^n + c^2*dt*( curlBg^{n+theta} - mu0*Jg^{n+1/2} ) + * Bg^{n+1} = Bg^n - dt*curlEg^{n+theta} + * xp^{n+1} = xp^n + dt*up^{n+1/2}/(0.5*(gammap^n + gammap^{n+1})) + * up^{n+1} = up^n + dt*qp/mp*(Ep^{n+theta} + up^{n+1/2}/gammap^{n+1/2} x Bp^{n+theta}) + * where f^{n+theta} = (1.0-theta)*f^{n} + theta*f^{n+1} with 0.5 <= theta <= 1.0 + * + * The user-specified time-biasing parameter theta used for the fields on the RHS is bound + * between 0.5 and 1.0. The algorithm is exactly energy conserving for theta = 0.5. + * Signifcant damping of high-k modes will occur as theta approaches 1.0. The algorithm is + * numerially stable for any time step. I.e., the CFL condition for light waves does not + * have to be satisifed and the time step is not limited by the plasma period. However, how + * efficiently the algorithm can use large time steps depends strongly on the nonlinear solver. + * Furthermore, the time step should always be such that particles do not travel outside the + * ghost region of the box they live in, which is an MPI-related limitation. The time step + * is always limited by the need to resolve the appropriate physics. + * + * See S. Markidis, G. Lapenta, "The energy conserving particle-in-cell method." JCP 230 (2011). + * + * See G. Chen, L. Chacon, D.C. Barnes, "An energy- and charge-conserving, implicit, + * elctrostatic particle-in-cell algorithm." JCP 230 (2011). + * + * See J.R. Angus, A. Link, A. Friedman, D. Ghosh, J. D. Johnson, "On numerical energy + * conservation for an implicit particle-in-cell method coupled with a binary Monte-Carlo + * algorithm for Coulomb collisions.", JCP 456 (2022). + * + * See J.R. Angus, W. Farmer, A. Friedman, D. Ghosh, D. Grote, D. Larson, A. Link, "An + * implicit particle code with exact energy and charge conservation for electromagnetic studies + * of dense plasmas.", JCP 491 (2023). + */ + +class ThetaImplicitEM : public ImplicitSolver +{ +public: + + ThetaImplicitEM() = default; + + ~ThetaImplicitEM() override = default; + + // Prohibit Move and Copy operations + ThetaImplicitEM(const ThetaImplicitEM&) = delete; + ThetaImplicitEM& operator=(const ThetaImplicitEM&) = delete; + ThetaImplicitEM(ThetaImplicitEM&&) = delete; + ThetaImplicitEM& operator=(ThetaImplicitEM&&) = delete; + + void Define ( WarpX* a_WarpX ) override; + + void PrintParameters () const override; + + void OneStep ( amrex::Real a_time, + amrex::Real a_dt, + int a_step ) override; + + void ComputeRHS ( WarpXSolverVec& a_RHS, + const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt, + int a_nl_iter, + bool a_from_jacobian ) override; + + [[nodiscard]] amrex::Real theta () const { return m_theta; } + +private: + + /** + * \brief Time-biasing parameter for fields used on RHS to advance system + */ + amrex::Real m_theta = 0.5; + + /** + * \brief Solver vectors to be used in the nonlinear solver to solve for the + * electric field E. The main logic for determining which variables should be + * WarpXSolverVec type is that it must have the same size and have the same + * centering of the data as the variable being solved for, which is E here. + * For example, if using a Yee grid then a container for curlB could be a + * WarpXSovlerVec, but magnetic field B should not be. + */ + WarpXSolverVec m_E, m_Eold; + + /** + * \brief B is a derived variable from E. Need to save Bold to update B during + * the iterative nonlinear solve for E. Bold is owned here, but only used by WarpX. + * It is not used directly by the nonlinear solver, nor is it the same size as the + * solver vector (size E), and so it should not be WarpXSolverVec type. + */ + amrex::Vector, 3 > > m_Bold; + + /** + * \brief Update the E and B fields owned by WarpX + */ + void UpdateWarpXFields ( const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt ); + + /** + * \brief Nonlinear solver is for the time-centered values of E. After + * the solver, need to use m_E and m_Eold to compute E^{n+1} + */ + void FinishFieldUpdate ( amrex::Real a_new_time ); + +}; + +#endif diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp new file mode 100644 index 00000000000..026c509c3ba --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp @@ -0,0 +1,171 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "FieldSolver/Fields.H" +#include "ThetaImplicitEM.H" +#include "WarpX.H" + +using namespace warpx::fields; +using namespace amrex::literals; + +void ThetaImplicitEM::Define ( WarpX* const a_WarpX ) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !m_is_defined, + "ThetaImplicitEM object is already defined!"); + + // Retain a pointer back to main WarpX class + m_WarpX = a_WarpX; + + // Define E and Eold vectors + m_E.Define( m_WarpX->getMultiLevelField(FieldType::Efield_fp) ); + m_Eold.Define( m_WarpX->getMultiLevelField(FieldType::Efield_fp) ); + + // Need to define the WarpXSolverVec owned dot_mask to do dot + // product correctly for linear and nonlinear solvers + const amrex::Vector& Geom = m_WarpX->Geom(); + m_E.SetDotMask(Geom); + + // Define Bold MultiFab + const int num_levels = 1; + m_Bold.resize(num_levels); // size is number of levels + for (int lev = 0; lev < num_levels; ++lev) { + for (int n=0; n<3; n++) { + const amrex::MultiFab& Bfp = m_WarpX->getField( FieldType::Bfield_fp,lev,n); + m_Bold[lev][n] = std::make_unique( Bfp.boxArray(), + Bfp.DistributionMap(), + Bfp.nComp(), + Bfp.nGrowVect() ); + } + } + + // Parse theta-implicit solver specific parameters + const amrex::ParmParse pp("implicit_evolve"); + pp.query("theta", m_theta); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + m_theta>=0.5 && m_theta<=1.0, + "theta parameter for theta implicit time solver must be between 0.5 and 1.0"); + + // Parse nonlinear solver parameters + parseNonlinearSolverParams( pp ); + + // Define the nonlinear solver + m_nlsolver->Define(m_E, this); + m_is_defined = true; + +} + +void ThetaImplicitEM::PrintParameters () const +{ + if (!m_WarpX->Verbose()) { return; } + amrex::Print() << std::endl; + amrex::Print() << "-----------------------------------------------------------" << std::endl; + amrex::Print() << "----------- THETA IMPLICIT EM SOLVER PARAMETERS -----------" << std::endl; + amrex::Print() << "-----------------------------------------------------------" << std::endl; + amrex::Print() << "Time-bias parameter theta: " << m_theta << std::endl; + amrex::Print() << "max particle iterations: " << m_max_particle_iterations << std::endl; + amrex::Print() << "particle tolerance: " << m_particle_tolerance << std::endl; + if (m_nlsolver_type==NonlinearSolverType::Picard) { + amrex::Print() << "Nonlinear solver type: Picard" << std::endl; + } + else if (m_nlsolver_type==NonlinearSolverType::Newton) { + amrex::Print() << "Nonlinear solver type: Newton" << std::endl; + } + m_nlsolver->PrintParams(); + amrex::Print() << "-----------------------------------------------------------" << std::endl; + amrex::Print() << std::endl; +} + +void ThetaImplicitEM::OneStep ( const amrex::Real a_time, + const amrex::Real a_dt, + const int a_step ) +{ + amrex::ignore_unused(a_step); + + // Fields have E^{n} and B^{n} + // Particles have p^{n} and x^{n}. + + // Save the values at the start of the time step, + m_WarpX->SaveParticlesAtImplicitStepStart ( ); + + // Save the fields at the start of the step + m_Eold.Copy( m_WarpX->getMultiLevelField(FieldType::Efield_fp) ); + m_E.Copy(m_Eold); // initial guess for E + + const int num_levels = static_cast(m_Bold.size()); + for (int lev = 0; lev < num_levels; ++lev) { + for (int n=0; n<3; n++) { + const amrex::MultiFab& Bfp = m_WarpX->getField(FieldType::Bfield_fp,lev,n); + amrex::MultiFab& Bold = *m_Bold[lev][n]; + amrex::MultiFab::Copy(Bold, Bfp, 0, 0, 1, Bold.nGrowVect()); + } + } + + const amrex::Real theta_time = a_time + m_theta*a_dt; + + // Solve nonlinear system for E at t_{n+theta} + // Particles will be advanced to t_{n+1/2} + m_nlsolver->Solve( m_E, m_Eold, theta_time, a_dt ); + + // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta} + UpdateWarpXFields( m_E, theta_time, a_dt ); + + // Advance particles from time n+1/2 to time n+1 + m_WarpX->FinishImplicitParticleUpdate(); + + // Advance E and B fields from time n+theta to time n+1 + const amrex::Real new_time = a_time + a_dt; + FinishFieldUpdate( new_time ); + +} + +void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS, + const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt, + int a_nl_iter, + bool a_from_jacobian ) +{ + // update WarpX-owned Efield_fp and Bfield_fp using current state of E from + // the nonlinear solver at time n+theta + UpdateWarpXFields( a_E, a_time, a_dt ); + + // Self consistently update particle positions and velocities using the + // current state of the fields E and B. Deposit current density at time n+1/2. + m_WarpX->ImplicitPreRHSOp( a_time, a_dt, a_nl_iter, a_from_jacobian ); + + // RHS = cvac^2*m_theta*dt*( curl(B^{n+theta}) - mu0*J^{n+1/2} ) + m_WarpX->ImplicitComputeRHSE(m_theta*a_dt, a_RHS); +} + +void ThetaImplicitEM::UpdateWarpXFields ( const WarpXSolverVec& a_E, + amrex::Real a_time, + amrex::Real a_dt ) +{ + amrex::ignore_unused(a_time); + + // Update Efield_fp owned by WarpX + m_WarpX->SetElectricFieldAndApplyBCs( a_E ); + + // Update Bfield_fp owned by WarpX + m_WarpX->UpdateMagneticFieldAndApplyBCs( m_Bold, m_theta*a_dt ); + +} + +void ThetaImplicitEM::FinishFieldUpdate ( amrex::Real a_new_time ) +{ + amrex::ignore_unused(a_new_time); + + // Eg^{n+1} = (1/theta)*E_g^{n+theta} + (1-1/theta)*E_g^n + // Bg^{n+1} = (1/theta)*B_g^{n+theta} + (1-1/theta)*B_g^n + + const amrex::Real c0 = 1._rt/m_theta; + const amrex::Real c1 = 1._rt - c0; + m_E.linComb( c0, m_E, c1, m_Eold ); + m_WarpX->SetElectricFieldAndApplyBCs( m_E ); + m_WarpX->FinishMagneticFieldAndApplyBCs( m_Bold, m_theta ); + +} diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp b/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp new file mode 100644 index 00000000000..7b95abade3e --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp @@ -0,0 +1,352 @@ +/* Copyright 2022 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "WarpX.H" + +#include "BoundaryConditions/PML.H" +#include "Diagnostics/MultiDiagnostics.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" +#include "Evolve/WarpXDtType.H" +#include "Evolve/WarpXPushType.H" +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" +#ifdef WARPX_USE_PSATD +# ifdef WARPX_DIM_RZ +# include "FieldSolver/SpectralSolver/SpectralSolverRZ.H" +# else +# include "FieldSolver/SpectralSolver/SpectralSolver.H" +# endif +#endif +#include "Parallelization/GuardCellManager.H" +#include "Particles/MultiParticleContainer.H" +#include "Particles/ParticleBoundaryBuffer.H" +#include "Python/callbacks.H" +#include "Utils/TextMsg.H" +#include "Utils/WarpXAlgorithmSelection.H" +#include "Utils/WarpXUtil.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +void +WarpX::ImplicitPreRHSOp ( amrex::Real a_cur_time, + amrex::Real a_full_dt, + int a_nl_iter, + bool a_from_jacobian ) +{ + using namespace amrex::literals; + amrex::ignore_unused( a_full_dt, a_nl_iter, a_from_jacobian ); + + // Advance the particle positions by 1/2 dt, + // particle velocities by dt, then take average of old and new v, + // deposit currents, giving J at n+1/2 + // This uses Efield_fp and Bfield_fp, the field at n+1/2 from the previous iteration. + const bool skip_current = false; + const PushType push_type = PushType::Implicit; + PushParticlesandDeposit(a_cur_time, skip_current, push_type); + + SyncCurrentAndRho(); + +} + +void +WarpX::SetElectricFieldAndApplyBCs ( const WarpXSolverVec& a_E ) +{ + const amrex::Vector, 3 > >& Evec = a_E.getVec(); + amrex::MultiFab::Copy(*Efield_fp[0][0], *Evec[0][0], 0, 0, ncomps, Evec[0][0]->nGrowVect()); + amrex::MultiFab::Copy(*Efield_fp[0][1], *Evec[0][1], 0, 0, ncomps, Evec[0][1]->nGrowVect()); + amrex::MultiFab::Copy(*Efield_fp[0][2], *Evec[0][2], 0, 0, ncomps, Evec[0][2]->nGrowVect()); + FillBoundaryE(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); + ApplyEfieldBoundary(0, PatchType::fine); +} + +void +WarpX::UpdateMagneticFieldAndApplyBCs( const amrex::Vector, 3 > >& a_Bn, + amrex::Real a_thetadt ) +{ + amrex::MultiFab::Copy(*Bfield_fp[0][0], *a_Bn[0][0], 0, 0, ncomps, a_Bn[0][0]->nGrowVect()); + amrex::MultiFab::Copy(*Bfield_fp[0][1], *a_Bn[0][1], 0, 0, ncomps, a_Bn[0][1]->nGrowVect()); + amrex::MultiFab::Copy(*Bfield_fp[0][2], *a_Bn[0][2], 0, 0, ncomps, a_Bn[0][2]->nGrowVect()); + EvolveB(a_thetadt, DtType::Full); + ApplyMagneticFieldBCs(); +} + +void +WarpX::FinishMagneticFieldAndApplyBCs( const amrex::Vector, 3 > >& a_Bn, + amrex::Real a_theta ) +{ + FinishImplicitField(Bfield_fp, a_Bn, a_theta); + ApplyMagneticFieldBCs(); +} + +void +WarpX::ApplyMagneticFieldBCs() +{ + FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points); + ApplyBfieldBoundary(0, PatchType::fine, DtType::Full); +} + +void +WarpX::SaveParticlesAtImplicitStepStart ( ) +{ + // The implicit advance routines require the particle velocity + // and position values at the beginning of the step to compute the + // time-centered position and velocity needed for the implicit stencil. + // Thus, we need to save this information. + + for (auto const& pc : *mypc) { + + for (int lev = 0; lev <= finest_level; ++lev) { +#ifdef AMREX_USE_OMP +#pragma omp parallel +#endif + { + + auto particle_comps = pc->getParticleComps(); + + for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) { + + const auto getPosition = GetParticlePosition(pti); + + auto& attribs = pti.GetAttribs(); + amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + +#if (AMREX_SPACEDIM >= 2) + amrex::ParticleReal* x_n = pti.GetAttribs(particle_comps["x_n"]).dataPtr(); +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + amrex::ParticleReal* y_n = pti.GetAttribs(particle_comps["y_n"]).dataPtr(); +#endif + amrex::ParticleReal* z_n = pti.GetAttribs(particle_comps["z_n"]).dataPtr(); + amrex::ParticleReal* ux_n = pti.GetAttribs(particle_comps["ux_n"]).dataPtr(); + amrex::ParticleReal* uy_n = pti.GetAttribs(particle_comps["uy_n"]).dataPtr(); + amrex::ParticleReal* uz_n = pti.GetAttribs(particle_comps["uz_n"]).dataPtr(); + + const long np = pti.numParticles(); + + amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) + { + amrex::ParticleReal xp, yp, zp; + getPosition(ip, xp, yp, zp); + +#if (AMREX_SPACEDIM >= 2) + x_n[ip] = xp; +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + y_n[ip] = yp; +#endif + z_n[ip] = zp; + + ux_n[ip] = ux[ip]; + uy_n[ip] = uy[ip]; + uz_n[ip] = uz[ip]; + + }); + + } + } + + } + + } + +} + +void +WarpX::FinishImplicitParticleUpdate () +{ + using namespace amrex::literals; + + // The implicit advance routines use the time-centered position and + // momentum to advance the system in time. Thus, at the end of the + // step we need to transform the particle postion and momentum from + // time n+1/2 to time n+1. This is done here. + + for (auto const& pc : *mypc) { + + for (int lev = 0; lev <= finest_level; ++lev) { +#ifdef AMREX_USE_OMP +#pragma omp parallel +#endif + { + + auto particle_comps = pc->getParticleComps(); + + for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) { + + const auto getPosition = GetParticlePosition(pti); + const auto setPosition = SetParticlePosition(pti); + + auto& attribs = pti.GetAttribs(); + amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + +#if (AMREX_SPACEDIM >= 2) + amrex::ParticleReal* x_n = pti.GetAttribs(particle_comps["x_n"]).dataPtr(); +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + amrex::ParticleReal* y_n = pti.GetAttribs(particle_comps["y_n"]).dataPtr(); +#endif + amrex::ParticleReal* z_n = pti.GetAttribs(particle_comps["z_n"]).dataPtr(); + amrex::ParticleReal* ux_n = pti.GetAttribs(particle_comps["ux_n"]).dataPtr(); + amrex::ParticleReal* uy_n = pti.GetAttribs(particle_comps["uy_n"]).dataPtr(); + amrex::ParticleReal* uz_n = pti.GetAttribs(particle_comps["uz_n"]).dataPtr(); + + const long np = pti.numParticles(); + + amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) + { + amrex::ParticleReal xp, yp, zp; + getPosition(ip, xp, yp, zp); + +#if (AMREX_SPACEDIM >= 2) + xp = 2._rt*xp - x_n[ip]; +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + yp = 2._rt*yp - y_n[ip]; +#endif + zp = 2._rt*zp - z_n[ip]; + + ux[ip] = 2._rt*ux[ip] - ux_n[ip]; + uy[ip] = 2._rt*uy[ip] - uy_n[ip]; + uz[ip] = 2._rt*uz[ip] - uz_n[ip]; + + setPosition(ip, xp, yp, zp); + }); + + } + } + + } + + } + +} + +void +WarpX::FinishImplicitField( amrex::Vector, 3 > >& Field_fp, + const amrex::Vector, 3 > >& Field_n, + amrex::Real theta ) +{ + using namespace amrex::literals; + + // The implicit field advance routines use the fields at time n+theta + // with 0.5 <= theta <= 1.0. Thus, at the end of the step we need to + // transform the fields from time n+theta to time n+1. This is done here. + + for (int lev = 0; lev <= finest_level; ++lev) { + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( amrex::MFIter mfi(*Field_fp[lev][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + + amrex::Array4 const& Fx = Field_fp[lev][0]->array(mfi); + amrex::Array4 const& Fy = Field_fp[lev][1]->array(mfi); + amrex::Array4 const& Fz = Field_fp[lev][2]->array(mfi); + + amrex::Array4 const& Fx_n = Field_n[lev][0]->array(mfi); + amrex::Array4 const& Fy_n = Field_n[lev][1]->array(mfi); + amrex::Array4 const& Fz_n = Field_n[lev][2]->array(mfi); + + amrex::Box const& tbx = mfi.tilebox(Field_n[lev][0]->ixType().toIntVect()); + amrex::Box const& tby = mfi.tilebox(Field_n[lev][1]->ixType().toIntVect()); + amrex::Box const& tbz = mfi.tilebox(Field_n[lev][2]->ixType().toIntVect()); + + const amrex::Real c0 = 1._rt/theta; + const amrex::Real c1 = 1._rt - c0; + + amrex::ParallelFor( + tbx, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + Fx(i,j,k,n) = c0*Fx(i,j,k,n) + c1*Fx_n(i,j,k,n); + }, + tby, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + Fy(i,j,k,n) = c0*Fy(i,j,k,n) + c1*Fy_n(i,j,k,n); + }, + tbz, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + Fz(i,j,k,n) = c0*Fz(i,j,k,n) + c1*Fz_n(i,j,k,n); + }); + } + } +} + +void +WarpX::ImplicitComputeRHSE (amrex::Real a_dt, WarpXSolverVec& a_Erhs_vec) +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + ImplicitComputeRHSE(lev, a_dt, a_Erhs_vec); + } +} + +void +WarpX::ImplicitComputeRHSE (int lev, amrex::Real a_dt, WarpXSolverVec& a_Erhs_vec) +{ + WARPX_PROFILE("WarpX::ImplicitComputeRHSE()"); + ImplicitComputeRHSE(lev, PatchType::fine, a_dt, a_Erhs_vec); + if (lev > 0) + { + ImplicitComputeRHSE(lev, PatchType::coarse, a_dt, a_Erhs_vec); + } +} + +void +WarpX::ImplicitComputeRHSE (int lev, PatchType patch_type, amrex::Real a_dt, WarpXSolverVec& a_Erhs_vec) +{ + // set RHS to zero value + a_Erhs_vec.getVec()[lev][0]->setVal(0.0); + a_Erhs_vec.getVec()[lev][1]->setVal(0.0); + a_Erhs_vec.getVec()[lev][2]->setVal(0.0); + + // Compute Efield_rhs in regular cells by calling EvolveE. Because + // a_Erhs_vec is set to zero above, calling EvolveE below results in + // a_Erhs_vec storing only the RHS of the update equation. I.e., + // c^2*dt*(curl(B^{n+theta} - mu0*J^{n+1/2}) + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveE( a_Erhs_vec.getVec()[lev], Bfield_fp[lev], + current_fp[lev], m_edge_lengths[lev], + m_face_areas[lev], ECTRhofield[lev], + F_fp[lev], lev, a_dt ); + } else { + m_fdtd_solver_cp[lev]->EvolveE( a_Erhs_vec.getVec()[lev], Bfield_cp[lev], + current_cp[lev], m_edge_lengths[lev], + m_face_areas[lev], ECTRhofield[lev], + F_cp[lev], lev, a_dt ); + } + + // Compute Efield_rhs in PML cells by calling EvolveEPML + if (do_pml && pml[lev]->ok()) { + amrex::Abort("PML not yet implemented with implicit solvers."); + } + +} diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H new file mode 100644 index 00000000000..89a0b82b700 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.H @@ -0,0 +1,233 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WarpXSolverVec_H_ +#define WarpXSolverVec_H_ + +#include "Utils/TextMsg.H" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +/** + * \brief This is a wrapper class around a Vector of array of pointers to MultiFabs that + * contains basic math operators and functionality needed to interact with nonlinear + * solvers in WarpX and linear solvers in AMReX, such as GMRES. + * + * The size of the Vector is the number of amr levels. Hardcoded for 1 right now. The + * size of the array is the number of MultiFabs. It is hardcoded for 3 right now as it + * is only used for the electric field in the implicit electromagnetic time solvers. In + * the future, the array size can be made a template parameter so that this class can + * be used for other solver vectors, such as electrostatic (array size 1) or Darwin (array size 4). + */ + +class WarpXSolverVec +{ +public: + + WarpXSolverVec() = default; + + WarpXSolverVec(const WarpXSolverVec&) = delete; + + ~WarpXSolverVec() = default; + + using value_type = amrex::Real; + using RT = value_type; + + [[nodiscard]] inline bool IsDefined () const { return m_is_defined; } + + inline + void Define (const WarpXSolverVec& a_vec) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + a_vec.IsDefined(), + "WarpXSolverVec::Define(a_vec) called with undefined a_vec"); + Define( a_vec.getVec() ); + } + + inline + void Define ( const amrex::Vector, 3 > >& a_solver_vec ) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !IsDefined(), + "WarpXSolverVec::Define() called on undefined WarpXSolverVec"); + m_field_vec.resize(m_num_amr_levels); + const int lev = 0; + for (int n=0; n<3; n++) { + const amrex::MultiFab& mf_model = *a_solver_vec[lev][n]; + m_field_vec[lev][n] = std::make_unique( mf_model.boxArray(), mf_model.DistributionMap(), + mf_model.nComp(), amrex::IntVect::TheZeroVector() ); + } + m_is_defined = true; + } + + void SetDotMask( const amrex::Vector& a_Geom ); + [[nodiscard]] RT dotProduct( const WarpXSolverVec& a_X ) const; + + inline + void Copy ( const amrex::Vector, 3 > >& a_solver_vec ) + { + AMREX_ASSERT_WITH_MESSAGE( + IsDefined(), + "WarpXSolverVec::Copy() called on undefined WarpXSolverVec"); + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n = 0; n < 3; ++n) { + amrex::MultiFab::Copy(*m_field_vec[lev][n], *a_solver_vec[lev][n], 0, 0, m_ncomp, amrex::IntVect::TheZeroVector() ); + } + } + } + + inline + void Copy ( const WarpXSolverVec& a_vec ) + { + AMREX_ASSERT_WITH_MESSAGE( + a_vec.IsDefined(), + "WarpXSolverVec::Copy(a_vec) called with undefined a_vec"); + if (!IsDefined()) { Define(a_vec); } + const amrex::Vector, 3 > >& field_vec = a_vec.getVec(); + Copy(field_vec); + } + + // Prohibit Copy assignment operator + WarpXSolverVec& operator= ( const WarpXSolverVec& a_vec ) = delete; + + // Move assignment operator + WarpXSolverVec(WarpXSolverVec&&) noexcept = default; + WarpXSolverVec& operator= ( WarpXSolverVec&& a_vec ) noexcept + { + if (this != &a_vec) { + m_field_vec = std::move(a_vec.m_field_vec); + m_is_defined = true; + } + return *this; + } + + inline + void operator+= ( const WarpXSolverVec& a_vec ) + { + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n=0; n<3; n++) { + m_field_vec[lev][n]->plus(*(a_vec.getVec()[lev][n]), 0, 1, 0); + } + } + } + + inline + void operator-= (const WarpXSolverVec& a_vec) + { + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n=0; n<3; n++) { + m_field_vec[lev][n]->minus(*(a_vec.getVec()[lev][n]), 0, 1, 0); + } + } + } + + /** + * \brief Y = a*X + b*Y + */ + inline + void linComb (const RT a, const WarpXSolverVec& X, const RT b, const WarpXSolverVec& Y) + { + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n=0; n<3; n++) { + amrex::MultiFab::LinComb(*m_field_vec[lev][n], a, *X.getVec()[lev][n], 0, + b, *Y.getVec()[lev][n], 0, + 0, 1, 0); + } + } + } + + /** + * \brief Increment Y by a*X (Y += a*X) + */ + void increment (const WarpXSolverVec& X, const RT a) + { + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n=0; n<3; n++) { + amrex::MultiFab::Saxpy( *m_field_vec[lev][n], a, *X.getVec()[lev][n], + 0, 0, 1, amrex::IntVect::TheZeroVector() ); + } + } + } + + /** + * \brief Scale Y by a (Y *= a) + */ + inline + void scale (RT a_a) + { + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n=0; n<3; n++) { + m_field_vec[lev][n]->mult(a_a, 0, 1); + } + } + } + + inline + void zero () { setVal(0.0); } + + inline + void setVal ( const RT a_val ) + { + AMREX_ASSERT_WITH_MESSAGE( + IsDefined(), + "WarpXSolverVec::ones() called on undefined WarpXSolverVec"); + for (int lev = 0; lev < m_num_amr_levels; ++lev) { + for (int n=0; n<3; n++) { + m_field_vec[lev][n]->setVal(a_val); + } + } + } + + [[nodiscard]] inline RT norm2 () const + { + auto const norm = dotProduct(*this); + return std::sqrt(norm); + } + + [[nodiscard]] const amrex::Vector, 3 > >& getVec() const {return m_field_vec;} + amrex::Vector, 3 > >& getVec() {return m_field_vec;} + + // clearDotMask() must be called by the highest class that owns WarpXSolverVec() + // after it is done being used ( typically in the destructor ) to avoid the + // following error message after the simulation finishes: + // malloc_consolidate(): unaligned fastbin chunk detected + static void clearDotMask() { m_dotMask.clear(); } + +private: + + bool m_is_defined = false; + amrex::Vector, 3 > > m_field_vec; + + static constexpr int m_ncomp = 1; + static constexpr int m_num_amr_levels = 1; + + inline static bool m_dot_mask_defined = false; + inline static amrex::Vector,3>> m_dotMask; + +}; + +#endif diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp new file mode 100644 index 00000000000..b181f038fb5 --- /dev/null +++ b/Source/FieldSolver/ImplicitSolvers/WarpXSolverVec.cpp @@ -0,0 +1,51 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" + +void WarpXSolverVec::SetDotMask( const amrex::Vector& a_Geom ) +{ + if (m_dot_mask_defined) { return; } + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + IsDefined(), + "WarpXSolverVec::SetDotMask() called from undefined instance "); + + m_dotMask.resize(m_num_amr_levels); + for ( int n = 0; n < 3; n++) { + const amrex::BoxArray& grids = m_field_vec[0][n]->boxArray(); + const amrex::MultiFab tmp( grids, m_field_vec[0][n]->DistributionMap(), + 1, 0, amrex::MFInfo().SetAlloc(false) ); + const amrex::Periodicity& period = a_Geom[0].periodicity(); + m_dotMask[0][n] = tmp.OwnerMask(period); + } + m_dot_mask_defined = true; + + // If the function below is not called, then the following + // error message occurs after the simulation finishes: + // malloc_consolidate(): unaligned fastbin chunk detected + amrex::ExecOnFinalize(WarpXSolverVec::clearDotMask); +} + +[[nodiscard]] amrex::Real WarpXSolverVec::dotProduct ( const WarpXSolverVec& a_X ) const +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + m_dot_mask_defined, + "WarpXSolverVec::dotProduct called with m_dotMask not yet defined"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + a_X.IsDefined(), + "WarpXSolverVec::dotProduct(a_X) called with undefined a_X"); + amrex::Real result = 0.0; + const int lev = 0; + const bool local = true; + for (int n = 0; n < 3; ++n) { + auto rtmp = amrex::MultiFab::Dot( *m_dotMask[lev][n], + *m_field_vec[lev][n], 0, + *a_X.getVec()[lev][n], 0, 1, 0, local); + result += rtmp; + } + amrex::ParallelAllReduce::Sum(result, amrex::ParallelContext::CommunicatorSub()); + return result; +} diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package index 3b88f4627e5..a8af4c2de97 100644 --- a/Source/FieldSolver/Make.package +++ b/Source/FieldSolver/Make.package @@ -7,5 +7,6 @@ ifeq ($(USE_FFT),TRUE) endif include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/Make.package include $(WARPX_HOME)/Source/FieldSolver/MagnetostaticSolver/Make.package +include $(WARPX_HOME)/Source/FieldSolver/ImplicitSolvers/Make.package VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index bd40c5e0cd4..68b38153f8d 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -501,6 +501,9 @@ WarpX::InitData () CheckGuardCells(); PrintMainPICparameters(); + if (m_implicit_solver) { + m_implicit_solver->PrintParameters(); + } WriteUsedInputsFile(); if (restart_chkfile.empty()) @@ -568,10 +571,33 @@ WarpX::InitFromScratch () AmrCore::InitFromScratch(time); // This will call MakeNewLevelFromScratch + if (m_implicit_solver) { + + m_implicit_solver->Define(this); + m_implicit_solver->GetParticleSolverParams( max_particle_its_in_implicit_scheme, + particle_tol_in_implicit_scheme ); + + // Add space to save the positions and velocities at the start of the time steps + for (auto const& pc : *mypc) { +#if (AMREX_SPACEDIM >= 2) + pc->AddRealComp("x_n"); +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + pc->AddRealComp("y_n"); +#endif + pc->AddRealComp("z_n"); + pc->AddRealComp("ux_n"); + pc->AddRealComp("uy_n"); + pc->AddRealComp("uz_n"); + } + + } + mypc->AllocData(); mypc->InitData(); InitPML(); + } void diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 29f047c14e6..519eed52820 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -76,6 +76,7 @@ include $(WARPX_HOME)/Source/BoundaryConditions/Make.package include $(WARPX_HOME)/Source/Diagnostics/Make.package include $(WARPX_HOME)/Source/EmbeddedBoundary/Make.package include $(WARPX_HOME)/Source/FieldSolver/Make.package +include $(WARPX_HOME)/Source/NonlinearSolvers/Make.package include $(WARPX_HOME)/Source/Filter/Make.package include $(WARPX_HOME)/Source/Fluids/Make.package include $(WARPX_HOME)/Source/Initialization/Make.package @@ -90,7 +91,7 @@ include $(AMREX_HOME)/Src/Base/Make.package include $(AMREX_HOME)/Src/Particle/Make.package include $(AMREX_HOME)/Src/Boundary/Make.package include $(AMREX_HOME)/Src/AmrCore/Make.package -include $(AMREX_HOME)/Src/LinearSolvers/MLMG/Make.package +include $(AMREX_HOME)/Src/LinearSolvers/Make.package ifeq ($(USE_SENSEI_INSITU),TRUE) include $(AMREX_HOME)/Src/Amr/Make.package diff --git a/Source/NonlinearSolvers/CMakeLists.txt b/Source/NonlinearSolvers/CMakeLists.txt new file mode 100644 index 00000000000..938cdaad547 --- /dev/null +++ b/Source/NonlinearSolvers/CMakeLists.txt @@ -0,0 +1,6 @@ +foreach(D IN LISTS WarpX_DIMS) + warpx_set_suffix_dims(SD ${D}) + target_sources(lib_${SD} + PRIVATE + ) +endforeach() diff --git a/Source/NonlinearSolvers/JacobianFunctionMF.H b/Source/NonlinearSolvers/JacobianFunctionMF.H new file mode 100644 index 00000000000..823523df23c --- /dev/null +++ b/Source/NonlinearSolvers/JacobianFunctionMF.H @@ -0,0 +1,234 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef JacobianFunctionMF_H_ +#define JacobianFunctionMF_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. + * This class has all of the required functions to be used as the + * linear operator template parameter in AMReX_GMRES. + */ + +template +class JacobianFunctionMF +{ + public: + + using RT = typename T::value_type; + + JacobianFunctionMF() = default; + ~JacobianFunctionMF() = default; + + // Default move and copy operations + JacobianFunctionMF(const JacobianFunctionMF&) = default; + JacobianFunctionMF& operator=(const JacobianFunctionMF&) = default; + JacobianFunctionMF(JacobianFunctionMF&&) noexcept = default; + JacobianFunctionMF& operator=(JacobianFunctionMF&&) noexcept = default; + + void apply ( T& a_dF, const T& a_dU ); + + inline + void precond ( T& a_U, const T& a_X ) + { + if (m_usePreCond) { a_U.zero(); } + else { a_U.Copy(a_X); } + } + + inline + void updatePreCondMat ( const T& a_X ) + { + amrex::ignore_unused(a_X); + } + + inline + void create ( T& a_Z, const T& a_U ) + { + a_Z.define(a_U); + } + + T makeVecLHS () const; + T makeVecRHS () const; + + inline + void assign( T& a_Z, const T& a_U ) { + a_Z.Copy(a_U); + } + + inline + void increment( T& a_Z, const T& a_U, RT a_scale ) + { + a_Z.increment(a_U,a_scale); + } + + inline + void scale ( T& a_U, RT a_scale ) + { + a_U.scale(a_scale); + } + + inline + void linComb ( T& a_U, RT a, const T& X, RT b, const T& Y ) + { + a_U.linComb( a, X, b, Y ); + } + + inline + void setToZero ( T& a_U ) + { + a_U.zero(); + } + + inline + void setVal ( T& a_U, RT a_val ) + { + a_U.setVal(a_val); + } + + inline + RT dotProduct( const T& a_X, const T& a_Y ) + { + return( a_X.dotProduct(a_Y) ); + } + + inline + RT norm2( const T& a_U ) + { + return ( a_U.norm2() ); + } + + [[nodiscard]] inline + bool isDefined() const { return m_is_defined; } + + inline + void setBaseSolution ( const T& a_U ) + { + m_Y0.Copy(a_U); + m_normY0 = norm2(m_Y0); + } + + inline + void setBaseRHS ( const T& a_R ) + { + m_R0.Copy(a_R); + } + + inline + void setJFNKEps ( RT a_eps ) + { + m_epsJFNK = a_eps; + } + + inline + void setIsLinear ( bool a_isLinear ) + { + m_is_linear = a_isLinear; + } + + inline + void curTime ( RT a_time ) + { + m_cur_time = a_time; + } + + inline + void curTimeStep ( RT a_dt ) + { + m_dt = a_dt; + } + + void define( const T&, Ops* ); + + private: + + bool m_is_defined = false; + bool m_is_linear = false; + bool m_usePreCond = false; + 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; + +}; + +template +void JacobianFunctionMF::define ( const T& a_U, + Ops* a_ops ) +{ + m_Z.Define(a_U); + m_Y0.Define(a_U); + m_R0.Define(a_U); + m_R.Define(a_U); + + m_ops = a_ops; + + m_is_defined = true; +} + +template +auto JacobianFunctionMF::makeVecRHS () const -> T +{ + T Vec; + Vec.Define(m_R); + return Vec; +} + +template +auto JacobianFunctionMF::makeVecLHS () const -> T +{ + T Vec; + Vec.Define(m_R); + return Vec; +} + +template +void JacobianFunctionMF::apply (T& a_dF, const T& a_dU) +{ + BL_PROFILE("JacobianFunctionMF::apply()"); + using namespace amrex::literals; + RT const normY = norm2(a_dU); // always 1 when called from GMRES + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + isDefined(), + "JacobianFunction::apply() called on undefined JacobianFunction"); + + if (normY < 1.0e-15) { a_dF.zero(); } + else { + + RT eps; + if (m_is_linear) { + eps = 1.0_rt; + } else { + /* eps = error_rel * sqrt(1 + ||Y0||) / ||dU|| + * M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative Solver for + * Nonlinear Systems", SIAM J. Sci. Stat. Comput., 1998, vol 19, + * pp. 302--318. */ + if (m_normY0==0.0) { eps = m_epsJFNK * norm2(m_R0) / normY; } + else { + // m_epsJFNK * sqrt(1.0 + m_normY0) / normY + // above commonly used form not recommend for poorly scaled Y0 + eps = m_epsJFNK * m_normY0 / normY; + } + } + const RT eps_inv = 1.0_rt/eps; + + m_Z.linComb( 1.0, m_Y0, eps, a_dU ); // Z = Y0 + eps*dU + m_ops->ComputeRHS(m_R, m_Z, m_cur_time, m_dt, -1, true ); + + // F(Y) = Y - b - R(Y) ==> dF = dF/dY*dU = [1 - dR/dY]*dU + // = dU - (R(Z)-R(Y0))/eps + a_dF.linComb( 1.0, a_dU, eps_inv, m_R0 ); + a_dF.increment(m_R,-eps_inv); + + } + +} + +#endif diff --git a/Source/NonlinearSolvers/Make.package b/Source/NonlinearSolvers/Make.package new file mode 100644 index 00000000000..4f2dcfa9617 --- /dev/null +++ b/Source/NonlinearSolvers/Make.package @@ -0,0 +1 @@ +VPATH_LOCATIONS += $(WARPX_HOME)/Source/NonlinearSolvers diff --git a/Source/NonlinearSolvers/NewtonSolver.H b/Source/NonlinearSolvers/NewtonSolver.H new file mode 100644 index 00000000000..93ad432208a --- /dev/null +++ b/Source/NonlinearSolvers/NewtonSolver.H @@ -0,0 +1,345 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef NEWTON_SOLVER_H_ +#define NEWTON_SOLVER_H_ + +#include "NonlinearSolver.H" +#include "JacobianFunctionMF.H" + +#include +#include +#include "Utils/TextMsg.H" + +#include + +/** + * \brief Newton method to solve nonlinear equation of form: + * F(U) = U - b - R(U) = 0. U is the solution vector, b is a constant, + * and R(U) is some nonlinear function of U, which is computed in the + * ComputeRHS() Ops function. + */ + +template +class NewtonSolver : public NonlinearSolver +{ +public: + + NewtonSolver() = default; + + ~NewtonSolver() override = default; + + // Prohibit Move and Copy operations + NewtonSolver(const NewtonSolver&) = delete; + NewtonSolver& operator=(const NewtonSolver&) = delete; + NewtonSolver(NewtonSolver&&) noexcept = delete; + NewtonSolver& operator=(NewtonSolver&&) noexcept = delete; + + void Define ( const Vec& a_U, + Ops* a_ops ) override; + + void Solve ( Vec& a_U, + const Vec& a_b, + amrex::Real a_time, + amrex::Real a_dt ) const override; + + void GetSolverParams ( amrex::Real& a_rtol, + amrex::Real& a_atol, + int& a_maxits ) override + { + a_rtol = m_rtol; + a_atol = m_atol; + a_maxits = m_maxits; + } + + inline void CurTime ( amrex::Real a_time ) const + { + m_cur_time = a_time; + m_linear_function->curTime( a_time ); + } + + inline void CurTimeStep ( amrex::Real a_dt ) const + { + m_dt = a_dt; + m_linear_function->curTimeStep( a_dt ); + } + + void PrintParams () const override + { + amrex::Print() << "Newton verbose: " << (this->m_verbose?"true":"false") << std::endl; + amrex::Print() << "Newton max iterations: " << m_maxits << std::endl; + amrex::Print() << "Newton relative tolerance: " << m_rtol << std::endl; + amrex::Print() << "Newton absolute tolerance: " << m_atol << std::endl; + amrex::Print() << "Newton require convergence: " << (m_require_convergence?"true":"false") << std::endl; + amrex::Print() << "GMRES verbose: " << m_gmres_verbose_int << std::endl; + amrex::Print() << "GMRES restart length: " << m_gmres_restart_length << std::endl; + amrex::Print() << "GMRES max iterations: " << m_gmres_maxits << std::endl; + amrex::Print() << "GMRES relative tolerance: " << m_gmres_rtol << std::endl; + amrex::Print() << "GMRES absolute tolerance: " << m_gmres_atol << std::endl; + } + +private: + + /** + * \brief Intermediate Vec containers used by the solver. + */ + mutable Vec m_dU, m_F, m_R; + + /** + * \brief Pointer to Ops class. + */ + Ops* m_ops = nullptr; + + /** + * \brief Flag to determine whether convergence is required. + */ + bool m_require_convergence = true; + + /** + * \brief Relative tolerance for the Newton solver. + */ + amrex::Real m_rtol = 1.0e-6; + + /** + * \brief Absolute tolerance for the Newton solver. + */ + amrex::Real m_atol = 0.; + + /** + * \brief Maximum iterations for the Newton solver. + */ + int m_maxits = 100; + + /** + * \brief Relative tolerance for GMRES. + */ + amrex::Real m_gmres_rtol = 1.0e-4; + + /** + * \brief Absolute tolerance for GMRES. + */ + amrex::Real m_gmres_atol = 0.; + + /** + * \brief Maximum iterations for GMRES. + */ + int m_gmres_maxits = 1000; + + /** + * \brief Verbosity flag for GMRES. + */ + int m_gmres_verbose_int = 2; + + /** + * \brief Restart iteration for GMRES. + */ + int m_gmres_restart_length = 30; + + 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. + * In the contect of JFNK, A = dF/dU (i.e., system Jacobian) + */ + std::unique_ptr> m_linear_function; + + /** + * \brief The linear solver (GMRES) object. + */ + std::unique_ptr>> m_linear_solver; + + void ParseParameters (); + + /** + * \brief Compute the nonlinear residual: F(U) = U - b - R(U). + */ + void EvalResidual ( Vec& a_F, + const Vec& a_U, + const Vec& a_b, + amrex::Real a_time, + amrex::Real a_dt, + int a_iter ) const; + +}; + +template +void NewtonSolver::Define ( const Vec& a_U, + Ops* a_ops ) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !this->m_is_defined, + "Newton nonlinear solver object is already defined!"); + + ParseParameters(); + + m_dU.Define(a_U); + m_F.Define(a_U); // residual function F(U) = U - b - R(U) = 0 + m_R.Define(a_U); // right hand side function R(U) + + m_ops = a_ops; + + m_linear_function = std::make_unique>(); + m_linear_function->define(m_F, m_ops); + + m_linear_solver = std::make_unique>>(); + m_linear_solver->define(*m_linear_function); + m_linear_solver->setVerbose( m_gmres_verbose_int ); + m_linear_solver->setRestartLength( m_gmres_restart_length ); + m_linear_solver->setMaxIters( m_gmres_maxits ); + + this->m_is_defined = true; + +} + +template +void NewtonSolver::ParseParameters () +{ + const amrex::ParmParse pp_newton("newton"); + pp_newton.query("verbose", this->m_verbose); + pp_newton.query("absolute_tolerance", m_atol); + pp_newton.query("relative_tolerance", m_rtol); + pp_newton.query("max_iterations", m_maxits); + pp_newton.query("require_convergence", m_require_convergence); + + const amrex::ParmParse pp_gmres("gmres"); + pp_gmres.query("verbose_int", m_gmres_verbose_int); + pp_gmres.query("restart_length", m_gmres_restart_length); + pp_gmres.query("absolute_tolerance", m_gmres_atol); + pp_gmres.query("relative_tolerance", m_gmres_rtol); + pp_gmres.query("max_iterations", m_gmres_maxits); +} + +template +void NewtonSolver::Solve ( Vec& a_U, + const Vec& a_b, + amrex::Real a_time, + amrex::Real a_dt ) const +{ + BL_PROFILE("NewtonSolver::Solve()"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + this->m_is_defined, + "NewtonSolver::Solve() called on undefined object"); + using namespace amrex::literals; + + // + // Newton routine to solve nonlinear equation of form: + // F(U) = U - b - R(U) = 0 + // + + CurTime(a_time); + CurTimeStep(a_dt); + + amrex::Real norm_abs = 0.; + amrex::Real norm0 = 1._rt; + amrex::Real norm_rel = 0.; + + int iter; + for (iter = 0; iter < m_maxits;) { + + // Compute residual: F(U) = U - b - R(U) + EvalResidual(m_F, a_U, a_b, a_time, a_dt, iter); + + // Compute norm of the residual + norm_abs = m_F.norm2(); + if (iter == 0) { + if (norm_abs > 0.) { norm0 = norm_abs; } + else { norm0 = 1._rt; } + } + norm_rel = norm_abs/norm0; + + // Check for convergence criteria + if (this->m_verbose || iter == m_maxits) { + amrex::Print() << "Newton: iteration = " << std::setw(3) << iter << ", norm = " + << std::scientific << std::setprecision(5) << norm_abs << " (abs.), " + << std::scientific << std::setprecision(5) << norm_rel << " (rel.)" << "\n"; + } + + if (norm_abs < m_rtol) { + amrex::Print() << "Newton: exiting at iteration = " << std::setw(3) << iter + << ". Satisfied absolute tolerance " << m_atol << std::endl; + break; + } + + if (norm_rel < m_rtol) { + amrex::Print() << "Newton: exiting at iteration = " << std::setw(3) << iter + << ". Satisfied relative tolerance " << m_rtol << std::endl; + break; + } + + if (norm_abs > 100._rt*norm0) { + amrex::Print() << "Newton: exiting at iteration = " << std::setw(3) << iter + << ". SOLVER DIVERGED! relative tolerance = " << m_rtol << std::endl; + std::stringstream convergenceMsg; + convergenceMsg << "Newton: exiting at iteration " << std::setw(3) << iter << + ". SOLVER DIVERGED! absolute norm = " << norm_abs << + " has increased by 100X from that after first iteration."; + WARPX_ABORT_WITH_MESSAGE(convergenceMsg.str()); + } + + // Solve linear system for Newton step [Jac]*dU = F + m_dU.zero(); + m_linear_solver->solve( m_dU, m_F, m_gmres_rtol, m_gmres_atol ); + + // Update solution + a_U -= m_dU; + + iter++; + if (iter >= m_maxits) { + amrex::Print() << "Newton: exiting at iter = " << std::setw(3) << iter + << ". Maximum iteration reached: iter = " << m_maxits << std::endl; + break; + } + + } + + if (m_rtol > 0. && iter == m_maxits) { + std::stringstream convergenceMsg; + convergenceMsg << "Newton solver failed to converge after " << iter << + " iterations. Relative norm is " << norm_rel << + " and the relative tolerance is " << m_rtol << + ". Absolute norm is " << norm_abs << + " and the absolute tolerance is " << m_atol; + if (this->m_verbose) { amrex::Print() << convergenceMsg.str() << std::endl; } + if (m_require_convergence) { + WARPX_ABORT_WITH_MESSAGE(convergenceMsg.str()); + } else { + ablastr::warn_manager::WMRecordWarning("NewtonSolver", convergenceMsg.str()); + } + } + +} + +template +void NewtonSolver::EvalResidual ( Vec& a_F, + const Vec& a_U, + const Vec& a_b, + amrex::Real a_time, + amrex::Real a_dt, + int a_iter ) const +{ + + m_ops->ComputeRHS( m_R, a_U, a_time, a_dt, a_iter, false ); + + // set base U and R(U) for matrix-free Jacobian action calculation + m_linear_function->setBaseSolution(a_U); + 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; + + // Compute residual: F(U) = U - b - R(U) + a_F.Copy(a_U); + a_F -= m_R; + a_F -= a_b; + +} + +#endif diff --git a/Source/NonlinearSolvers/NonlinearSolver.H b/Source/NonlinearSolvers/NonlinearSolver.H new file mode 100644 index 00000000000..5587826474c --- /dev/null +++ b/Source/NonlinearSolvers/NonlinearSolver.H @@ -0,0 +1,87 @@ +#ifndef WARPX_NONLINEAR_SOLVER_H_ +#define WARPX_NONLINEAR_SOLVER_H_ + +#include +#include + +#include + +#include +#include +#include + +/** + * \brief Top-level class for the nonlinear solver + * + * This class is templated on a vector class Vec, and an operator class Ops. + * + * The Ops class must have the following function: + * ComputeRHS( R_vec, U_vec, time, dt, nl_iter, from_jacobian ), + * where U_vec and R_vec are of type Vec. + * + * The Vec class must have basic math operators, such as Copy, +=, -=, + * increment(), linComb(), scale(), etc.. See WarpXSolverVec.H for an example. + */ + +template +class NonlinearSolver +{ +public: + + NonlinearSolver() = default; + + virtual ~NonlinearSolver() = default; + + // Prohibit Move and Copy operations + NonlinearSolver(const NonlinearSolver&) = delete; + NonlinearSolver& operator=(const NonlinearSolver&) = delete; + NonlinearSolver(NonlinearSolver&&) noexcept = delete; + NonlinearSolver& operator=(NonlinearSolver&&) noexcept = delete; + + /** + * \brief Read user-provided parameters that control the nonlinear solver. + * Allocate intermediate data containers needed by the solver. For Newton, + * setup the linear solver for computing the Newton step. + */ + virtual void Define ( const Vec&, + Ops* ) = 0; + + /** + * \brief Check if the nonlinear solver has been defined. + */ + [[nodiscard]] bool IsDefined () const { return m_is_defined; } + + /** + * \brief Solve the specified nonlinear equation for U. + * Picard: U = b + R(U). + * Newton: F(U) = U - b - R(U) = 0. + */ + virtual void Solve ( Vec&, + const Vec&, + amrex::Real, + amrex::Real ) const = 0; + + /** + * \brief Print parameters used by the nonlinear solver. + */ + virtual void PrintParams () const = 0; + + /** + * \brief Return the convergence parameters used by the nonlinear solver. + */ + virtual void GetSolverParams (amrex::Real&, amrex::Real&, int&) = 0; + + /** + * \brief Allow caller to dynamically change the verbosity flag. For + * example, one may want to only print solver information every 100 steps. + */ + void Verbose ( bool a_verbose ) { m_verbose = a_verbose; } + +protected: + + bool m_is_defined = false; + mutable bool m_verbose = true; + +}; + +#endif diff --git a/Source/NonlinearSolvers/NonlinearSolverLibrary.H b/Source/NonlinearSolvers/NonlinearSolverLibrary.H new file mode 100644 index 00000000000..5803d6ca753 --- /dev/null +++ b/Source/NonlinearSolvers/NonlinearSolverLibrary.H @@ -0,0 +1,15 @@ +#ifndef NONLINEAR_SOLVER_LIBRARY_H_ +#define NONLINEAR_SOLVER_LIBRARY_H_ + +#include "PicardSolver.H" // IWYU pragma: export +#include "NewtonSolver.H" // IWYU pragma: export + +/** + * \brief struct to select the nonlinear solver for implicit schemes + */ +enum NonlinearSolverType { + Picard = 0, + Newton = 1 +}; + +#endif diff --git a/Source/NonlinearSolvers/PicardSolver.H b/Source/NonlinearSolvers/PicardSolver.H new file mode 100644 index 00000000000..f05b9a106e6 --- /dev/null +++ b/Source/NonlinearSolvers/PicardSolver.H @@ -0,0 +1,217 @@ +/* Copyright 2024 Justin Angus + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef PICARD_SOLVER_H_ +#define PICARD_SOLVER_H_ + +#include "NonlinearSolver.H" + +#include +#include "Utils/TextMsg.H" + +#include + +/** + * \brief Picard fixed-point iteration method to solve nonlinear + * equation of form: U = b + R(U). U is the solution vector. b + * is a constant. R(U) is some nonlinear function of U, which + * is computed in the Ops function ComputeRHS(). + */ + +template +class PicardSolver : public NonlinearSolver +{ +public: + + PicardSolver() = default; + + ~PicardSolver() override = default; + + // Prohibit Move and Copy operations + PicardSolver(const PicardSolver&) = delete; + PicardSolver& operator=(const PicardSolver&) = delete; + PicardSolver(PicardSolver&&) noexcept = delete; + PicardSolver& operator=(PicardSolver&&) noexcept = delete; + + void Define ( const Vec& a_U, + Ops* a_ops ) override; + + void Solve ( Vec& a_U, + const Vec& a_b, + amrex::Real a_time, + amrex::Real a_dt ) const override; + + void GetSolverParams ( amrex::Real& a_rtol, + amrex::Real& a_atol, + int& a_maxits ) override + { + a_rtol = m_rtol; + a_atol = m_atol; + a_maxits = m_maxits; + } + + void PrintParams () const override + { + amrex::Print() << "Picard max iterations: " << m_maxits << std::endl; + amrex::Print() << "Picard relative tolerance: " << m_rtol << std::endl; + amrex::Print() << "Picard absolute tolerance: " << m_atol << std::endl; + amrex::Print() << "Picard require convergence: " << (m_require_convergence?"true":"false") << std::endl; + } + +private: + + /** + * \brief Intermediate Vec containers used by the solver. + */ + mutable Vec m_Usave, m_R; + + /** + * \brief Pointer to Ops class. + */ + Ops* m_ops = nullptr; + + /** + * \brief Flag to determine whether convergence is required. + */ + bool m_require_convergence = true; + + /** + * \brief Relative tolerance for the Picard nonlinear solver + */ + amrex::Real m_rtol = 1.0e-6; + + /** + * \brief Absolute tolerance for the Picard nonlinear solver + */ + amrex::Real m_atol = 0.; + + /** + * \brief Maximum iterations for the Picard nonlinear solver + */ + int m_maxits = 100; + + void ParseParameters( ); + +}; + +template +void PicardSolver::Define ( const Vec& a_U, + Ops* a_ops ) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !this->m_is_defined, + "Picard nonlinear solver object is already defined!"); + + ParseParameters(); + + m_Usave.Define(a_U); + m_R.Define(a_U); + + m_ops = a_ops; + + this->m_is_defined = true; + +} + +template +void PicardSolver::ParseParameters () +{ + const amrex::ParmParse pp_picard("picard"); + pp_picard.query("verbose", this->m_verbose); + pp_picard.query("absolute_tolerance", m_atol); + pp_picard.query("relative_tolerance", m_rtol); + pp_picard.query("max_iterations", m_maxits); + pp_picard.query("require_convergence", m_require_convergence); + +} + +template +void PicardSolver::Solve ( Vec& a_U, + const Vec& a_b, + amrex::Real a_time, + amrex::Real a_dt ) const +{ + BL_PROFILE("PicardSolver::Solve()"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + this->m_is_defined, + "PicardSolver::Solve() called on undefined object"); + using namespace amrex::literals; + + // + // Picard fixed-point iteration method to solve nonlinear + // equation of form: U = b + R(U) + // + + amrex::Real norm_abs = 0.; + amrex::Real norm0 = 1._rt; + amrex::Real norm_rel = 0.; + + int iter; + for (iter = 0; iter < m_maxits;) { + + // Save previous state for norm calculation + m_Usave.Copy(a_U); + + // Update the solver state (a_U = a_b + m_R) + m_ops->ComputeRHS( m_R, a_U, a_time, a_dt, iter, false ); + a_U.Copy(a_b); + a_U += m_R; + + // Compute the step norm and update iter + m_Usave -= a_U; + norm_abs = m_Usave.norm2(); + if (iter == 0) { + if (norm_abs > 0.) { norm0 = norm_abs; } + else { norm0 = 1._rt; } + } + norm_rel = norm_abs/norm0; + iter++; + + // Check for convergence criteria + if (this->m_verbose || iter == m_maxits) { + amrex::Print() << "Picard: iter = " << std::setw(3) << iter << ", norm = " + << std::scientific << std::setprecision(5) << norm_abs << " (abs.), " + << std::scientific << std::setprecision(5) << norm_rel << " (rel.)" << "\n"; + } + + if (norm_abs < m_atol) { + amrex::Print() << "Picard: exiting at iter = " << std::setw(3) << iter + << ". Satisfied absolute tolerance " << m_atol << std::endl; + break; + } + + if (norm_rel < m_rtol) { + amrex::Print() << "Picard: exiting at iter = " << std::setw(3) << iter + << ". Satisfied relative tolerance " << m_rtol << std::endl; + break; + } + + if (iter >= m_maxits) { + amrex::Print() << "Picard: exiting at iter = " << std::setw(3) << iter + << ". Maximum iteration reached: iter = " << m_maxits << std::endl; + break; + } + + } + + if (m_rtol > 0. && iter == m_maxits) { + std::stringstream convergenceMsg; + convergenceMsg << "Picard solver failed to converge after " << iter << + " iterations. Relative norm is " << norm_rel << + " and the relative tolerance is " << m_rtol << + ". Absolute norm is " << norm_abs << + " and the absolute tolerance is " << m_atol; + if (this->m_verbose) { amrex::Print() << convergenceMsg.str() << std::endl; } + if (m_require_convergence) { + WARPX_ABORT_WITH_MESSAGE(convergenceMsg.str()); + } else { + ablastr::warn_manager::WMRecordWarning("PicardSolver", convergenceMsg.str()); + } + } + +} + +#endif diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H index fac94ff20a3..197cb897602 100644 --- a/Source/Particles/LaserParticleContainer.H +++ b/Source/Particles/LaserParticleContainer.H @@ -94,7 +94,7 @@ public: amrex::ParticleReal * AMREX_RESTRICT puzp, amrex::ParticleReal const * AMREX_RESTRICT pwp, amrex::Real const * AMREX_RESTRICT amplitude, - amrex::Real dt); + amrex::Real dt, PushType push_type=PushType::Explicit); protected: diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 7b735053f0b..0693a13c1f9 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -653,7 +653,7 @@ LaserParticleContainer::Evolve (int lev, // Calculate the corresponding momentum and position for the particles update_laser_particle(pti, static_cast(np), uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), wp.dataPtr(), - amplitude_E.dataPtr(), dt); + amplitude_E.dataPtr(), dt, push_type ); WARPX_PROFILE_VAR_STOP(blp_pp); // Current Deposition @@ -851,7 +851,7 @@ LaserParticleContainer::update_laser_particle (WarpXParIter& pti, ParticleReal * AMREX_RESTRICT const puzp, ParticleReal const * AMREX_RESTRICT const pwp, Real const * AMREX_RESTRICT const amplitude, - const Real dt) + const Real dt, PushType push_type) { const auto GetPosition = GetParticlePosition(pti); auto SetPosition = SetParticlePosition(pti); @@ -863,6 +863,26 @@ LaserParticleContainer::update_laser_particle (WarpXParIter& pti, const Real tmp_nvec_1 = m_nvec[1]; const Real tmp_nvec_2 = m_nvec[2]; + // When using the implicit solver, this function is called multiple times per timestep + // (within the linear and nonlinear solver). Thus, the position of the particles needs to be reset + // to the initial position (at the beginning of the timestep), before updating the particle position +#if (AMREX_SPACEDIM >= 2) + ParticleReal* x_n = nullptr; + if (push_type == PushType::Implicit) { + x_n = pti.GetAttribs(particle_comps["x_n"]).dataPtr(); + } +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + ParticleReal* y_n = nullptr; + if (push_type == PushType::Implicit) { + y_n = pti.GetAttribs(particle_comps["y_n"]).dataPtr(); + } +#endif + ParticleReal* z_n = nullptr; + if (push_type == PushType::Implicit) { + z_n = pti.GetAttribs(particle_comps["z_n"]).dataPtr(); + } + // Copy member variables to tmp copies for GPU runs. const Real tmp_mobility = m_mobility; const Real gamma_boost = WarpX::gamma_boost; @@ -894,15 +914,33 @@ LaserParticleContainer::update_laser_particle (WarpXParIter& pti, puzp[i] = gamma * vz; // Push the the particle positions - ParticleReal x, y, z; - GetPosition(i, x, y, z); + + // When using the implicit solver, this function is called multiple times per timestep + // (within the linear and nonlinear solver). Thus, the position of the particles needs to be reset + // to the initial position (at the beginning of the timestep), before updating the particle position + + ParticleReal x=0., y=0., z=0.; + if (push_type == PushType::Explicit) { + GetPosition(i, x, y, z); + } + #if !defined(WARPX_DIM_1D_Z) + if (push_type == PushType::Implicit) { + x = x_n[i]; + } x += vx * dt; #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + if (push_type == PushType::Implicit) { + y = y_n[i]; + } y += vy * dt; #endif + if (push_type == PushType::Implicit) { + z = z_n[i]; + } z += vz * dt; + SetPosition(i, x, y, z); } ); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3d134393f06..e2be4f948ca 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -61,6 +61,7 @@ #include #include #include +#include #include #include #include @@ -3118,6 +3119,12 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, const int qed_runtime_flag = no_qed; #endif + const int max_iterations = WarpX::max_particle_its_in_implicit_scheme; + const amrex::ParticleReal particle_tolerance = WarpX::particle_tol_in_implicit_scheme; + + amrex::Gpu::Buffer unconverged_particles({0}); + amrex::Long* unconverged_particles_ptr = unconverged_particles.data(); + // Using this version of ParallelFor with compile time options // improves performance when qed or external EB are not used by reducing // register pressure. @@ -3129,107 +3136,163 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, { // Position advance starts from the position at the start of the step // but uses the most recent velocity. + #if (AMREX_SPACEDIM >= 2) amrex::ParticleReal xp = x_n[ip]; const amrex::ParticleReal xp_n = x_n[ip]; #else - amrex::ParticleReal xp = 0._rt; + const amrex::ParticleReal xp = 0._rt; const amrex::ParticleReal xp_n = 0._rt; #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) amrex::ParticleReal yp = y_n[ip]; const amrex::ParticleReal yp_n = y_n[ip]; #else - amrex::ParticleReal yp = 0._rt; + const amrex::ParticleReal yp = 0._rt; const amrex::ParticleReal yp_n = 0._rt; #endif amrex::ParticleReal zp = z_n[ip]; const amrex::ParticleReal zp_n = z_n[ip]; - UpdatePositionImplicit(xp, yp, zp, ux_n[ip], uy_n[ip], uz_n[ip], ux[ip], uy[ip], uz[ip], 0.5_rt*dt); - setPosition(ip, xp, yp, zp); + amrex::ParticleReal dxp, dxp_save; + amrex::ParticleReal dyp, dyp_save; + amrex::ParticleReal dzp, dzp_save; + auto idxg2 = static_cast(1._rt/(dx[0]*dx[0])); + auto idyg2 = static_cast(1._rt/(dx[1]*dx[1])); + auto idzg2 = static_cast(1._rt/(dx[2]*dx[2])); - amrex::ParticleReal Exp = Ex_external_particle; - amrex::ParticleReal Eyp = Ey_external_particle; - amrex::ParticleReal Ezp = Ez_external_particle; - amrex::ParticleReal Bxp = Bx_external_particle; - amrex::ParticleReal Byp = By_external_particle; - amrex::ParticleReal Bzp = Bz_external_particle; + amrex::ParticleReal step_norm = 1._prt; + for (int iter=0; iter(ux[ip], uy[ip], uz[ip], - Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ion_lev ? ion_lev[ip] : 1, - m, q, pusher_algo, do_crr, -#ifdef WARPX_QED - t_chi_max, + if (!do_sync) #endif - dt); - } -#ifdef WARPX_QED - else { - if constexpr (qed_control == has_qed) { - doParticleMomentumPush<1>(ux[ip], uy[ip], uz[ip], + { + doParticleMomentumPush<0>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, ion_lev ? ion_lev[ip] : 1, m, q, pusher_algo, do_crr, +#ifdef WARPX_QED t_chi_max, +#endif dt); } - } +#ifdef WARPX_QED + else { + if constexpr (qed_control == has_qed) { + doParticleMomentumPush<1>(ux[ip], uy[ip], uz[ip], + Exp, Eyp, Ezp, Bxp, Byp, Bzp, + ion_lev ? ion_lev[ip] : 1, + m, q, pusher_algo, do_crr, + t_chi_max, + dt); + } + } #endif #ifdef WARPX_QED - [[maybe_unused]] auto foo_local_has_quantum_sync = local_has_quantum_sync; - [[maybe_unused]] auto *foo_podq = p_optical_depth_QSR; - [[maybe_unused]] const auto& foo_evolve_opt = evolve_opt; // have to do all these for nvcc - if constexpr (qed_control == has_qed) { - if (local_has_quantum_sync) { - evolve_opt(ux[ip], uy[ip], uz[ip], - Exp, Eyp, Ezp,Bxp, Byp, Bzp, - dt, p_optical_depth_QSR[ip]); + [[maybe_unused]] auto foo_local_has_quantum_sync = local_has_quantum_sync; + [[maybe_unused]] auto *foo_podq = p_optical_depth_QSR; + [[maybe_unused]] const auto& foo_evolve_opt = evolve_opt; // have to do all these for nvcc + if constexpr (qed_control == has_qed) { + if (local_has_quantum_sync) { + evolve_opt(ux[ip], uy[ip], uz[ip], + Exp, Eyp, Ezp,Bxp, Byp, Bzp, + dt, p_optical_depth_QSR[ip]); + } } - } #else amrex::ignore_unused(qed_control); #endif - // Take average to get the time centered value - ux[ip] = 0.5_rt*(ux[ip] + ux_n[ip]); - uy[ip] = 0.5_rt*(uy[ip] + uy_n[ip]); - uz[ip] = 0.5_rt*(uz[ip] + uz_n[ip]); + // Take average to get the time centered value + ux[ip] = 0.5_rt*(ux[ip] + ux_n[ip]); + uy[ip] = 0.5_rt*(uy[ip] + uy_n[ip]); + uz[ip] = 0.5_rt*(uz[ip] + uz_n[ip]); + + iter++; + + // particle did not converge + if ( iter > 1 && iter == max_iterations ) { +#if !defined(AMREX_USE_GPU) + std::stringstream convergenceMsg; + convergenceMsg << "Picard solver for particle failed to converge after " << + iter << " iterations. " << std::endl; + convergenceMsg << "Position step norm is " << step_norm << + " and the tolerance is " << particle_tolerance << std::endl; + convergenceMsg << " ux = " << ux[ip] << ", uy = " << uy[ip] << ", uz = " << uz[ip] << std::endl; + convergenceMsg << " xp = " << xp << ", yp = " << yp << ", zp = " << zp; + ablastr::warn_manager::WMRecordWarning("ImplicitPushXP", convergenceMsg.str()); +#endif + + // write signaling flag: how many particles did not converge? + amrex::Gpu::Atomic::Add(unconverged_particles_ptr, amrex::Long(1)); + } + + } // end Picard iterations }); + + auto const num_unconverged_particles = *(unconverged_particles.copyToHost()); + if (num_unconverged_particles > 0) { + ablastr::warn_manager::WMRecordWarning("ImplicitPushXP", + "Picard solver for " + + std::to_string(num_unconverged_particles) + + " particles failed to converge after " + + std::to_string(max_iterations) + " iterations." + ); + } } void diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 6e2e82632ee..89c2de88e47 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -85,4 +85,46 @@ void UpdatePositionImplicit(amrex::ParticleReal& x, amrex::ParticleReal& y, amre z += uz * inv_gamma * dt; } +/** \brief Check particle position for convergence. This is used by the theta-implicit + * and semi-implicit time solvers to obtain a self-consistent time-centered update + * of the particles for given electric and magnetic fields on the grid. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void PositionNorm( amrex::ParticleReal dxp, amrex::ParticleReal dyp, amrex::ParticleReal dzp, + amrex::ParticleReal& dxp_save, amrex::ParticleReal& dyp_save, amrex::ParticleReal& dzp_save, + amrex::ParticleReal idxg2, amrex::ParticleReal idyg2, amrex::ParticleReal idzg2, + amrex::ParticleReal& step_norm, const int iter ) +{ + using namespace amrex::literals; + +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(dxp, dxp_save, idxg2); +#endif +#if !defined(WARPX_DIM_3D) + amrex::ignore_unused(dyp, dyp_save, idyg2); +#endif + + if (iter==0) { step_norm = 1.0_prt; } + else { + step_norm = (dzp - dzp_save)*(dzp - dzp_save)*idzg2; +#if !defined(WARPX_DIM_1D_Z) + step_norm += (dxp - dxp_save)*(dxp - dxp_save)*idxg2; +#endif +#if defined(WARPX_DIM_3D) + step_norm += (dyp - dyp_save)*(dyp - dyp_save)*idyg2; +#elif defined(WARPX_DIM_RZ) + step_norm += (dyp - dyp_save)*(dyp - dyp_save)*idxg2; +#endif + step_norm = std::sqrt(step_norm); + } + dzp_save = dzp; +#if !defined(WARPX_DIM_1D_Z) + dxp_save = dxp; +#endif +#if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) + dyp_save = dyp; +#endif + +} + #endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index e067ed03a80..527b7766895 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -29,8 +29,8 @@ struct MediumForEM { struct EvolveScheme { enum { Explicit = 0, - ImplicitPicard = 1, - SemiImplicitPicard = 2 + ThetaImplicitEM = 1, + SemiImplicitEM = 2 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 2449ea2f8dd..89784c3b86c 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -24,10 +24,10 @@ // and corresponding integer for use inside the code const std::map evolve_scheme_to_int = { - {"explicit", EvolveScheme::Explicit }, - {"implicit_picard", EvolveScheme::ImplicitPicard }, - {"semi_implicit_picard", EvolveScheme::SemiImplicitPicard }, - {"default", EvolveScheme::Explicit } + {"explicit", EvolveScheme::Explicit }, + {"theta_implicit_em", EvolveScheme::ThetaImplicitEM }, + {"semi_implicit_em", EvolveScheme::SemiImplicitEM }, + {"default", EvolveScheme::Explicit } }; const std::map grid_to_int = { diff --git a/Source/WarpX.H b/Source/WarpX.H index df06a1723ca..a0a1379d9e2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -41,6 +41,8 @@ #include "FieldSolver/Fields.H" #include "FieldSolver/ElectrostaticSolver.H" #include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H" +#include "FieldSolver/ImplicitSolvers/ImplicitSolver.H" +#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H" #include "Filter/BilinearFilter.H" #include "Parallelization/GuardCellManager.H" #include "Utils/Parser/IntervalsParser.H" @@ -113,11 +115,27 @@ public: void Evolve (int numsteps = -1); - void EvolveImplicitPicardInit (int lev); - void SaveParticlesAtImplicitStepStart (WarpXParticleContainer& pc, int lev); - void FinishImplicitParticleUpdate (WarpXParticleContainer& pc, int lev); - void FinishImplicitFieldUpdate(amrex::Vector, 3 > >& Efield_fp, - amrex::Vector, 3 > >& Efield_n); + // + // Functions used by implicit solvers + // + void ImplicitPreRHSOp ( amrex::Real cur_time, + amrex::Real a_full_dt, + int a_nl_iter, + bool a_from_jacobian ); + void SaveParticlesAtImplicitStepStart (); + void FinishImplicitParticleUpdate (); + void SetElectricFieldAndApplyBCs ( const WarpXSolverVec& a_E ); + void UpdateMagneticFieldAndApplyBCs ( const amrex::Vector, 3 > >& a_Bn, + amrex::Real a_thetadt ); + void ApplyMagneticFieldBCs (); + void FinishMagneticFieldAndApplyBCs ( const amrex::Vector, 3 > >& a_Bn, + amrex::Real a_theta ); + void FinishImplicitField ( amrex::Vector, 3 > >& Field_fp, + const amrex::Vector, 3 > >& Field_n, + amrex::Real theta ); + void ImplicitComputeRHSE ( amrex::Real dt, WarpXSolverVec& a_Erhs_vec); + void ImplicitComputeRHSE (int lev, amrex::Real dt, WarpXSolverVec& a_Erhs_vec); + void ImplicitComputeRHSE (int lev, PatchType patch_type, amrex::Real dt, WarpXSolverVec& a_Erhs_vec); MultiParticleContainer& GetPartContainer () { return *mypc; } MultiFluidContainer& GetFluidContainer () { return *myfl; } @@ -159,14 +177,12 @@ public: static short particle_pusher_algo; //! Integer that corresponds to the type of Maxwell solver (Yee, CKC, PSATD, ECT) static short electromagnetic_solver_id; - //! Integer that corresponds to the evolve scheme (explicit, implicit_picard, semi_implicit_picard) + //! Integer that corresponds to the evolve scheme (explicit, semi_implicit_em, theta_implicit_em) static short evolve_scheme; - //! The maximum number of Picard iterations to do each time step - static int max_picard_iterations; - //! The tolerance for the Picard iteration convergence - static amrex::Real picard_iteration_tolerance; - //! Flags whether the Picard iterations are required to converge - static bool require_picard_convergence; + //! Maximum iterations used for self-consistent particle update in implicit particle-suppressed evolve schemes + static int max_particle_its_in_implicit_scheme; + //! Relative tolerance used for self-consistent particle update in implicit particle-suppressed evolve schemes + static amrex::ParticleReal particle_tol_in_implicit_scheme; /** Records a number corresponding to the load balance cost update strategy * being used (0 or 1 corresponding to timers or heuristic). */ @@ -506,6 +522,17 @@ public: [[nodiscard]] const amrex::MultiFab& getField(warpx::fields::FieldType field_type, int lev, int direction = 0) const; + /** + * \brief + * Get a constant reference to the specified vector field on the different MR levels + * + * \param field_type[in] the field type + * + * \return a vector (which one element per MR level) of arrays of three pointers (for 3 vector components) amrex::MultiFab* containing the field data + */ + [[nodiscard]] const amrex::Vector,3>>& + getMultiLevelField(warpx::fields::FieldType field_type) const; + [[nodiscard]] bool DoPML () const {return do_pml;} [[nodiscard]] bool DoFluidSpecies () const {return do_fluid_species;} @@ -1268,8 +1295,6 @@ private: void OneStep_nosub (amrex::Real cur_time); void OneStep_sub1 (amrex::Real cur_time); - void OneStep_ImplicitPicard(amrex::Real cur_time); - /** * \brief Perform one PIC iteration, with the multiple J deposition per time step */ @@ -1477,12 +1502,6 @@ private: amrex::Vector, 3 > > Efield_avg_fp; amrex::Vector, 3 > > Bfield_avg_fp; - // Implicit, fields at start of step and from the previous iteration - amrex::Vector, 3 > > Efield_n; - amrex::Vector, 3 > > Bfield_n; - amrex::Vector, 3 > > Efield_save; - amrex::Vector, 3 > > Bfield_save; - // Memory buffers for computing magnetostatic fields // Vector Potential A and previous step. Time buffer needed for computing dA/dt to first order amrex::Vector, 3 > > vector_potential_fp_nodal; @@ -1936,6 +1955,10 @@ private: amrex::Vector> m_fdtd_solver_fp; amrex::Vector> m_fdtd_solver_cp; + + // implicit solver object + std::unique_ptr m_implicit_solver; + }; #endif diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 512ed7fe9d0..d3f91002ef4 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -41,6 +41,8 @@ #include "Utils/WarpXProfilerWrapper.H" #include "Utils/WarpXUtil.H" +#include "FieldSolver/ImplicitSolvers/ImplicitSolverLibrary.H" + #include #include @@ -113,9 +115,8 @@ short WarpX::field_gathering_algo; short WarpX::particle_pusher_algo; short WarpX::electromagnetic_solver_id; short WarpX::evolve_scheme; -int WarpX::max_picard_iterations = 10; -Real WarpX::picard_iteration_tolerance = 1.e-7; -bool WarpX::require_picard_convergence = true; +int WarpX::max_particle_its_in_implicit_scheme = 21; +ParticleReal WarpX::particle_tol_in_implicit_scheme = 1.e-10; short WarpX::psatd_solution_type; short WarpX::J_in_time; short WarpX::rho_in_time; @@ -1195,6 +1196,21 @@ WarpX::ReadParameters () particle_pusher_algo = static_cast(GetAlgorithmInteger(pp_algo, "particle_pusher")); evolve_scheme = static_cast(GetAlgorithmInteger(pp_algo, "evolve_scheme")); + // check for implicit evolve scheme + if (evolve_scheme == EvolveScheme::SemiImplicitEM) { + m_implicit_solver = std::make_unique(); + } + else if (evolve_scheme == EvolveScheme::ThetaImplicitEM) { + m_implicit_solver = std::make_unique(); + } + + // implicit evolve schemes not setup to use mirrors + if (evolve_scheme == EvolveScheme::SemiImplicitEM || + evolve_scheme == EvolveScheme::ThetaImplicitEM) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( num_mirrors == 0, + "Mirrors cannot be used with Implicit evolve schemes."); + } + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo != CurrentDepositionAlgo::Esirkepov || !do_current_centering, @@ -1231,8 +1247,8 @@ WarpX::ReadParameters () if (current_deposition_algo == CurrentDepositionAlgo::Villasenor) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - evolve_scheme == EvolveScheme::ImplicitPicard || - evolve_scheme == EvolveScheme::SemiImplicitPicard, + evolve_scheme == EvolveScheme::SemiImplicitEM || + evolve_scheme == EvolveScheme::ThetaImplicitEM, "Villasenor current deposition can only" "be used with Implicit evolve schemes."); } @@ -1293,11 +1309,8 @@ WarpX::ReadParameters () macroscopic_solver_algo = GetAlgorithmInteger(pp_algo,"macroscopic_sigma_method"); } - if (evolve_scheme == EvolveScheme::ImplicitPicard || - evolve_scheme == EvolveScheme::SemiImplicitPicard) { - utils::parser::queryWithParser(pp_algo, "max_picard_iterations", max_picard_iterations); - utils::parser::queryWithParser(pp_algo, "picard_iteration_tolerance", picard_iteration_tolerance); - utils::parser::queryWithParser(pp_algo, "require_picard_convergence", require_picard_convergence); + if (evolve_scheme == EvolveScheme::SemiImplicitEM || + evolve_scheme == EvolveScheme::ThetaImplicitEM) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo == CurrentDepositionAlgo::Esirkepov || @@ -2177,11 +2190,6 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d AllocLevelMFs(lev, ba, dm, guard_cells.ng_alloc_EB, guard_cells.ng_alloc_J, guard_cells.ng_alloc_Rho, guard_cells.ng_alloc_F, guard_cells.ng_alloc_G, aux_is_nodal); - if (evolve_scheme == EvolveScheme::ImplicitPicard || - evolve_scheme == EvolveScheme::SemiImplicitPicard) { - EvolveImplicitPicardInit(lev); - } - m_accelerator_lattice[lev] = std::make_unique(); m_accelerator_lattice[lev]->InitElementFinder(lev, ba, dm); @@ -3513,3 +3521,32 @@ WarpX::getField(FieldType field_type, const int lev, const int direction) const { return *getFieldPointer(field_type, lev, direction); } + +const amrex::Vector,3>>& +WarpX::getMultiLevelField(warpx::fields::FieldType field_type) const +{ + switch(field_type) + { + case FieldType::Efield_aux : + return Efield_aux; + case FieldType::Bfield_aux : + return Bfield_aux; + case FieldType::Efield_fp : + return Efield_fp; + case FieldType::Bfield_fp : + return Bfield_fp; + case FieldType::current_fp : + return current_fp; + case FieldType::current_fp_nodal : + return current_fp_nodal; + case FieldType::Efield_cp : + return Efield_cp; + case FieldType::Bfield_cp : + return Bfield_cp; + case FieldType::current_cp : + return current_cp; + default: + WARPX_ABORT_WITH_MESSAGE("Invalid field type"); + return Efield_fp; + } +}