From 7020e9cee4f7b4894b8c489e43fcd47a2f91c794 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 11 Dec 2023 18:04:37 -0800 Subject: [PATCH 01/28] WIP implementation of the Darwin ES solver --- Source/Diagnostics/Diagnostics.cpp | 3 +- Source/Evolve/WarpXEvolve.cpp | 2 +- Source/FieldSolver/CMakeLists.txt | 1 + Source/FieldSolver/ElectrostaticSolver.cpp | 8 + .../ElectrostaticSolver/CMakeLists.txt | 7 + .../ImplicitDarwinSolver.H | 76 +++++ .../ImplicitDarwinSolver.cpp | 296 ++++++++++++++++++ .../ImplicitDarwinSolver_fwd.H | 15 + .../ElectrostaticSolver/Make.package | 3 + Source/FieldSolver/Make.package | 1 + Source/Initialization/WarpXInitData.cpp | 4 + Source/Utils/WarpXAlgorithmSelection.H | 3 +- Source/Utils/WarpXAlgorithmSelection.cpp | 1 + Source/WarpX.H | 4 + Source/WarpX.cpp | 27 +- 15 files changed, 446 insertions(+), 5 deletions(-) create mode 100644 Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt create mode 100644 Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H create mode 100644 Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp create mode 100644 Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H create mode 100644 Source/FieldSolver/ElectrostaticSolver/Make.package diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index f6e2da74127..80a61564148 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -76,7 +76,8 @@ Diagnostics::BaseReadParameters () if (utils::algorithms::is_in(m_varnames_fields, "phi")){ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrame || - warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic, + warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameDarwinImplicit, "plot phi only works if do_electrostatic = labframe or do_electrostatic = labframe-electromagnetostatic"); } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 32f7a493916..b557836d878 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -238,7 +238,7 @@ WarpX::Evolve (int numsteps) // loop (i.e. immediately after a `Redistribute` and before particle // positions are next pushed) so that the particles do not deposit out of bounds // and so that the fields are at the correct time in the output. - bool const reset_fields = true; + bool const reset_fields = (electrostatic_solver_id != ElectrostaticSolverAlgo::LabFrameDarwinImplicit); ComputeSpaceChargeField( reset_fields ); if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { // Call Magnetostatic Solver to solve for the vector potential A and compute the diff --git a/Source/FieldSolver/CMakeLists.txt b/Source/FieldSolver/CMakeLists.txt index 859117eb214..25928f1d86c 100644 --- a/Source/FieldSolver/CMakeLists.txt +++ b/Source/FieldSolver/CMakeLists.txt @@ -9,6 +9,7 @@ foreach(D IN LISTS WarpX_DIMS) ) endforeach() +add_subdirectory(ElectrostaticSolver) add_subdirectory(FiniteDifferenceSolver) add_subdirectory(MagnetostaticSolver) add_subdirectory(ImplicitSolvers) diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 189f2f2bb0a..35bf77238b1 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -7,6 +7,7 @@ #include "WarpX.H" #include "FieldSolver/ElectrostaticSolver.H" +#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H" #include "Fluids/MultiFluidContainer.H" #include "Fluids/WarpXFluidContainer.H" #include "Parallelization/GuardCellManager.H" @@ -76,6 +77,13 @@ WarpX::ComputeSpaceChargeField (bool const reset_fields) electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { AddSpaceChargeFieldLabFrame(); } + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { + m_implicit_darwin_solver->AddSpaceChargeField( + rho_fp, phi_fp, Efield_fp, + self_fields_required_precision, self_fields_absolute_tolerance, + self_fields_max_iters, self_fields_verbosity + ); + } else { // Loop over the species and add their space-charge contribution to E and B. // Note that the fields calculated here does not include the E field diff --git a/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt b/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt new file mode 100644 index 00000000000..c3596fe4fdd --- /dev/null +++ b/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt @@ -0,0 +1,7 @@ +foreach(D IN LISTS WarpX_DIMS) + warpx_set_suffix_dims(SD ${D}) + target_sources(lib_${SD} + PRIVATE + ImplicitDarwinSolver.cpp + ) +endforeach() diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H new file mode 100644 index 00000000000..7ce651d0820 --- /dev/null +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H @@ -0,0 +1,76 @@ +/* Copyright 2023 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Roelof Groenewald (TAE Technologies) + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_IMPLICITDARWINSOLVER_H_ +#define WARPX_IMPLICITDARWINSOLVER_H_ + +#include "ImplicitDarwinSolver_fwd.H" + +#include "Particles/MultiParticleContainer.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include +#include +#include + +#include +#include + +/** + * \brief This class contains the parameters needed to evaluate the implicit + * electrostatic Darwin model. + */ +class ImplicitDarwinSolver +{ +public: + ImplicitDarwinSolver (int nlevs_max); // constructor + + /** Allocate multifabs needed for the Darwin model. Called in constructor. */ + void AllocateMFs (int nlevs_max); + void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, + int ncomps, const amrex::IntVect& ngRho, const amrex::IntVect& rho_nodal_flag); + + /** Helper function to clear values from multifabs. */ + void ClearLevel (int lev); + + /** + * \brief + * Function to update the E-field using the implicit Darwin model. + */ + void AddSpaceChargeField ( + amrex::Vector>& rhofield, + amrex::Vector>& phifield, + amrex::Vector, 3>>& Efield, + amrex::Real required_precision, amrex::Real absolute_tolerance, + int max_iters, int verbosity + ); + + template + void ComputePhi ( + const amrex::Vector >& rho, + amrex::Vector >& phi, + amrex::Real relative_tolerance, amrex::Real absolute_tolerance, + int max_iters, int verbosity, + amrex::Geometry const geom, + T_BoundaryHandler const boundary_handler, + std::optional > eb_farray_box_factory = std::nullopt + ) const; + + void ComputeSigma (); + + // Declare multifabs specifically needed for the implicit Darwin model + // amrex::Vector> sigma_fp; + amrex::Vector, AMREX_SPACEDIM > > sigma; + // amrex::Vector > sigma; + +}; + +#endif // WARPX_IMPLICITDARWINSOLVER_H_ diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp new file mode 100644 index 00000000000..23f56094111 --- /dev/null +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp @@ -0,0 +1,296 @@ +/* Copyright 2023 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Roelof Groenewald (TAE Technologies) + * + * License: BSD-3-Clause-LBNL + */ + +#include "ImplicitDarwinSolver.H" + +using namespace amrex; + +ImplicitDarwinSolver::ImplicitDarwinSolver ( int nlevs_max ) +{ + AllocateMFs(nlevs_max); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (nlevs_max == 1), + "Implicit Darwin solver only support one level at present" + ); +} + +void ImplicitDarwinSolver::AllocateMFs (int nlevs_max) +{ + sigma.resize(nlevs_max); +} + +void ImplicitDarwinSolver::AllocateLevelMFs ( + int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, + int ncomps, const amrex::IntVect& ngRho, const amrex::IntVect& rho_nodal_flag) +{ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (rho_nodal_flag == IntVect::TheNodeVector()), + "Implicit Darwin solver only support a nodal charge density field at present" + ); + + // sigma[lev][0] = MultiFab( + // amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), + // dm, ncomps, ngRho + // ); + // sigma[lev][0].setVal(0.0_rt); + // sigma[lev][1] = MultiFab( + // amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), + // dm, ncomps, ngRho + // ); + // sigma[lev][1].setVal(0.0_rt); + + // The "sigma" multifabs stores the dressing of the Poisson equation that + // eliminates plasma frequency effects. + WarpX::AllocInitMultiFab( + sigma[lev][0], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), + dm, ncomps, ngRho, lev, "sigma[x]", 1.0_rt + ); + WarpX::AllocInitMultiFab( + sigma[lev][1], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), + dm, ncomps, ngRho, lev, "sigma[y]", 1.0_rt + ); + // WarpX::AllocInitMultiFab( + // sigma[lev][2], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), + // dm, ncomps, ngRho, lev, "sigma[z]", 0.0_rt + // ); + // sigma[lev][1] = std::make_unique( + // sigma[lev][0], amrex::make_alias, 0, ncomps + // ); + // sigma[lev][2] = std::make_unique( + // sigma[lev][0], amrex::make_alias, 0, ncomps + // ); + +#ifdef WARPX_DIM_RZ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (ncomps == 1), + "Implicit Darwin solver only support m = 0 azimuthal mode at present."); +#endif +} + +void ImplicitDarwinSolver::ClearLevel (int lev) +{ + // sigma[lev].reset(); + for (int i = 0; i < 3; ++i) { + sigma[lev][i].reset(); + } +} + +void +ImplicitDarwinSolver::AddSpaceChargeField ( + amrex::Vector>& rhofield, + amrex::Vector>& phifield, + amrex::Vector, 3>>& Efield, + Real const required_precision, Real absolute_tolerance, + int const max_iters, int const verbosity +) +{ + // Reset all E fields to 0, before calculating space-charge fields + WARPX_PROFILE("WarpX::ComputeSpaceChargeField::reset_fields"); + for (int lev = 0; lev <= 0; lev++) { + for (int comp=0; comp<3; comp++) { + Efield[lev][comp]->setVal(0); + // Bfield_fp[lev][comp]->setVal(0); + } + } + + WARPX_PROFILE("ImplicitDarwinSolver::AddSpaceChargeField"); + auto& warpx = WarpX::GetInstance(); + + // Store the boundary conditions for the field solver if they haven't been + // stored yet + if (!warpx.m_poisson_boundary_handler.bcs_set) + warpx.m_poisson_boundary_handler.definePhiBCs(warpx.Geom(0)); + + // Deposit particle charge density (source of Poisson solver) + auto& mypc = warpx.GetPartContainer(); + mypc.DepositCharge(rhofield, 0.0_rt); + // if (warpx.do_fluid_species) { + // int const lev = 0; + // myfl->DepositCharge( lev, *rho_fp[lev] ); + // } + + // Reflect density over PEC boundaries, if needed. + // filter (if used), exchange guard cells, interpolate across MR levels + // and apply boundary conditions + warpx.SyncCurrentAndRho(); + + // sigma dresses the Poisson equation to filter out plasma frequency effects + ComputeSigma(); + +#if defined(AMREX_USE_EB) + std::optional > eb_farray_box_factory; + amrex::Vector< + amrex::EBFArrayBoxFactory const * + > factories; + for (int lev = 0; lev <= finest_level; ++lev) { + factories.push_back(&warpx.fieldEBFactory(lev)); + } + eb_farray_box_factory = factories; +#else + const std::optional post_phi_calculation; + const std::optional > eb_farray_box_factory; +#endif + + // Use the AMREX MLMG solver + ComputePhi( + rhofield, phifield, required_precision, + absolute_tolerance, max_iters, verbosity, + warpx.Geom(0), + warpx.m_poisson_boundary_handler, + eb_farray_box_factory + ); + + // Compute the electric field. Note that if an EB is used the electric + // field will be calculated in the computePhi call above. +#ifndef AMREX_USE_EB + // beta is the moving frame velocity i.e. zero in lab frame + const std::array beta = {0._rt}; + + warpx.computeE( Efield, phifield, beta ); +#endif + +// // Compute the magnetic field +// computeB( Bfield_fp, phi_fp, beta ); +} + +void ImplicitDarwinSolver::ComputeSigma () { + + int const lev = 0; + + auto& warpx = WarpX::GetInstance(); + auto& mypc = warpx.GetPartContainer(); + auto& pc_electron = mypc.GetParticleContainerFromName("electrons"); + + // GetChargeDensity returns a cell-centered multifab + auto rho_electron = pc_electron.GetChargeDensity(lev, false); + + // multiply charge density by C_SI * dt**2 /4.0 q/(m \varepsilon_0) + rho_electron->mult( + (2.0_rt / 4.0_rt) * warpx.getdt(lev)*warpx.getdt(lev) + * pc_electron.getCharge() + / (pc_electron.getMass()) + ); + + // set sigma multifabs to 1 then add rho_electron to it + for (int ii = 0; ii < AMREX_SPACEDIM; ii++) + { + sigma[lev][ii]->setVal(1.0_rt); + MultiFab::Add( + *sigma[lev][ii], *rho_electron, 0, 0, 1, + sigma[lev][ii]->nGrowVect() + ); + } + + + // // Temporary cell-centered, single-component MultiFab for storing particles per cell. + // amrex::MultiFab ppc_mf(warpx.boxArray(lev), warpx.DistributionMap(lev), 1, 1); + // // Set value to 0, and increment the value in each cell with ppc. + // ppc_mf.setVal(0._rt); + // // Compute ppc which includes a summation over all species. + // pc_electron.Increment(ppc_mf, lev); +} + +template +void ImplicitDarwinSolver::ComputePhi ( + const amrex::Vector >& rho, + amrex::Vector >& phi, + Real const relative_tolerance, + Real absolute_tolerance, + int const max_iters, + int const verbosity, + amrex::Geometry const geom, + T_BoundaryHandler const boundary_handler, + std::optional > eb_farray_box_factory ) const +{ + using namespace amrex::literals; + + int const lev = 0; + + // scale rho appropriately; also determine if rho is zero everywhere + amrex::Real max_norm_b = 0.0; + rho[lev]->mult(-1._rt/ablastr::constant::SI::ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! + max_norm_b = amrex::max(max_norm_b, rho[lev]->norm0()); + + amrex::ParallelDescriptor::ReduceRealMax(max_norm_b); + + const bool always_use_bnorm = (max_norm_b > 0); + if (!always_use_bnorm) { + if (absolute_tolerance == 0.0) absolute_tolerance = amrex::Real(1e-6); + ablastr::warn_manager::WMRecordWarning( + "ElectrostaticSolver", + "Max norm of rho is 0", + ablastr::warn_manager::WarnPriority::low + ); + } + + amrex::LPInfo info; + +#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) + // In the presence of EB or RZ: the solver assumes that the beam is + // propagating along one of the axes of the grid, i.e. that only *one* + // of the components of `beta` is non-negligible. + amrex::MLABecLaplacian linop( {geom}, {rho.boxArray()}, {rho[lev].dm()}, info +#if defined(AMREX_USE_EB) + , {eb_farray_box_factory.value()[lev]} +#endif + ); + + // Note: this assumes that the beam is propagating along + // one of the axes of the grid, i.e. that only *one* of the + // components of `beta` is non-negligible. // we use this +#if defined(WARPX_DIM_RZ) + linop.setSigma({0._rt, 1._rt-beta_solver[1]*beta_solver[1]}); +#else + linop.setSigma({AMREX_D_DECL( + 1._rt-beta_solver[0]*beta_solver[0], + 1._rt-beta_solver[1]*beta_solver[1], + 1._rt-beta_solver[2]*beta_solver[2])}); +#endif + +#if defined(AMREX_USE_EB) + // if the EB potential only depends on time, the potential can be passed + // as a float instead of a callable + if (boundary_handler.phi_EB_only_t) { + linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value())); + } + else + linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); +#endif +#else + // In the absence of EB and RZ: use a more generic solver + // that can handle beams propagating in any direction + amrex::MLABecLaplacian linop( + {geom}, {rho[lev]->boxArray()}, {rho[lev]->DistributionMap()}, info + ); + linop.setScalars(0.0, 1.0); + linop.setBCoeffs(lev, amrex::GetArrOfConstPtrs(sigma[lev]) ); // for the non-axis-aligned solver +#endif + + // Solve the Poisson equation + linop.setDomainBC( boundary_handler.lobc, boundary_handler.hibc ); +#ifdef WARPX_DIM_RZ + linop.setRZ(true); +#endif + amrex::MLMG mlmg(linop); // actual solver defined here + mlmg.setVerbose(verbosity); + mlmg.setMaxIter(max_iters); + mlmg.setAlwaysUseBNorm(always_use_bnorm); + + // Solve Poisson equation at lev + mlmg.solve( + GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), + relative_tolerance, absolute_tolerance + ); + + // Run additional operations, such as calculation of the E field for embedded boundaries + // if constexpr (!std::is_same::value) + // if (post_phi_calculation.has_value()) + // post_phi_calculation.value()(mlmg, lev); + +} diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H new file mode 100644 index 00000000000..80b0ad22858 --- /dev/null +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H @@ -0,0 +1,15 @@ +/* Copyright 2023 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Roelof Groenewald (TAE Technologies) + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_IMPLICITDARWINSOLVER_FWD_H +#define WARPX_IMPLICITDARWINSOLVER_FWD_H + +class ImplicitDarwinSolver; + +#endif /* WARPX_IMPLICITDARWINSOLVER_FWD_H */ diff --git a/Source/FieldSolver/ElectrostaticSolver/Make.package b/Source/FieldSolver/ElectrostaticSolver/Make.package new file mode 100644 index 00000000000..c087c171c05 --- /dev/null +++ b/Source/FieldSolver/ElectrostaticSolver/Make.package @@ -0,0 +1,3 @@ +CEXE_sources += ImplicitDarwinSolver.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/ElectrostaticSolver diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package index a8af4c2de97..167b1f69bfd 100644 --- a/Source/FieldSolver/Make.package +++ b/Source/FieldSolver/Make.package @@ -5,6 +5,7 @@ CEXE_sources += WarpX_QED_Field_Pushers.cpp ifeq ($(USE_FFT),TRUE) include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/Make.package endif +include $(WARPX_HOME)/Source/FieldSolver/ElectrostaticSolver/Make.package include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/Make.package include $(WARPX_HOME)/Source/FieldSolver/MagnetostaticSolver/Make.package include $(WARPX_HOME)/Source/FieldSolver/ImplicitSolvers/Make.package diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index eb93d121313..bf636df01bf 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -267,6 +267,10 @@ WarpX::PrintMainPICparameters () amrex::Print() << "Operation mode: | Electrostatic" << "\n"; amrex::Print() << " | - laboratory frame, electrostatic + magnetostatic" << "\n"; } + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit){ + amrex::Print() << "Operation mode: | Electrostatic" << "\n"; + amrex::Print() << " | - laboratory frame, Darwin implicit scheme" << "\n"; + } else{ amrex::Print() << "Operation mode: | Electromagnetic" << "\n"; } diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index b9105557462..06429e64ca5 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -67,7 +67,8 @@ struct ElectrostaticSolverAlgo { None = 0, Relativistic = 1, LabFrameElectroMagnetostatic = 2, - LabFrame = 3 // Non relativistic + LabFrame = 3, // Non relativistic + LabFrameDarwinImplicit = 4 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 9e2360315b8..c392defd35c 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -53,6 +53,7 @@ const std::map electrostatic_solver_algo_to_int = { {"relativistic", ElectrostaticSolverAlgo::Relativistic}, {"labframe-electromagnetostatic", ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic}, {"labframe", ElectrostaticSolverAlgo::LabFrame}, + {"labframe-implicit", ElectrostaticSolverAlgo::LabFrameDarwinImplicit}, {"default", ElectrostaticSolverAlgo::None } }; diff --git a/Source/WarpX.H b/Source/WarpX.H index 2caea95e926..0c1cefbb9e6 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -16,6 +16,7 @@ #include "Diagnostics/MultiDiagnostics_fwd.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags_fwd.H" #include "EmbeddedBoundary/WarpXFaceInfoBox_fwd.H" +#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel_fwd.H" @@ -1635,6 +1636,9 @@ private: // Hybrid PIC algorithm parameters std::unique_ptr m_hybrid_pic_model; + // Implicit Darwin electrostatic algorithm + std::unique_ptr m_implicit_darwin_solver; + // Load balancing /** Load balancing intervals that reads the "load_balance_intervals" string int the input file * for getting steps at which load balancing is performed */ diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 9ed59c348b3..cd0562073b9 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -15,6 +15,7 @@ #include "Diagnostics/MultiDiagnostics.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "EmbeddedBoundary/WarpXFaceInfoBox.H" +#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H" @@ -333,6 +334,12 @@ WarpX::WarpX () vector_potential_grad_buf_b_stag.resize(nlevs_max); } + // Initialize the implicit Darwin electrostatic solver if required + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + { + m_implicit_darwin_solver = std::make_unique(nlevs_max); + } + if (fft_do_time_averaging) { Efield_avg_fp.resize(nlevs_max); @@ -743,7 +750,8 @@ WarpX::ReadParameters () } if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { // Note that with the relativistic version, these parameters would be // input for each species. @@ -2097,6 +2105,11 @@ WarpX::ClearLevel (int lev) current_buf[lev][i].reset(); } + if (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + { + m_implicit_darwin_solver->ClearLevel(lev); + } + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { m_hybrid_pic_model->ClearLevel(lev); @@ -2357,6 +2370,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm dm, ncomps, ngEB, lev, "vector_potential_grad_buf_b_stag[z]", 0.0_rt); } + // Allocate extra multifabs needed by the implicit Darwin electrostatic algorithm. + if (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + { + m_implicit_darwin_solver->AllocateLevelMFs( + lev, ba, dm, ncomps, ngRho, rho_nodal_flag + ); + } + // Allocate extra multifabs needed by the kinetic-fluid hybrid algorithm. if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { @@ -2441,6 +2462,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm int rho_ncomps = 0; if( (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) || (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) || + (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) || (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) ) { rho_ncomps = ncomps; } @@ -2459,7 +2481,8 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm } if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit ) { const IntVect ngPhi = IntVect( AMREX_D_DECL(1,1,1) ); AllocInitMultiFab(phi_fp[lev], amrex::convert(ba, phi_nodal_flag), dm, ncomps, ngPhi, lev, "phi_fp", 0.0_rt); From ac4bb0f3055e6862a717111cacb7ff1a75def243 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Fri, 22 Mar 2024 09:47:01 -0700 Subject: [PATCH 02/28] add picmi flag to use Darwin solver --- Python/pywarpx/picmi.py | 3 +++ Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H | 2 +- .../FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp | 2 +- Source/Python/WarpX.cpp | 1 + 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index a7ef3470c52..5b952eee9b6 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1551,6 +1551,7 @@ def init(self, kw): self.absolute_tolerance = kw.pop('warpx_absolute_tolerance', None) self.self_fields_verbosity = kw.pop('warpx_self_fields_verbosity', None) self.magnetostatic = kw.pop('warpx_magnetostatic', False) + self.darwin_labframe = kw.pop('warpx_darwin_labframe', False) def solver_initialize_inputs(self): @@ -1564,6 +1565,8 @@ def solver_initialize_inputs(self): else: if self.magnetostatic: pywarpx.warpx.do_electrostatic = 'labframe-electromagnetostatic' + elif self.darwin_labframe: + pywarpx.warpx.do_electrostatic = 'labframe-implicit' else: pywarpx.warpx.do_electrostatic = 'labframe' pywarpx.warpx.self_fields_required_precision = self.required_precision diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H index 7ce651d0820..733cd0735a6 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H @@ -1,4 +1,4 @@ -/* Copyright 2023 The WarpX Community +/* Copyright 2024 The WarpX Community * * This file is part of WarpX. * diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp index 23f56094111..51dff4bb1bd 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp @@ -1,4 +1,4 @@ -/* Copyright 2023 The WarpX Community +/* Copyright 2024 The WarpX Community * * This file is part of WarpX. * diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp index b3bd5414303..57bde115af6 100644 --- a/Source/Python/WarpX.cpp +++ b/Source/Python/WarpX.cpp @@ -11,6 +11,7 @@ #include #include #include +#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H" #include #include #include From df0dc3eb4e3aef094de09a7e8bafd22478d77d56 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Fri, 12 Apr 2024 17:40:51 -0700 Subject: [PATCH 03/28] slowly progressing --- .../ImplicitDarwinSolver.H | 7 +- .../ImplicitDarwinSolver.cpp | 83 ++++++++++--------- 2 files changed, 48 insertions(+), 42 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H index 733cd0735a6..b501a667df0 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H @@ -21,7 +21,8 @@ #include #include -#include +#include +#include #include /** @@ -68,9 +69,9 @@ public: // Declare multifabs specifically needed for the implicit Darwin model // amrex::Vector> sigma_fp; - amrex::Vector, AMREX_SPACEDIM > > sigma; + // amrex::Vector, AMREX_SPACEDIM > > sigma; // amrex::Vector > sigma; - + amrex::Vector > sigma; }; #endif // WARPX_IMPLICITDARWINSOLVER_H_ diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp index 51dff4bb1bd..8cf9c90ef86 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp @@ -45,26 +45,12 @@ void ImplicitDarwinSolver::AllocateLevelMFs ( // ); // sigma[lev][1].setVal(0.0_rt); - // The "sigma" multifabs stores the dressing of the Poisson equation that - // eliminates plasma frequency effects. + // The "sigma" multifabs stores the dressing of the Poisson equation. It + // is a cell-centered multifab. WarpX::AllocInitMultiFab( - sigma[lev][0], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), - dm, ncomps, ngRho, lev, "sigma[x]", 1.0_rt + sigma[lev], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), + dm, ncomps, ngRho, lev, "sigma" ); - WarpX::AllocInitMultiFab( - sigma[lev][1], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), - dm, ncomps, ngRho, lev, "sigma[y]", 1.0_rt - ); - // WarpX::AllocInitMultiFab( - // sigma[lev][2], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), - // dm, ncomps, ngRho, lev, "sigma[z]", 0.0_rt - // ); - // sigma[lev][1] = std::make_unique( - // sigma[lev][0], amrex::make_alias, 0, ncomps - // ); - // sigma[lev][2] = std::make_unique( - // sigma[lev][0], amrex::make_alias, 0, ncomps - // ); #ifdef WARPX_DIM_RZ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -75,10 +61,10 @@ void ImplicitDarwinSolver::AllocateLevelMFs ( void ImplicitDarwinSolver::ClearLevel (int lev) { - // sigma[lev].reset(); - for (int i = 0; i < 3; ++i) { - sigma[lev][i].reset(); - } + sigma[lev].reset(); + // for (int i = 0; i < 3; ++i) { + // sigma[lev][i].reset(); + // } } void @@ -165,7 +151,7 @@ void ImplicitDarwinSolver::ComputeSigma () { auto& warpx = WarpX::GetInstance(); auto& mypc = warpx.GetPartContainer(); - auto& pc_electron = mypc.GetParticleContainerFromName("electrons"); + auto& pc_electron = mypc.GetParticleContainerFromName("electron"); // GetChargeDensity returns a cell-centered multifab auto rho_electron = pc_electron.GetChargeDensity(lev, false); @@ -177,17 +163,24 @@ void ImplicitDarwinSolver::ComputeSigma () { / (pc_electron.getMass()) ); - // set sigma multifabs to 1 then add rho_electron to it - for (int ii = 0; ii < AMREX_SPACEDIM; ii++) - { - sigma[lev][ii]->setVal(1.0_rt); - MultiFab::Add( - *sigma[lev][ii], *rho_electron, 0, 0, 1, - sigma[lev][ii]->nGrowVect() - ); - } + // interpolate rho_electron to cell-centered multifab + // set sigma multifabs to 1 then add rho_electron to it + // for (int ii = 0; ii < AMREX_SPACEDIM; ii++) + // { + // sigma[lev][ii]->setVal(1.0_rt); + // MultiFab::Add( + // *sigma[lev][ii], *rho_electron, 0, 0, 1, + // sigma[lev][ii]->nGrowVect() + // ); + // } + sigma[lev]->setVal(1.0_rt); + // MultiFab::Add( + // *sigma[lev][ii], *rho_electron, 0, 0, 1, + // sigma[lev][ii]->nGrowVect() + // ); + // // Temporary cell-centered, single-component MultiFab for storing particles per cell. // amrex::MultiFab ppc_mf(warpx.boxArray(lev), warpx.DistributionMap(lev), 1, 1); // // Set value to 0, and increment the value in each cell with ppc. @@ -223,19 +216,21 @@ void ImplicitDarwinSolver::ComputePhi ( if (!always_use_bnorm) { if (absolute_tolerance == 0.0) absolute_tolerance = amrex::Real(1e-6); ablastr::warn_manager::WMRecordWarning( - "ElectrostaticSolver", + "ImplicitDarwinSolver", "Max norm of rho is 0", ablastr::warn_manager::WarnPriority::low ); } + // auto& warpx = WarpX::GetInstance(); + // amrex::LPInfo info = LPInfo().setMaxCoarseningLevel(0); amrex::LPInfo info; #if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) // In the presence of EB or RZ: the solver assumes that the beam is // propagating along one of the axes of the grid, i.e. that only *one* // of the components of `beta` is non-negligible. - amrex::MLABecLaplacian linop( {geom}, {rho.boxArray()}, {rho[lev].dm()}, info + amrex::MLEBNodeLaplacian linop( {geom}, {rho.boxArray()}, {rho[lev].dm()}, info #if defined(AMREX_USE_EB) , {eb_farray_box_factory.value()[lev]} #endif @@ -251,6 +246,8 @@ void ImplicitDarwinSolver::ComputePhi ( 1._rt-beta_solver[0]*beta_solver[0], 1._rt-beta_solver[1]*beta_solver[1], 1._rt-beta_solver[2]*beta_solver[2])}); + // linop.setScalars(0.0, 1.0); + // linop.setBCoeffs(lev, amrex::GetArrOfConstPtrs(sigma[lev]) ); // for the non-axis-aligned solver #endif #if defined(AMREX_USE_EB) @@ -264,12 +261,12 @@ void ImplicitDarwinSolver::ComputePhi ( #endif #else // In the absence of EB and RZ: use a more generic solver - // that can handle beams propagating in any direction - amrex::MLABecLaplacian linop( + // that can handle sigma in all directions + amrex::MLNodeLaplacian linop( {geom}, {rho[lev]->boxArray()}, {rho[lev]->DistributionMap()}, info + // {}, 1.0 ); - linop.setScalars(0.0, 1.0); - linop.setBCoeffs(lev, amrex::GetArrOfConstPtrs(sigma[lev]) ); // for the non-axis-aligned solver + linop.setSigma(lev, *sigma[lev]); #endif // Solve the Poisson equation @@ -282,9 +279,17 @@ void ImplicitDarwinSolver::ComputePhi ( mlmg.setMaxIter(max_iters); mlmg.setAlwaysUseBNorm(always_use_bnorm); + // // create a vector to our fields, sorted by level + // amrex::Vector sorted_rho; + // amrex::Vector sorted_phi; + // // for (int nlev = 0; nlev <= finest_level; ++nlev) { + // sorted_rho.emplace_back(rho[lev].get()); + // sorted_phi.emplace_back(phi[lev].get()); + // // } + // Solve Poisson equation at lev mlmg.solve( - GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), + { phi[lev].get() }, { rho[lev].get() }, // GetVecOfPtrs(phi), GetVecOfConstPtrs(rho) // relative_tolerance, absolute_tolerance ); From 115bc692b941019b67593a0cca00318e32e4147e Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 15 Apr 2024 09:06:13 -0700 Subject: [PATCH 04/28] make C_SI a runtime parameter --- Python/pywarpx/picmi.py | 2 + .../ImplicitDarwinSolver.H | 4 + .../ImplicitDarwinSolver.cpp | 81 ++++++++++++------- 3 files changed, 57 insertions(+), 30 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 5b952eee9b6..2c6fe8b3bb0 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1552,6 +1552,7 @@ def init(self, kw): self.self_fields_verbosity = kw.pop('warpx_self_fields_verbosity', None) self.magnetostatic = kw.pop('warpx_magnetostatic', False) self.darwin_labframe = kw.pop('warpx_darwin_labframe', False) + self.semi_implicit_factor = kw.pop('warpx_semi_implicit_factor', None) def solver_initialize_inputs(self): @@ -1567,6 +1568,7 @@ def solver_initialize_inputs(self): pywarpx.warpx.do_electrostatic = 'labframe-electromagnetostatic' elif self.darwin_labframe: pywarpx.warpx.do_electrostatic = 'labframe-implicit' + pywarpx.warpx.semi_implicit_factor = self.semi_implicit_factor else: pywarpx.warpx.do_electrostatic = 'labframe' pywarpx.warpx.self_fields_required_precision = self.required_precision diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H index b501a667df0..bc4e3b35ccc 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H @@ -13,10 +13,12 @@ #include "ImplicitDarwinSolver_fwd.H" #include "Particles/MultiParticleContainer.H" +#include "Utils/Parser/ParserUtils.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" +#include #include #include #include @@ -72,6 +74,8 @@ public: // amrex::Vector, AMREX_SPACEDIM > > sigma; // amrex::Vector > sigma; amrex::Vector > sigma; + + amrex::Real m_C_SI = 6.0; }; #endif // WARPX_IMPLICITDARWINSOLVER_H_ diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp index 8cf9c90ef86..04f7edb73d7 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp @@ -18,6 +18,10 @@ ImplicitDarwinSolver::ImplicitDarwinSolver ( int nlevs_max ) (nlevs_max == 1), "Implicit Darwin solver only support one level at present" ); + + // read input parameters + const ParmParse pp_warpx("warpx"); + utils::parser::queryWithParser(pp_warpx, "semi_implicit_factor", m_C_SI); } void ImplicitDarwinSolver::AllocateMFs (int nlevs_max) @@ -148,45 +152,62 @@ ImplicitDarwinSolver::AddSpaceChargeField ( void ImplicitDarwinSolver::ComputeSigma () { int const lev = 0; + sigma[lev]->setVal(1.0_rt); + + // GetChargeDensity returns a nodal multifab + amrex::GpuArray nodal = {1, 1, 1}; + // sigma is a cell-centered array + amrex::GpuArray const& cell_centered = {0, 0, 0}; + // The "coarsening is just 1 i.e. no coarsening" + amrex::GpuArray const& coarsen = {1, 1, 1}; + + // Below we set all the unused dimensions to have cell-centered values for + // rho since these values will be interpolated onto a cell-centered grid + // - if this is not done the Interp function returns nonsense values. +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_1D_Z) + nodal[2] = 0; +#endif +#if defined(WARPX_DIM_1D_Z) + nodal[1] = 0; +#endif auto& warpx = WarpX::GetInstance(); auto& mypc = warpx.GetPartContainer(); - auto& pc_electron = mypc.GetParticleContainerFromName("electron"); - - // GetChargeDensity returns a cell-centered multifab - auto rho_electron = pc_electron.GetChargeDensity(lev, false); - // multiply charge density by C_SI * dt**2 /4.0 q/(m \varepsilon_0) - rho_electron->mult( - (2.0_rt / 4.0_rt) * warpx.getdt(lev)*warpx.getdt(lev) - * pc_electron.getCharge() - / (pc_electron.getMass()) + auto m_mult_factor = ( + m_C_SI * warpx.getdt(lev) * warpx.getdt(lev) + / (16._rt * MathConst::pi * MathConst::pi * PhysConst::ep0) ); - // interpolate rho_electron to cell-centered multifab + // Loop over each species to calculate the Poisson equation dressing + for (auto const& pc : mypc) { + auto rho = pc->GetChargeDensity(lev, false); + // Handle the parallel transfers of guard cells and + // apply the filtering if requested - might not be needed... + warpx.ApplyFilterandSumBoundaryRho(lev, lev, *rho, 0, rho->nComp()); + // multiply charge density by C_SI * dt**2 /4.0 q/(m \varepsilon_0) + rho->mult(m_mult_factor * pc->getCharge() / pc->getMass()); - // set sigma multifabs to 1 then add rho_electron to it - // for (int ii = 0; ii < AMREX_SPACEDIM; ii++) - // { - // sigma[lev][ii]->setVal(1.0_rt); - // MultiFab::Add( - // *sigma[lev][ii], *rho_electron, 0, 0, 1, - // sigma[lev][ii]->nGrowVect() - // ); - // } - sigma[lev]->setVal(1.0_rt); - // MultiFab::Add( - // *sigma[lev][ii], *rho_electron, 0, 0, 1, - // sigma[lev][ii]->nGrowVect() - // ); + // update sigma +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*sigma[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Array4 const& sigma_arr = sigma[lev]->array(mfi); + Array4 const& rho_arr = rho->const_array(mfi); + + // Loop over the cells and update the sigma field + amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int k){ + // Interpolate rho to cell-centered multifab to add to sigma. + auto const rho_interp = ablastr::coarsen::sample::Interp( + rho_arr, nodal, cell_centered, coarsen, i, j, k, 0 + ); + sigma_arr(i, j, k, 0) += rho_interp; + }); - // // Temporary cell-centered, single-component MultiFab for storing particles per cell. - // amrex::MultiFab ppc_mf(warpx.boxArray(lev), warpx.DistributionMap(lev), 1, 1); - // // Set value to 0, and increment the value in each cell with ppc. - // ppc_mf.setVal(0._rt); - // // Compute ppc which includes a summation over all species. - // pc_electron.Increment(ppc_mf, lev); + } + } } template From 8549b505661b4142af7d2eea97995efbe0a69edf Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Sat, 8 Jun 2024 23:04:42 -0700 Subject: [PATCH 05/28] use new sigma support in MLEBNodeLaplacian --- Source/FieldSolver/ElectrostaticSolver.cpp | 77 +++-- .../ImplicitDarwinSolver.H | 30 -- .../ImplicitDarwinSolver.cpp | 205 +----------- .../fields/SemiImplicitPoissonSolver.H | 303 ++++++++++++++++++ 4 files changed, 352 insertions(+), 263 deletions(-) create mode 100644 Source/ablastr/fields/SemiImplicitPoissonSolver.H diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 35bf77238b1..a5c803f45de 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -21,6 +21,7 @@ #include "Utils/WarpXProfilerWrapper.H" #include +#include #include #include @@ -74,16 +75,10 @@ WarpX::ComputeSpaceChargeField (bool const reset_fields) } if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit ) { AddSpaceChargeFieldLabFrame(); } - else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { - m_implicit_darwin_solver->AddSpaceChargeField( - rho_fp, phi_fp, Efield_fp, - self_fields_required_precision, self_fields_absolute_tolerance, - self_fields_max_iters, self_fields_verbosity - ); - } else { // Loop over the species and add their space-charge contribution to E and B. // Note that the fields calculated here does not include the E field @@ -284,6 +279,9 @@ WarpX::AddSpaceChargeFieldLabFrame () } else { #if defined(WARPX_DIM_1D_Z) + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { + amrex::Abort("Cannot use SI solver in 1d."); + } // Use the tridiag solver with 1D computePhiTriDiagonal(rho_fp, phi_fp); #else @@ -393,27 +391,48 @@ WarpX::computePhi (const amrex::Vector >& rho, bool const is_solver_igf_on_lev0 = WarpX::poisson_solver_id == PoissonSolverAlgo::IntegratedGreenFunction; - ablastr::fields::computePhi( - sorted_rho, - sorted_phi, - beta, - required_precision, - absolute_tolerance, - max_iters, - verbosity, - this->geom, - this->dmap, - this->grids, - WarpX::grid_type, - this->m_poisson_boundary_handler, - is_solver_igf_on_lev0, - WarpX::do_single_precision_comms, - this->ref_ratio, - post_phi_calculation, - gett_new(0), - eb_farray_box_factory - ); - + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { + m_implicit_darwin_solver->ComputeSigma(); + ablastr::fields::computeSemiImplicitPhi( + sorted_rho, + sorted_phi, + m_implicit_darwin_solver->sigma, + required_precision, + absolute_tolerance, + max_iters, + verbosity, + this->geom, + this->dmap, + this->grids, + this->m_poisson_boundary_handler, + is_solver_igf_on_lev0, + WarpX::do_single_precision_comms, + this->ref_ratio, + post_phi_calculation, + gett_new(0), + eb_farray_box_factory + ); + } else { + ablastr::fields::computePhi( + sorted_rho, + sorted_phi, + beta, + required_precision, + absolute_tolerance, + max_iters, + verbosity, + this->geom, + this->dmap, + this->grids, + this->m_poisson_boundary_handler, + is_solver_igf_on_lev0, + WarpX::do_single_precision_comms, + this->ref_ratio, + post_phi_calculation, + gett_new(0), + eb_farray_box_factory + ); + } } diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H index bc4e3b35ccc..7f86c1dc4da 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H @@ -19,13 +19,9 @@ #include "WarpX.H" #include -#include #include #include -#include -#include -#include /** * \brief This class contains the parameters needed to evaluate the implicit @@ -44,35 +40,9 @@ public: /** Helper function to clear values from multifabs. */ void ClearLevel (int lev); - /** - * \brief - * Function to update the E-field using the implicit Darwin model. - */ - void AddSpaceChargeField ( - amrex::Vector>& rhofield, - amrex::Vector>& phifield, - amrex::Vector, 3>>& Efield, - amrex::Real required_precision, amrex::Real absolute_tolerance, - int max_iters, int verbosity - ); - - template - void ComputePhi ( - const amrex::Vector >& rho, - amrex::Vector >& phi, - amrex::Real relative_tolerance, amrex::Real absolute_tolerance, - int max_iters, int verbosity, - amrex::Geometry const geom, - T_BoundaryHandler const boundary_handler, - std::optional > eb_farray_box_factory = std::nullopt - ) const; - void ComputeSigma (); // Declare multifabs specifically needed for the implicit Darwin model - // amrex::Vector> sigma_fp; - // amrex::Vector, AMREX_SPACEDIM > > sigma; - // amrex::Vector > sigma; amrex::Vector > sigma; amrex::Real m_C_SI = 6.0; diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp index 04f7edb73d7..b6be7edb646 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp @@ -38,18 +38,7 @@ void ImplicitDarwinSolver::AllocateLevelMFs ( "Implicit Darwin solver only support a nodal charge density field at present" ); - // sigma[lev][0] = MultiFab( - // amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), - // dm, ncomps, ngRho - // ); - // sigma[lev][0].setVal(0.0_rt); - // sigma[lev][1] = MultiFab( - // amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), - // dm, ncomps, ngRho - // ); - // sigma[lev][1].setVal(0.0_rt); - - // The "sigma" multifabs stores the dressing of the Poisson equation. It + // The "sigma" multifab stores the dressing of the Poisson equation. It // is a cell-centered multifab. WarpX::AllocInitMultiFab( sigma[lev], amrex::convert(ba, IntVect( AMREX_D_DECL(0,0,0) )), @@ -66,87 +55,6 @@ void ImplicitDarwinSolver::AllocateLevelMFs ( void ImplicitDarwinSolver::ClearLevel (int lev) { sigma[lev].reset(); - // for (int i = 0; i < 3; ++i) { - // sigma[lev][i].reset(); - // } -} - -void -ImplicitDarwinSolver::AddSpaceChargeField ( - amrex::Vector>& rhofield, - amrex::Vector>& phifield, - amrex::Vector, 3>>& Efield, - Real const required_precision, Real absolute_tolerance, - int const max_iters, int const verbosity -) -{ - // Reset all E fields to 0, before calculating space-charge fields - WARPX_PROFILE("WarpX::ComputeSpaceChargeField::reset_fields"); - for (int lev = 0; lev <= 0; lev++) { - for (int comp=0; comp<3; comp++) { - Efield[lev][comp]->setVal(0); - // Bfield_fp[lev][comp]->setVal(0); - } - } - - WARPX_PROFILE("ImplicitDarwinSolver::AddSpaceChargeField"); - auto& warpx = WarpX::GetInstance(); - - // Store the boundary conditions for the field solver if they haven't been - // stored yet - if (!warpx.m_poisson_boundary_handler.bcs_set) - warpx.m_poisson_boundary_handler.definePhiBCs(warpx.Geom(0)); - - // Deposit particle charge density (source of Poisson solver) - auto& mypc = warpx.GetPartContainer(); - mypc.DepositCharge(rhofield, 0.0_rt); - // if (warpx.do_fluid_species) { - // int const lev = 0; - // myfl->DepositCharge( lev, *rho_fp[lev] ); - // } - - // Reflect density over PEC boundaries, if needed. - // filter (if used), exchange guard cells, interpolate across MR levels - // and apply boundary conditions - warpx.SyncCurrentAndRho(); - - // sigma dresses the Poisson equation to filter out plasma frequency effects - ComputeSigma(); - -#if defined(AMREX_USE_EB) - std::optional > eb_farray_box_factory; - amrex::Vector< - amrex::EBFArrayBoxFactory const * - > factories; - for (int lev = 0; lev <= finest_level; ++lev) { - factories.push_back(&warpx.fieldEBFactory(lev)); - } - eb_farray_box_factory = factories; -#else - const std::optional post_phi_calculation; - const std::optional > eb_farray_box_factory; -#endif - - // Use the AMREX MLMG solver - ComputePhi( - rhofield, phifield, required_precision, - absolute_tolerance, max_iters, verbosity, - warpx.Geom(0), - warpx.m_poisson_boundary_handler, - eb_farray_box_factory - ); - - // Compute the electric field. Note that if an EB is used the electric - // field will be calculated in the computePhi call above. -#ifndef AMREX_USE_EB - // beta is the moving frame velocity i.e. zero in lab frame - const std::array beta = {0._rt}; - - warpx.computeE( Efield, phifield, beta ); -#endif - -// // Compute the magnetic field -// computeB( Bfield_fp, phi_fp, beta ); } void ImplicitDarwinSolver::ComputeSigma () { @@ -209,114 +117,3 @@ void ImplicitDarwinSolver::ComputeSigma () { } } } - -template -void ImplicitDarwinSolver::ComputePhi ( - const amrex::Vector >& rho, - amrex::Vector >& phi, - Real const relative_tolerance, - Real absolute_tolerance, - int const max_iters, - int const verbosity, - amrex::Geometry const geom, - T_BoundaryHandler const boundary_handler, - std::optional > eb_farray_box_factory ) const -{ - using namespace amrex::literals; - - int const lev = 0; - - // scale rho appropriately; also determine if rho is zero everywhere - amrex::Real max_norm_b = 0.0; - rho[lev]->mult(-1._rt/ablastr::constant::SI::ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! - max_norm_b = amrex::max(max_norm_b, rho[lev]->norm0()); - - amrex::ParallelDescriptor::ReduceRealMax(max_norm_b); - - const bool always_use_bnorm = (max_norm_b > 0); - if (!always_use_bnorm) { - if (absolute_tolerance == 0.0) absolute_tolerance = amrex::Real(1e-6); - ablastr::warn_manager::WMRecordWarning( - "ImplicitDarwinSolver", - "Max norm of rho is 0", - ablastr::warn_manager::WarnPriority::low - ); - } - - // auto& warpx = WarpX::GetInstance(); - // amrex::LPInfo info = LPInfo().setMaxCoarseningLevel(0); - amrex::LPInfo info; - -#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) - // In the presence of EB or RZ: the solver assumes that the beam is - // propagating along one of the axes of the grid, i.e. that only *one* - // of the components of `beta` is non-negligible. - amrex::MLEBNodeLaplacian linop( {geom}, {rho.boxArray()}, {rho[lev].dm()}, info -#if defined(AMREX_USE_EB) - , {eb_farray_box_factory.value()[lev]} -#endif - ); - - // Note: this assumes that the beam is propagating along - // one of the axes of the grid, i.e. that only *one* of the - // components of `beta` is non-negligible. // we use this -#if defined(WARPX_DIM_RZ) - linop.setSigma({0._rt, 1._rt-beta_solver[1]*beta_solver[1]}); -#else - linop.setSigma({AMREX_D_DECL( - 1._rt-beta_solver[0]*beta_solver[0], - 1._rt-beta_solver[1]*beta_solver[1], - 1._rt-beta_solver[2]*beta_solver[2])}); - // linop.setScalars(0.0, 1.0); - // linop.setBCoeffs(lev, amrex::GetArrOfConstPtrs(sigma[lev]) ); // for the non-axis-aligned solver -#endif - -#if defined(AMREX_USE_EB) - // if the EB potential only depends on time, the potential can be passed - // as a float instead of a callable - if (boundary_handler.phi_EB_only_t) { - linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value())); - } - else - linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); -#endif -#else - // In the absence of EB and RZ: use a more generic solver - // that can handle sigma in all directions - amrex::MLNodeLaplacian linop( - {geom}, {rho[lev]->boxArray()}, {rho[lev]->DistributionMap()}, info - // {}, 1.0 - ); - linop.setSigma(lev, *sigma[lev]); -#endif - - // Solve the Poisson equation - linop.setDomainBC( boundary_handler.lobc, boundary_handler.hibc ); -#ifdef WARPX_DIM_RZ - linop.setRZ(true); -#endif - amrex::MLMG mlmg(linop); // actual solver defined here - mlmg.setVerbose(verbosity); - mlmg.setMaxIter(max_iters); - mlmg.setAlwaysUseBNorm(always_use_bnorm); - - // // create a vector to our fields, sorted by level - // amrex::Vector sorted_rho; - // amrex::Vector sorted_phi; - // // for (int nlev = 0; nlev <= finest_level; ++nlev) { - // sorted_rho.emplace_back(rho[lev].get()); - // sorted_phi.emplace_back(phi[lev].get()); - // // } - - // Solve Poisson equation at lev - mlmg.solve( - { phi[lev].get() }, { rho[lev].get() }, // GetVecOfPtrs(phi), GetVecOfConstPtrs(rho) // - relative_tolerance, absolute_tolerance - ); - - // Run additional operations, such as calculation of the E field for embedded boundaries - // if constexpr (!std::is_same::value) - // if (post_phi_calculation.has_value()) - // post_phi_calculation.value()(mlmg, lev); - -} diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H new file mode 100644 index 00000000000..177e9f25a21 --- /dev/null +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -0,0 +1,303 @@ +/* Copyright 2024 Roelof Groenewald (TAE Technologies) + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H +#define ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H + +#include +#include +#include +#include +#include +#include +#include + +#if defined(WARPX_USE_FFT) && defined(WARPX_DIM_3D) +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) +# include +#endif +#ifdef AMREX_USE_EB +# include +#endif + +#include +#include + + +namespace ablastr::fields { + +/** Compute the potential `phi` by solving the Poisson equation + * + * Uses `rho` as a source. This uses the AMReX solver. + * + * More specifically, this solves the equation + * \f[ + * ... \phi = -\frac{r \rho}{\epsilon_0} + * \f] + * + * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler + * \tparam T_PostPhiCalculationFunctor a calculation per level directly after phi was calculated + * \tparam T_FArrayBoxFactory usually nothing or an amrex::EBFArrayBoxFactory (EB ONLY) + * \param[in] rho The charge density a given species + * \param[out] phi The potential to be computed by this function + * \param[in] sigma The semi-implicit matrix + * \param[in] relative_tolerance The relative convergence threshold for the MLMG solver + * \param[in] absolute_tolerance The absolute convergence threshold for the MLMG solver + * \param[in] max_iters The maximum number of iterations allowed for the MLMG solver + * \param[in] verbosity The verbosity setting for the MLMG solver + * \param[in] geom the geometry per level (e.g., from AmrMesh) + * \param[in] dmap the distribution mapping per level (e.g., from AmrMesh) + * \param[in] grids the grids per level (e.g., from AmrMesh) + * \param[in] boundary_handler a handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler + * \param[in] is_solver_igf_on_lev0 boolean to select the Poisson solver: 1 for FFT on level 0 & Multigrid on other levels, 0 for Multigrid on all levels + * \param[in] do_single_precision_comms perform communications in single precision + * \param[in] rel_ref_ratio mesh refinement ratio between levels (default: 1) + * \param[in] post_phi_calculation perform a calculation per level directly after phi was calculated; required for embedded boundaries (default: none) + * \param[in] current_time the current time; required for embedded boundaries (default: none) + * \param[in] eb_farray_box_factory a factory for field data, @see amrex::EBFArrayBoxFactory; required for embedded boundaries (default: none) + */ +template< + typename T_BoundaryHandler, + typename T_PostPhiCalculationFunctor = std::nullopt_t, + typename T_FArrayBoxFactory = void +> +void +computeSemiImplicitPhi (amrex::Vector const & rho, + amrex::Vector & phi, + amrex::Vector> const & sigma, + amrex::Real const relative_tolerance, + amrex::Real absolute_tolerance, + int const max_iters, + int const verbosity, + amrex::Vector const& geom, + amrex::Vector const& dmap, + amrex::Vector const& grids, + T_BoundaryHandler const boundary_handler, + bool is_solver_igf_on_lev0, + bool const do_single_precision_comms = false, + std::optional > rel_ref_ratio = std::nullopt, + [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation = std::nullopt, + [[maybe_unused]] std::optional current_time = std::nullopt, // only used for EB + [[maybe_unused]] std::optional > eb_farray_box_factory = std::nullopt // only used for EB +) +{ + using namespace amrex::literals; + + ABLASTR_PROFILE("computeSemiImplicitPhi"); + + if (!rel_ref_ratio.has_value()) { + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(rho.size() == 1u, + "rel_ref_ratio must be set if mesh-refinement is used"); + rel_ref_ratio = amrex::Vector{{amrex::IntVect(AMREX_D_DECL(1, 1, 1))}}; + } + + auto const finest_level = static_cast(rho.size() - 1); + + // determine if rho is zero everywhere + amrex::Real max_norm_b = 0.0; + for (int lev=0; lev<=finest_level; lev++) { + max_norm_b = amrex::max(max_norm_b, rho[lev]->norm0()); + } + amrex::ParallelDescriptor::ReduceRealMax(max_norm_b); + + const bool always_use_bnorm = (max_norm_b > 0); + if (!always_use_bnorm) { + if (absolute_tolerance == 0.0) { absolute_tolerance = amrex::Real(1e-6); } + ablastr::warn_manager::WMRecordWarning( + "ElectrostaticSolver", + "Max norm of rho is 0", + ablastr::warn_manager::WarnPriority::low + ); + } + +#if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) + amrex::LPInfo info; +#else + const amrex::LPInfo info; +#endif + + for (int lev=0; lev<=finest_level; lev++) { + +#if !defined(WARPX_USE_FFT) + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, + "Must compile with -DWarpX_FFT=ON to use the FFT solver!"); +#endif + +#if !defined(WARPX_DIM_3D) + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, + "The FFT Poisson solver is currently only implemented for 3D!"); +#endif + +#if (defined(WARPX_USE_FFT) && defined(WARPX_DIM_3D)) + amrex::Abort("Cannot use SI solver with WARPX_USE_FFT=ON."); + // Use the Integrated Green Function solver (FFT) on the coarsest level if it was selected + // if(is_solver_igf_on_lev0 && lev==0){ + // amrex::Array const dx_igf + // {AMREX_D_DECL(geom[lev].CellSize(0)/std::sqrt(1._rt-beta_solver[0]*beta_solver[0]), + // geom[lev].CellSize(1)/std::sqrt(1._rt-beta_solver[1]*beta_solver[1]), + // geom[lev].CellSize(2)/std::sqrt(1._rt-beta_solver[2]*beta_solver[2]))}; + // if ( max_norm_b == 0 ) { + // phi[lev]->setVal(0); + // } else { + // computePhiIGF( *rho[lev], *phi[lev], dx_igf, grids[lev] ); + // } + // continue; + // } +#endif + + // Use the Multigrid (MLMG) solver if selected or on refined patches + // but first scale rho appropriately + using namespace ablastr::constant::SI; + rho[lev]->mult(-1._rt/ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! + +#if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) + // Determine whether to use semi-coarsening + amrex::Array dx_scaled + {AMREX_D_DECL(geom[lev].CellSize(0), geom[lev].CellSize(1), geom[lev].CellSize(2))}; + int max_semicoarsening_level = 0; + int semicoarsening_direction = -1; + const auto min_dir = static_cast(std::distance(dx_scaled.begin(), + std::min_element(dx_scaled.begin(),dx_scaled.end()))); + const auto max_dir = static_cast(std::distance(dx_scaled.begin(), + std::max_element(dx_scaled.begin(),dx_scaled.end()))); + if (dx_scaled[max_dir] > dx_scaled[min_dir]) { + semicoarsening_direction = max_dir; + max_semicoarsening_level = static_cast + (std::log2(dx_scaled[max_dir]/dx_scaled[min_dir])); + } + if (max_semicoarsening_level > 0) { + info.setSemicoarsening(true); + info.setMaxSemicoarseningLevel(max_semicoarsening_level); + info.setSemicoarseningDirection(semicoarsening_direction); + } +#endif + +#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) + // In the presence of EB or RZ: use the EB enabled solver. + amrex::MLEBNodeFDLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info +#if defined(AMREX_USE_EB) + , {eb_farray_box_factory.value()[lev]} +#endif + ); + +#if defined(AMREX_USE_EB) + // if the EB potential only depends on time, the potential can be passed + // as a float instead of a callable + if (boundary_handler.phi_EB_only_t) { + linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value())); + } + else + linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); +#endif +#else + // In the absence of EB and RZ: use a more generic solver + amrex::MLNodeLaplacian linop( {geom[lev]}, {grids[lev]}, + {dmap[lev]}, info ); +#endif + linop.setSigma(lev, *sigma[lev]); + + // Solve the Poisson equation + linop.setDomainBC( boundary_handler.lobc, boundary_handler.hibc ); +#ifdef WARPX_DIM_RZ + linop.setRZ(true); +#endif + amrex::MLMG mlmg(linop); // actual solver defined here + mlmg.setVerbose(verbosity); + mlmg.setMaxIter(max_iters); + mlmg.setAlwaysUseBNorm(always_use_bnorm); + + // Solve Poisson equation at lev + mlmg.solve( {phi[lev]}, {rho[lev]}, + relative_tolerance, absolute_tolerance ); + + // needed for solving the levels by levels: + // - coarser level is initial guess for finer level + // - coarser level provides boundary values for finer level patch + // Interpolation from phi[lev] to phi[lev+1] + // (This provides both the boundary conditions and initial guess for phi[lev+1]) + if (lev < finest_level) { + + // Allocate phi_cp for lev+1 + amrex::BoxArray ba = phi[lev+1]->boxArray(); + const amrex::IntVect& refratio = rel_ref_ratio.value()[lev]; + ba.coarsen(refratio); + const int ncomp = linop.getNComp(); + amrex::MultiFab phi_cp(ba, phi[lev+1]->DistributionMap(), ncomp, 1); + + // Copy from phi[lev] to phi_cp (in parallel) + const amrex::IntVect& ng = amrex::IntVect::TheUnitVector(); + const amrex::Periodicity& crse_period = geom[lev].periodicity(); + + ablastr::utils::communication::ParallelCopy( + phi_cp, + *phi[lev], + 0, + 0, + 1, + ng, + ng, + do_single_precision_comms, + crse_period + ); + + // Local interpolation from phi_cp to phi[lev+1] +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(*phi[lev + 1], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + amrex::Array4 const phi_fp_arr = phi[lev + 1]->array(mfi); + amrex::Array4 const phi_cp_arr = phi_cp.array(mfi); + + details::PoissonInterpCPtoFP const interp(phi_fp_arr, phi_cp_arr, refratio); + + amrex::Box const b = mfi.tilebox(phi[lev + 1]->ixType().toIntVect()); + amrex::ParallelFor(b, interp); + } + + } + + // Run additional operations, such as calculation of the E field for embedded boundaries + if constexpr (!std::is_same::value) { + if (post_phi_calculation.has_value()) { + post_phi_calculation.value()(mlmg, lev); + } + } + + } // loop over lev(els) +} + +} // namespace ablastr::fields + +#endif // ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H From c29a6122532e71e14c27a996e254f517b9a888c1 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 22 Jul 2024 21:54:44 -0700 Subject: [PATCH 06/28] appropriately name new semi-implicit solver --- Python/pywarpx/picmi.py | 6 ++--- Source/Diagnostics/Diagnostics.cpp | 2 +- Source/Evolve/WarpXEvolve.cpp | 2 +- Source/FieldSolver/ElectrostaticSolver.cpp | 12 +++++----- .../ElectrostaticSolver/CMakeLists.txt | 2 +- .../ImplicitDarwinSolver_fwd.H | 15 ------------ .../ElectrostaticSolver/Make.package | 2 +- ...citDarwinSolver.H => SemiImplicitSolver.H} | 12 +++++----- ...arwinSolver.cpp => SemiImplicitSolver.cpp} | 12 +++++----- .../SemiImplicitSolver_fwd.H | 15 ++++++++++++ Source/Initialization/WarpXInitData.cpp | 4 ++-- Source/Python/WarpX.cpp | 2 +- Source/Utils/WarpXAlgorithmSelection.H | 2 +- Source/Utils/WarpXAlgorithmSelection.cpp | 2 +- Source/WarpX.H | 6 ++--- Source/WarpX.cpp | 24 +++++++++---------- 16 files changed, 60 insertions(+), 60 deletions(-) delete mode 100644 Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H rename Source/FieldSolver/ElectrostaticSolver/{ImplicitDarwinSolver.H => SemiImplicitSolver.H} (83%) rename Source/FieldSolver/ElectrostaticSolver/{ImplicitDarwinSolver.cpp => SemiImplicitSolver.cpp} (92%) create mode 100644 Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver_fwd.H diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 2c6fe8b3bb0..cca6de13e1d 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1551,7 +1551,7 @@ def init(self, kw): self.absolute_tolerance = kw.pop('warpx_absolute_tolerance', None) self.self_fields_verbosity = kw.pop('warpx_self_fields_verbosity', None) self.magnetostatic = kw.pop('warpx_magnetostatic', False) - self.darwin_labframe = kw.pop('warpx_darwin_labframe', False) + self.semi_implicit = kw.pop('warpx_semi_implicit', False) self.semi_implicit_factor = kw.pop('warpx_semi_implicit_factor', None) def solver_initialize_inputs(self): @@ -1566,8 +1566,8 @@ def solver_initialize_inputs(self): else: if self.magnetostatic: pywarpx.warpx.do_electrostatic = 'labframe-electromagnetostatic' - elif self.darwin_labframe: - pywarpx.warpx.do_electrostatic = 'labframe-implicit' + elif self.semi_implicit: + pywarpx.warpx.do_electrostatic = 'labframe-semi-implicit' pywarpx.warpx.semi_implicit_factor = self.semi_implicit_factor else: pywarpx.warpx.do_electrostatic = 'labframe' diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index 80a61564148..cb11f8782b4 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -77,7 +77,7 @@ Diagnostics::BaseReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE( warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrame || warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameDarwinImplicit, + warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameSemiImplicit, "plot phi only works if do_electrostatic = labframe or do_electrostatic = labframe-electromagnetostatic"); } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b557836d878..24a165853b3 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -238,7 +238,7 @@ WarpX::Evolve (int numsteps) // loop (i.e. immediately after a `Redistribute` and before particle // positions are next pushed) so that the particles do not deposit out of bounds // and so that the fields are at the correct time in the output. - bool const reset_fields = (electrostatic_solver_id != ElectrostaticSolverAlgo::LabFrameDarwinImplicit); + bool const reset_fields = (electrostatic_solver_id != ElectrostaticSolverAlgo::LabFrameSemiImplicit); ComputeSpaceChargeField( reset_fields ); if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { // Call Magnetostatic Solver to solve for the vector potential A and compute the diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index a5c803f45de..d0774a57b57 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -7,7 +7,7 @@ #include "WarpX.H" #include "FieldSolver/ElectrostaticSolver.H" -#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H" +#include "FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H" #include "Fluids/MultiFluidContainer.H" #include "Fluids/WarpXFluidContainer.H" #include "Parallelization/GuardCellManager.H" @@ -76,7 +76,7 @@ WarpX::ComputeSpaceChargeField (bool const reset_fields) if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit ) { + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit ) { AddSpaceChargeFieldLabFrame(); } else { @@ -279,7 +279,7 @@ WarpX::AddSpaceChargeFieldLabFrame () } else { #if defined(WARPX_DIM_1D_Z) - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { amrex::Abort("Cannot use SI solver in 1d."); } // Use the tridiag solver with 1D @@ -391,12 +391,12 @@ WarpX::computePhi (const amrex::Vector >& rho, bool const is_solver_igf_on_lev0 = WarpX::poisson_solver_id == PoissonSolverAlgo::IntegratedGreenFunction; - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) { - m_implicit_darwin_solver->ComputeSigma(); + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { + m_semi_implicit_solver->ComputeSigma(); ablastr::fields::computeSemiImplicitPhi( sorted_rho, sorted_phi, - m_implicit_darwin_solver->sigma, + m_semi_implicit_solver->sigma, required_precision, absolute_tolerance, max_iters, diff --git a/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt b/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt index c3596fe4fdd..6486f3e96a0 100644 --- a/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt +++ b/Source/FieldSolver/ElectrostaticSolver/CMakeLists.txt @@ -2,6 +2,6 @@ foreach(D IN LISTS WarpX_DIMS) warpx_set_suffix_dims(SD ${D}) target_sources(lib_${SD} PRIVATE - ImplicitDarwinSolver.cpp + SemiImplicitSolver.cpp ) endforeach() diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H b/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H deleted file mode 100644 index 80b0ad22858..00000000000 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H +++ /dev/null @@ -1,15 +0,0 @@ -/* Copyright 2023 The WarpX Community - * - * This file is part of WarpX. - * - * Authors: Roelof Groenewald (TAE Technologies) - * - * License: BSD-3-Clause-LBNL - */ - -#ifndef WARPX_IMPLICITDARWINSOLVER_FWD_H -#define WARPX_IMPLICITDARWINSOLVER_FWD_H - -class ImplicitDarwinSolver; - -#endif /* WARPX_IMPLICITDARWINSOLVER_FWD_H */ diff --git a/Source/FieldSolver/ElectrostaticSolver/Make.package b/Source/FieldSolver/ElectrostaticSolver/Make.package index c087c171c05..19e1f13d8d4 100644 --- a/Source/FieldSolver/ElectrostaticSolver/Make.package +++ b/Source/FieldSolver/ElectrostaticSolver/Make.package @@ -1,3 +1,3 @@ -CEXE_sources += ImplicitDarwinSolver.cpp +CEXE_sources += SemiImplicitSolver.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/ElectrostaticSolver diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H similarity index 83% rename from Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H rename to Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H index 7f86c1dc4da..fdbebe08aef 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H @@ -7,10 +7,10 @@ * License: BSD-3-Clause-LBNL */ -#ifndef WARPX_IMPLICITDARWINSOLVER_H_ -#define WARPX_IMPLICITDARWINSOLVER_H_ +#ifndef WARPX_SEMIIMPLICITSOLVER_H_ +#define WARPX_SEMIIMPLICITSOLVER_H_ -#include "ImplicitDarwinSolver_fwd.H" +#include "SemiImplicitSolver_fwd.H" #include "Particles/MultiParticleContainer.H" #include "Utils/Parser/ParserUtils.H" @@ -27,10 +27,10 @@ * \brief This class contains the parameters needed to evaluate the implicit * electrostatic Darwin model. */ -class ImplicitDarwinSolver +class SemiImplicitSolver { public: - ImplicitDarwinSolver (int nlevs_max); // constructor + SemiImplicitSolver (int nlevs_max); // constructor /** Allocate multifabs needed for the Darwin model. Called in constructor. */ void AllocateMFs (int nlevs_max); @@ -48,4 +48,4 @@ public: amrex::Real m_C_SI = 6.0; }; -#endif // WARPX_IMPLICITDARWINSOLVER_H_ +#endif // WARPX_SEMIIMPLICITSOLVER_H_ diff --git a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp similarity index 92% rename from Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp rename to Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp index b6be7edb646..79009356d13 100644 --- a/Source/FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp @@ -7,11 +7,11 @@ * License: BSD-3-Clause-LBNL */ -#include "ImplicitDarwinSolver.H" +#include "SemiImplicitSolver.H" using namespace amrex; -ImplicitDarwinSolver::ImplicitDarwinSolver ( int nlevs_max ) +SemiImplicitSolver::SemiImplicitSolver ( int nlevs_max ) { AllocateMFs(nlevs_max); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -24,12 +24,12 @@ ImplicitDarwinSolver::ImplicitDarwinSolver ( int nlevs_max ) utils::parser::queryWithParser(pp_warpx, "semi_implicit_factor", m_C_SI); } -void ImplicitDarwinSolver::AllocateMFs (int nlevs_max) +void SemiImplicitSolver::AllocateMFs (int nlevs_max) { sigma.resize(nlevs_max); } -void ImplicitDarwinSolver::AllocateLevelMFs ( +void SemiImplicitSolver::AllocateLevelMFs ( int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, int ncomps, const amrex::IntVect& ngRho, const amrex::IntVect& rho_nodal_flag) { @@ -52,12 +52,12 @@ void ImplicitDarwinSolver::AllocateLevelMFs ( #endif } -void ImplicitDarwinSolver::ClearLevel (int lev) +void SemiImplicitSolver::ClearLevel (int lev) { sigma[lev].reset(); } -void ImplicitDarwinSolver::ComputeSigma () { +void SemiImplicitSolver::ComputeSigma () { int const lev = 0; sigma[lev]->setVal(1.0_rt); diff --git a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver_fwd.H b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver_fwd.H new file mode 100644 index 00000000000..8ef3b73af18 --- /dev/null +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver_fwd.H @@ -0,0 +1,15 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Roelof Groenewald (TAE Technologies) + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_SEMIIMPLICITSOLVER_FWD_H +#define WARPX_SEMIIMPLICITSOLVER_FWD_H + +class SemiImplicitSolver; + +#endif /* WARPX_SEMIIMPLICITSOLVER_FWD_H */ diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index bf636df01bf..1018cf6297e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -267,9 +267,9 @@ WarpX::PrintMainPICparameters () amrex::Print() << "Operation mode: | Electrostatic" << "\n"; amrex::Print() << " | - laboratory frame, electrostatic + magnetostatic" << "\n"; } - else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit){ + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit){ amrex::Print() << "Operation mode: | Electrostatic" << "\n"; - amrex::Print() << " | - laboratory frame, Darwin implicit scheme" << "\n"; + amrex::Print() << " | - laboratory frame, semi-implicit scheme" << "\n"; } else{ amrex::Print() << "Operation mode: | Electromagnetic" << "\n"; diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp index 57bde115af6..4d314689835 100644 --- a/Source/Python/WarpX.cpp +++ b/Source/Python/WarpX.cpp @@ -11,7 +11,7 @@ #include #include #include -#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H" +#include "FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H" #include #include #include diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 06429e64ca5..e1bea71738c 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -68,7 +68,7 @@ struct ElectrostaticSolverAlgo { Relativistic = 1, LabFrameElectroMagnetostatic = 2, LabFrame = 3, // Non relativistic - LabFrameDarwinImplicit = 4 + LabFrameSemiImplicit = 4 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index c392defd35c..7942e8f4f46 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -53,7 +53,7 @@ const std::map electrostatic_solver_algo_to_int = { {"relativistic", ElectrostaticSolverAlgo::Relativistic}, {"labframe-electromagnetostatic", ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic}, {"labframe", ElectrostaticSolverAlgo::LabFrame}, - {"labframe-implicit", ElectrostaticSolverAlgo::LabFrameDarwinImplicit}, + {"labframe-semi-implicit", ElectrostaticSolverAlgo::LabFrameSemiImplicit}, {"default", ElectrostaticSolverAlgo::None } }; diff --git a/Source/WarpX.H b/Source/WarpX.H index 0c1cefbb9e6..d38193952ba 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -16,7 +16,7 @@ #include "Diagnostics/MultiDiagnostics_fwd.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags_fwd.H" #include "EmbeddedBoundary/WarpXFaceInfoBox_fwd.H" -#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver_fwd.H" +#include "FieldSolver/ElectrostaticSolver/SemiImplicitSolver_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties_fwd.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel_fwd.H" @@ -1636,8 +1636,8 @@ private: // Hybrid PIC algorithm parameters std::unique_ptr m_hybrid_pic_model; - // Implicit Darwin electrostatic algorithm - std::unique_ptr m_implicit_darwin_solver; + // Semi-implicit electrostatic algorithm + std::unique_ptr m_semi_implicit_solver; // Load balancing /** Load balancing intervals that reads the "load_balance_intervals" string int the input file diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index cd0562073b9..c4bbf2ddfc2 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -15,7 +15,7 @@ #include "Diagnostics/MultiDiagnostics.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "EmbeddedBoundary/WarpXFaceInfoBox.H" -#include "FieldSolver/ElectrostaticSolver/ImplicitDarwinSolver.H" +#include "FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H" @@ -334,10 +334,10 @@ WarpX::WarpX () vector_potential_grad_buf_b_stag.resize(nlevs_max); } - // Initialize the implicit Darwin electrostatic solver if required - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + // Initialize the semi-implicit electrostatic solver if required + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { - m_implicit_darwin_solver = std::make_unique(nlevs_max); + m_semi_implicit_solver = std::make_unique(nlevs_max); } if (fft_do_time_averaging) @@ -751,7 +751,7 @@ WarpX::ReadParameters () if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { // Note that with the relativistic version, these parameters would be // input for each species. @@ -2105,9 +2105,9 @@ WarpX::ClearLevel (int lev) current_buf[lev][i].reset(); } - if (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + if (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { - m_implicit_darwin_solver->ClearLevel(lev); + m_semi_implicit_solver->ClearLevel(lev); } if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) @@ -2370,10 +2370,10 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm dm, ncomps, ngEB, lev, "vector_potential_grad_buf_b_stag[z]", 0.0_rt); } - // Allocate extra multifabs needed by the implicit Darwin electrostatic algorithm. - if (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) + // Allocate extra multifabs needed by the semi-implicit electrostatic algorithm. + if (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { - m_implicit_darwin_solver->AllocateLevelMFs( + m_semi_implicit_solver->AllocateLevelMFs( lev, ba, dm, ncomps, ngRho, rho_nodal_flag ); } @@ -2462,7 +2462,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm int rho_ncomps = 0; if( (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) || (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) || - (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit) || + (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) || (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) ) { rho_ncomps = ncomps; } @@ -2482,7 +2482,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameDarwinImplicit ) + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit ) { const IntVect ngPhi = IntVect( AMREX_D_DECL(1,1,1) ); AllocInitMultiFab(phi_fp[lev], amrex::convert(ba, phi_nodal_flag), dm, ncomps, ngPhi, lev, "phi_fp", 0.0_rt); From c0aac0398a709c1ffae7fff22afa7044c3cc15e7 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Tue, 23 Jul 2024 11:38:34 -0700 Subject: [PATCH 07/28] use AMReX to calculate E-field directly --- Source/FieldSolver/ElectrostaticSolver.cpp | 40 +++++++-- .../ElectrostaticSolver/SemiImplicitSolver.H | 8 +- .../SemiImplicitSolver.cpp | 6 +- .../fields/SemiImplicitPoissonSolver.H | 89 ++++++------------- 4 files changed, 68 insertions(+), 75 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index d0774a57b57..f2bef2b37b8 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -279,9 +279,8 @@ WarpX::AddSpaceChargeFieldLabFrame () } else { #if defined(WARPX_DIM_1D_Z) - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { - amrex::Abort("Cannot use SI solver in 1d."); - } + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( electrostatic_solver_id != ElectrostaticSolverAlgo::LabFrameSemiImplicit, + "Cannot use SI solver in 1d."); // Use the tridiag solver with 1D computePhiTriDiagonal(rho_fp, phi_fp); #else @@ -339,14 +338,14 @@ WarpX::computePhi (const amrex::Vector >& rho, sorted_phi.emplace_back(phi[lev].get()); } -#if defined(AMREX_USE_EB) - std::optional post_phi_calculation; +#if defined(AMREX_USE_EB) // EB: use AMReX to directly calculate the electric field since with EB's the // simple finite difference scheme in WarpX::computeE sometimes fails if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) { // TODO: maybe make this a helper function or pass Efield_fp directly amrex::Vector< @@ -384,7 +383,34 @@ WarpX::computePhi (const amrex::Vector >& rho, } eb_farray_box_factory = factories; #else - const std::optional post_phi_calculation; + if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) + { + // TODO: maybe make this a helper function or pass Efield_fp directly + amrex::Vector< + amrex::Array + > e_field; + for (int lev = 0; lev <= finest_level; ++lev) { + e_field.push_back( +# if defined(WARPX_DIM_1D_Z) + amrex::Array{ + getFieldPointer(FieldType::Efield_fp, lev, 2) + } +# elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::Array{ + getFieldPointer(FieldType::Efield_fp, lev, 0), + getFieldPointer(FieldType::Efield_fp, lev, 2) + } +# elif defined(WARPX_DIM_3D) + amrex::Array{ + getFieldPointer(FieldType::Efield_fp, lev, 0), + getFieldPointer(FieldType::Efield_fp, lev, 1), + getFieldPointer(FieldType::Efield_fp, lev, 2) + } +# endif + ); + } + post_phi_calculation = ElectrostaticSolver::EBCalcEfromPhiPerLevel(e_field); + } const std::optional > eb_farray_box_factory; #endif diff --git a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H index fdbebe08aef..45b96b7377d 100644 --- a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H @@ -24,15 +24,15 @@ /** - * \brief This class contains the parameters needed to evaluate the implicit - * electrostatic Darwin model. + * \brief This class contains the parameters needed to evaluate the semi-implicit + * electrostatic model. */ class SemiImplicitSolver { public: SemiImplicitSolver (int nlevs_max); // constructor - /** Allocate multifabs needed for the Darwin model. Called in constructor. */ + /** Allocate multifabs needed for the semi-implicit solver. Called in constructor. */ void AllocateMFs (int nlevs_max); void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, int ncomps, const amrex::IntVect& ngRho, const amrex::IntVect& rho_nodal_flag); @@ -42,7 +42,7 @@ public: void ComputeSigma (); - // Declare multifabs specifically needed for the implicit Darwin model + // Declare multifabs specifically needed for the semi-implicit solver amrex::Vector > sigma; amrex::Real m_C_SI = 6.0; diff --git a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp index 79009356d13..3256341e801 100644 --- a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp @@ -16,7 +16,7 @@ SemiImplicitSolver::SemiImplicitSolver ( int nlevs_max ) AllocateMFs(nlevs_max); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (nlevs_max == 1), - "Implicit Darwin solver only support one level at present" + "Semi-implicit electrostatic solver only supports one level at present" ); // read input parameters @@ -35,7 +35,7 @@ void SemiImplicitSolver::AllocateLevelMFs ( { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (rho_nodal_flag == IntVect::TheNodeVector()), - "Implicit Darwin solver only support a nodal charge density field at present" + "Semi-implicit electrostatic solver only supports a nodal charge density field" ); // The "sigma" multifab stores the dressing of the Poisson equation. It @@ -48,7 +48,7 @@ void SemiImplicitSolver::AllocateLevelMFs ( #ifdef WARPX_DIM_RZ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (ncomps == 1), - "Implicit Darwin solver only support m = 0 azimuthal mode at present."); + "Semi-implicit electrostatic solver only supports m = 0 azimuthal mode at present."); #endif } diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index 177e9f25a21..ae74de88233 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -45,9 +45,7 @@ #include #include #include -#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) -# include -#endif +#include #ifdef AMREX_USE_EB # include #endif @@ -148,64 +146,38 @@ computeSemiImplicitPhi (amrex::Vector const & rho, const amrex::LPInfo info; #endif - for (int lev=0; lev<=finest_level; lev++) { - -#if !defined(WARPX_USE_FFT) - ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, - "Must compile with -DWarpX_FFT=ON to use the FFT solver!"); -#endif - -#if !defined(WARPX_DIM_3D) - ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, - "The FFT Poisson solver is currently only implemented for 3D!"); -#endif + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, + "FFT solver cannot be used with semi-implicit Poisson solve"); -#if (defined(WARPX_USE_FFT) && defined(WARPX_DIM_3D)) - amrex::Abort("Cannot use SI solver with WARPX_USE_FFT=ON."); - // Use the Integrated Green Function solver (FFT) on the coarsest level if it was selected - // if(is_solver_igf_on_lev0 && lev==0){ - // amrex::Array const dx_igf - // {AMREX_D_DECL(geom[lev].CellSize(0)/std::sqrt(1._rt-beta_solver[0]*beta_solver[0]), - // geom[lev].CellSize(1)/std::sqrt(1._rt-beta_solver[1]*beta_solver[1]), - // geom[lev].CellSize(2)/std::sqrt(1._rt-beta_solver[2]*beta_solver[2]))}; - // if ( max_norm_b == 0 ) { - // phi[lev]->setVal(0); - // } else { - // computePhiIGF( *rho[lev], *phi[lev], dx_igf, grids[lev] ); - // } - // continue; - // } -#endif + for (int lev=0; lev<=finest_level; lev++) { - // Use the Multigrid (MLMG) solver if selected or on refined patches - // but first scale rho appropriately + // Use the Multigrid (MLMG) solver but first scale rho appropriately using namespace ablastr::constant::SI; rho[lev]->mult(-1._rt/ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! -#if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) - // Determine whether to use semi-coarsening - amrex::Array dx_scaled - {AMREX_D_DECL(geom[lev].CellSize(0), geom[lev].CellSize(1), geom[lev].CellSize(2))}; - int max_semicoarsening_level = 0; - int semicoarsening_direction = -1; - const auto min_dir = static_cast(std::distance(dx_scaled.begin(), - std::min_element(dx_scaled.begin(),dx_scaled.end()))); - const auto max_dir = static_cast(std::distance(dx_scaled.begin(), - std::max_element(dx_scaled.begin(),dx_scaled.end()))); - if (dx_scaled[max_dir] > dx_scaled[min_dir]) { - semicoarsening_direction = max_dir; - max_semicoarsening_level = static_cast - (std::log2(dx_scaled[max_dir]/dx_scaled[min_dir])); - } - if (max_semicoarsening_level > 0) { - info.setSemicoarsening(true); - info.setMaxSemicoarseningLevel(max_semicoarsening_level); - info.setSemicoarseningDirection(semicoarsening_direction); - } -#endif - -#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) - // In the presence of EB or RZ: use the EB enabled solver. +// #if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) +// // Determine whether to use semi-coarsening +// amrex::Array dx_scaled +// {AMREX_D_DECL(geom[lev].CellSize(0), geom[lev].CellSize(1), geom[lev].CellSize(2))}; +// int max_semicoarsening_level = 0; +// int semicoarsening_direction = -1; +// const auto min_dir = static_cast(std::distance(dx_scaled.begin(), +// std::min_element(dx_scaled.begin(),dx_scaled.end()))); +// const auto max_dir = static_cast(std::distance(dx_scaled.begin(), +// std::max_element(dx_scaled.begin(),dx_scaled.end()))); +// if (dx_scaled[max_dir] > dx_scaled[min_dir]) { +// semicoarsening_direction = max_dir; +// max_semicoarsening_level = static_cast +// (std::log2(dx_scaled[max_dir]/dx_scaled[min_dir])); +// } +// if (max_semicoarsening_level > 0) { +// info.setSemicoarsening(true); +// info.setMaxSemicoarseningLevel(max_semicoarsening_level); +// info.setSemicoarseningDirection(semicoarsening_direction); +// } +// #endif + + // Use the EB enabled solver. amrex::MLEBNodeFDLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info #if defined(AMREX_USE_EB) , {eb_farray_box_factory.value()[lev]} @@ -220,11 +192,6 @@ computeSemiImplicitPhi (amrex::Vector const & rho, } else linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); -#endif -#else - // In the absence of EB and RZ: use a more generic solver - amrex::MLNodeLaplacian linop( {geom[lev]}, {grids[lev]}, - {dmap[lev]}, info ); #endif linop.setSigma(lev, *sigma[lev]); From 4ad711239f52886cd1933d3ae55a6b1868049429 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Tue, 30 Jul 2024 16:41:12 -0700 Subject: [PATCH 08/28] add example script --- .../PICMI_inputs_2d.py | 234 ++++++++++++++++++ Regression/WarpX-tests.ini | 16 ++ 2 files changed, 250 insertions(+) create mode 100644 Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py diff --git a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py new file mode 100644 index 00000000000..be5b7915c55 --- /dev/null +++ b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 +# +# --- Test script for the semi-implicit Poisson solver in which a plasma column +# --- expands inside a conducting cylinder. + +import argparse +import sys + +import numpy as np +from mpi4py import MPI as mpi + +from pywarpx import picmi + +constants = picmi.constants + +comm = mpi.COMM_WORLD + +simulation = picmi.Simulation( + warpx_serialize_initial_conditions=True, + verbose=0 +) + + +class PlasmaExpansionSimulation(object): + '''Simulation setup for an expanding plasma.''' + # Plasma parameters + m_ion = 250.0 # Ion mass (electron masses) + n_plasma = 1e18 # Plasma density (m^-3) + T_e = 2.0 # Electron temperature (eV) + T_i = 1.0 # Ion temperature (eV) + + scale_fac = 4.0 # scaling factor used to move between explicit and SI + + # Plasma size + plasma_radius = 20 # Plasma column radius (Debye lengths) + + # Spatial domain + NZ = 256 / scale_fac # Number of cells in x and z directions + + # Temporal domain (if not run as a CI test) + LT = 0.1 # Simulation temporal length in ion crossing times + + # Numerical parameters + NPPC = 750 # Seed number of particles per cell + DZ = 1.0 * scale_fac # Cell size (Debye lengths) + DT = 0.75 # Time step (electron CFL) + + # Solver parameter + C_SI = 20 # Semi-implicit factor + + def __init__(self, test, verbose): + """Get input parameters for the specific case desired.""" + self.test = test + self.verbose = verbose or self.test + + # calculate various plasma parameters based on the simulation input + self.get_plasma_quantities() + + self.dz = self.DZ * self.lambda_e + self.Lz = self.NZ * self.dz + + self.dt = self.DT * self.dz / self.v_te + + self.plasma_radius *= self.lambda_e + + self.total_steps = int(self.LT * self.Lz / (self.v_ti * self.dt)) + self.diag_steps = max(50, self.total_steps // 10) + + # print out plasma parameters + if comm.rank == 0: + print( + f"Initializing simulation with input parameters:\n" + f"\tT_e = {self.T_e:.1f} eV\n" + f"\tT_i = {self.T_i:.1f} eV\n" + f"\tn = {self.n_plasma:.1e} m^-3\n" + f"\tM/m = {self.m_ion:.0f}\n" + ) + print( + f"Plasma parameters:\n" + f"\tlambda_e = {self.lambda_e:.1e} m\n" + f"\tt_pe = {1.0/self.f_pe:.1e} s\n" + f"\tv_ti = {self.v_ti:.1e} m/s\n" + ) + print( + f"Numerical parameters:\n" + f"\tdz = {self.dz:.1e} m\n" + f"\tdt*f_pe = {self.dt*self.f_pe:.2f}\n" + f"\tdiag steps = {self.diag_steps:d}\n" + f"\ttotal steps = {self.total_steps:d}\n" + ) + self.setup_run() + + def get_plasma_quantities(self): + """Calculate various plasma parameters based on the simulation input.""" + # Ion mass (kg) + self.M = self.m_ion * constants.m_e + + # Electron plasma frequency (Hz) + self.f_pe = np.sqrt( + constants.q_e**2 * self.n_plasma / (constants.m_e * constants.ep0) + ) / (2.0 * np.pi) + + # Debye length (m) + self.lambda_e = np.sqrt( + constants.ep0 * self.T_e / (self.n_plasma * constants.q_e) + ) + + # Thermal velocities (m/s) from v_th = np.sqrt(kT / m) + self.v_ti = np.sqrt(self.T_i * constants.q_e / self.M) + self.v_te = np.sqrt(self.T_e * constants.q_e / constants.m_e) + + def setup_run(self): + """Setup simulation components.""" + + ####################################################################### + # Set geometry and boundary conditions # + ####################################################################### + + self.grid = picmi.Cartesian2DGrid( + number_of_cells=[self.NZ, self.NZ], + lower_bound=[-self.Lz/2.0, -self.Lz/2.0], + upper_bound=[self.Lz/2.0, self.Lz/2.0], + lower_boundary_conditions=['dirichlet']*2, + upper_boundary_conditions=['dirichlet']*2, + lower_boundary_conditions_particles=['absorbing']*2, + upper_boundary_conditions_particles=['absorbing']*2, + warpx_max_grid_size=self.NZ + ) + simulation.time_step_size = self.dt + simulation.max_steps = self.total_steps + simulation.current_deposition_algo = 'direct' + simulation.particle_shape = 1 + simulation.verbose = self.verbose + + ####################################################################### + # Insert probe as embedded boundary # + ####################################################################### + + embedded_boundary = picmi.EmbeddedBoundary( + implicit_function=f"(x**2+z**2-{self.Lz/2.-8*self.lambda_e}**2)", + potential=0.0, + ) + simulation.embedded_boundary = embedded_boundary + + ####################################################################### + # Field solver and external field # + ####################################################################### + + solver = picmi.ElectrostaticSolver( + grid=self.grid, method='Multigrid', + warpx_semi_implicit=True, warpx_semi_implicit_factor=self.C_SI, + warpx_self_fields_verbosity=0 + ) + simulation.solver = solver + + ####################################################################### + # Particle types setup # + ####################################################################### + + density_expr = f"if(x**2+z**2<{self.plasma_radius**2},{self.n_plasma},0)" + self.electrons = picmi.Species( + name='electrons', particle_type='electron', + initial_distribution=picmi.AnalyticDistribution( + density_expression=density_expr, + warpx_momentum_spread_expressions=[self.v_te]*3, + warpx_density_min=self.n_plasma/10.0, + ) + ) + simulation.add_species( + self.electrons, + layout=picmi.PseudoRandomLayout( + grid=self.grid, n_macroparticles_per_cell=self.NPPC + ) + ) + + self.ions = picmi.Species( + name='ions', charge='q_e', mass=self.M, + initial_distribution=picmi.AnalyticDistribution( + density_expression=density_expr, + warpx_momentum_spread_expressions=[self.v_ti]*3, + warpx_density_min=self.n_plasma/10.0, + ) + ) + simulation.add_species( + self.ions, + layout=picmi.PseudoRandomLayout( + grid=self.grid, n_macroparticles_per_cell=self.NPPC + ) + ) + + ####################################################################### + # Add diagnostics # + ####################################################################### + + field_diag = picmi.FieldDiagnostic( + name='field_diag', + grid=self.grid, + period=self.diag_steps, + write_dir='.', + data_list=[ + 'E', 'J', 'T_electrons', 'T_ions', 'phi', + 'rho_electrons', 'rho_ions' + ], + warpx_file_prefix='Python_si_electrostatic_plt', + warpx_format='openpmd', + warpx_openpmd_backend='h5' + ) + simulation.add_diagnostic(field_diag) + + ####################################################################### + # Initialize simulation # + ####################################################################### + + simulation.initialize_inputs() + simulation.initialize_warpx() + + +########################## +# parse input parameters +########################## + +parser = argparse.ArgumentParser() +parser.add_argument( + '-t', '--test', help='toggle whether this script is run as a short CI test', + action='store_true', +) +parser.add_argument( + '-v', '--verbose', help='Verbose output', action='store_true', +) +args, left = parser.parse_known_args() +sys.argv = sys.argv[:1]+left + +run = PlasmaExpansionSimulation(test=args.test, verbose=args.verbose) +simulation.step() diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index dcf53b5cc98..b85eaf057f4 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2735,6 +2735,22 @@ useOMP = 1 numthreads = 1 analysisRoutine = Examples/analysis_default_regression.py +[Python_si_electrostatic] +buildDir = . +inputFile = Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py +runtime_params = +customRunCmd = python3 PICMI_inputs_2d.py +dim = 2 +addToCompileString = USE_EB=TRUE USE_PYTHON_MAIN=TRUE +cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_EB=ON +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +analysisRoutine = Examples/analysis_default_regression.py + [Python_LaserAcceleration] buildDir = . inputFile = Examples/Physics_applications/laser_acceleration/PICMI_inputs_3d.py From 5e9b96faecfaee6d48588256e4fb386d5dad8c96 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Sun, 4 Aug 2024 14:39:18 -0700 Subject: [PATCH 09/28] fix clang-tidy issues --- .../PICMI_inputs_2d.py | 4 +- .../Python_si_electrostatic.json | 15 ++++++++ .../fields/SemiImplicitPoissonSolver.H | 37 +++---------------- 3 files changed, 21 insertions(+), 35 deletions(-) create mode 100644 Regression/Checksum/benchmarks_json/Python_si_electrostatic.json diff --git a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py index be5b7915c55..64f8afe2212 100644 --- a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py +++ b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py @@ -201,9 +201,7 @@ def setup_run(self): 'E', 'J', 'T_electrons', 'T_ions', 'phi', 'rho_electrons', 'rho_ions' ], - warpx_file_prefix='Python_si_electrostatic_plt', - warpx_format='openpmd', - warpx_openpmd_backend='h5' + warpx_file_prefix='Python_si_electrostatic_plt' ) simulation.add_diagnostic(field_diag) diff --git a/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json b/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json new file mode 100644 index 00000000000..c1e70d43629 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json @@ -0,0 +1,15 @@ +{ + "lev=0": { + "Ex": 1170672.5304594298, + "Ey": 0.0, + "Ez": 1181500.3223771458, + "T_electrons": 1699.5283920623333, + "T_ions": 854.404135786121, + "jx": 160011.99110381026, + "jy": 465704.0865225781, + "jz": 154829.447754826, + "phi": 523.2079772374389, + "rho_electrons": 12.0631083127128, + "rho_ions": 12.1004924341728 + } +} diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index ae74de88233..a0c4ecaae9b 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -56,13 +56,13 @@ namespace ablastr::fields { -/** Compute the potential `phi` by solving the Poisson equation +/** Compute the potential `phi` by solving the Poisson equation with a modifed dielectric function * * Uses `rho` as a source. This uses the AMReX solver. * * More specifically, this solves the equation * \f[ - * ... \phi = -\frac{r \rho}{\epsilon_0} + * ... \phi(r) = -\frac{r \rho(r)}{\epsilon(r)} * \f] * * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler @@ -102,7 +102,7 @@ computeSemiImplicitPhi (amrex::Vector const & rho, amrex::Vector const& geom, amrex::Vector const& dmap, amrex::Vector const& grids, - T_BoundaryHandler const boundary_handler, + T_BoundaryHandler const & boundary_handler, bool is_solver_igf_on_lev0, bool const do_single_precision_comms = false, std::optional > rel_ref_ratio = std::nullopt, @@ -120,6 +120,8 @@ computeSemiImplicitPhi (amrex::Vector const & rho, "rel_ref_ratio must be set if mesh-refinement is used"); rel_ref_ratio = amrex::Vector{{amrex::IntVect(AMREX_D_DECL(1, 1, 1))}}; } + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, + "FFT solver cannot be used with semi-implicit Poisson solve"); auto const finest_level = static_cast(rho.size() - 1); @@ -140,14 +142,7 @@ computeSemiImplicitPhi (amrex::Vector const & rho, ); } -#if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) - amrex::LPInfo info; -#else const amrex::LPInfo info; -#endif - - ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, - "FFT solver cannot be used with semi-implicit Poisson solve"); for (int lev=0; lev<=finest_level; lev++) { @@ -155,28 +150,6 @@ computeSemiImplicitPhi (amrex::Vector const & rho, using namespace ablastr::constant::SI; rho[lev]->mult(-1._rt/ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! -// #if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) -// // Determine whether to use semi-coarsening -// amrex::Array dx_scaled -// {AMREX_D_DECL(geom[lev].CellSize(0), geom[lev].CellSize(1), geom[lev].CellSize(2))}; -// int max_semicoarsening_level = 0; -// int semicoarsening_direction = -1; -// const auto min_dir = static_cast(std::distance(dx_scaled.begin(), -// std::min_element(dx_scaled.begin(),dx_scaled.end()))); -// const auto max_dir = static_cast(std::distance(dx_scaled.begin(), -// std::max_element(dx_scaled.begin(),dx_scaled.end()))); -// if (dx_scaled[max_dir] > dx_scaled[min_dir]) { -// semicoarsening_direction = max_dir; -// max_semicoarsening_level = static_cast -// (std::log2(dx_scaled[max_dir]/dx_scaled[min_dir])); -// } -// if (max_semicoarsening_level > 0) { -// info.setSemicoarsening(true); -// info.setMaxSemicoarseningLevel(max_semicoarsening_level); -// info.setSemicoarseningDirection(semicoarsening_direction); -// } -// #endif - // Use the EB enabled solver. amrex::MLEBNodeFDLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info #if defined(AMREX_USE_EB) From 21a2a707f022ef37dde0a93a13d54812f4d91684 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 5 Aug 2024 09:14:21 -0700 Subject: [PATCH 10/28] more CI fixes --- .../Python_si_electrostatic.json | 20 +++++++++---------- .../SemiImplicitSolver.cpp | 18 ++++++++--------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json b/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json index c1e70d43629..7edc638974d 100644 --- a/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json +++ b/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json @@ -1,15 +1,15 @@ { "lev=0": { - "Ex": 1170672.5304594298, + "Ex": 1166229.3217510823, "Ey": 0.0, - "Ez": 1181500.3223771458, - "T_electrons": 1699.5283920623333, - "T_ions": 854.404135786121, - "jx": 160011.99110381026, - "jy": 465704.0865225781, - "jz": 154829.447754826, - "phi": 523.2079772374389, - "rho_electrons": 12.0631083127128, - "rho_ions": 12.1004924341728 + "Ez": 1236947.9176095286, + "T_electrons": 1702.6191883360068, + "T_ions": 858.7295550843751, + "jx": 154473.7772622879, + "jy": 470386.6132005394, + "jz": 153606.65643944603, + "phi": 500.8578828191274, + "rho_electrons": 12.0613993243032, + "rho_ions": 12.083188926525597 } } diff --git a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp index 3256341e801..eb59b8a7ba0 100644 --- a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp @@ -62,21 +62,21 @@ void SemiImplicitSolver::ComputeSigma () { int const lev = 0; sigma[lev]->setVal(1.0_rt); - // GetChargeDensity returns a nodal multifab - amrex::GpuArray nodal = {1, 1, 1}; // sigma is a cell-centered array - amrex::GpuArray const& cell_centered = {0, 0, 0}; + amrex::GpuArray const cell_centered = {0, 0, 0}; // The "coarsening is just 1 i.e. no coarsening" - amrex::GpuArray const& coarsen = {1, 1, 1}; + amrex::GpuArray const coarsen = {1, 1, 1}; + // GetChargeDensity returns a nodal multifab // Below we set all the unused dimensions to have cell-centered values for // rho since these values will be interpolated onto a cell-centered grid // - if this is not done the Interp function returns nonsense values. -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_1D_Z) - nodal[2] = 0; -#endif -#if defined(WARPX_DIM_1D_Z) - nodal[1] = 0; +#if defined(WARPX_DIM_3D) + amrex::GpuArray const nodal = {1, 1, 1}; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::GpuArray const nodal = {1, 1, 0}; +#elif defined(WARPX_DIM_1D_Z) + amrex::GpuArray const nodal = {1, 0, 0}; #endif auto& warpx = WarpX::GetInstance(); From 56b40df24c071ec7302dd86cc75feda6fc4658ac Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 5 Aug 2024 14:30:39 -0700 Subject: [PATCH 11/28] change definition of C_SI to match Aleph --- .../PICMI_inputs_2d.py | 2 +- .../ElectrostaticSolver/SemiImplicitSolver.H | 2 +- .../SemiImplicitSolver.cpp | 26 ++++++++++++------- 3 files changed, 18 insertions(+), 12 deletions(-) diff --git a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py index 64f8afe2212..46ebcc2201b 100644 --- a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py +++ b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py @@ -46,7 +46,7 @@ class PlasmaExpansionSimulation(object): DT = 0.75 # Time step (electron CFL) # Solver parameter - C_SI = 20 # Semi-implicit factor + C_SI = 2 # Semi-implicit factor - marginal stability threshold = 1 def __init__(self, test, verbose): """Get input parameters for the specific case desired.""" diff --git a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H index 45b96b7377d..c574ef586fb 100644 --- a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.H @@ -45,7 +45,7 @@ public: // Declare multifabs specifically needed for the semi-implicit solver amrex::Vector > sigma; - amrex::Real m_C_SI = 6.0; + amrex::Real m_C_SI = 4.0; }; #endif // WARPX_SEMIIMPLICITSOLVER_H_ diff --git a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp index eb59b8a7ba0..5cd8908dd1e 100644 --- a/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver/SemiImplicitSolver.cpp @@ -82,20 +82,24 @@ void SemiImplicitSolver::ComputeSigma () { auto& warpx = WarpX::GetInstance(); auto& mypc = warpx.GetPartContainer(); - auto m_mult_factor = ( - m_C_SI * warpx.getdt(lev) * warpx.getdt(lev) - / (16._rt * MathConst::pi * MathConst::pi * PhysConst::ep0) + // The semi-implicit dielectric function is given by + // \varepsilon_{SI} = \varepsilon * (1 + \sum_{i in species} C_{SI}*(w_pi * dt)^2/4) + // Note the use of the plasma frequency in rad/s (not Hz) and the factor of 1/4, + // these choices make it so that C_SI = 1 is the marginal stability threshold. + auto mult_factor = ( + m_C_SI * warpx.getdt(lev) * warpx.getdt(lev) / (4._rt * PhysConst::ep0) ); // Loop over each species to calculate the Poisson equation dressing for (auto const& pc : mypc) { + // grab the charge density for this species auto rho = pc->GetChargeDensity(lev, false); - // Handle the parallel transfers of guard cells and - // apply the filtering if requested - might not be needed... + + // Handle the parallel transfer of guard cells and apply filtering warpx.ApplyFilterandSumBoundaryRho(lev, lev, *rho, 0, rho->nComp()); - // multiply charge density by C_SI * dt**2 /4.0 q/(m \varepsilon_0) - rho->mult(m_mult_factor * pc->getCharge() / pc->getMass()); + // get multiplication factor for this species + auto const mult_factor_pc = mult_factor * pc->getCharge() / pc->getMass(); // update sigma #ifdef AMREX_USE_OMP @@ -107,11 +111,13 @@ void SemiImplicitSolver::ComputeSigma () { // Loop over the cells and update the sigma field amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Interpolate rho to cell-centered multifab to add to sigma. - auto const rho_interp = ablastr::coarsen::sample::Interp( + // Interpolate rho to cell-centered value + auto const rho_cc = ablastr::coarsen::sample::Interp( rho_arr, nodal, cell_centered, coarsen, i, j, k, 0 ); - sigma_arr(i, j, k, 0) += rho_interp; + // add species term to sigma: + // C_SI * w_p^2 * dt^2 / 4 = C_SI / 4 * q*rho/(m*eps0) * dt^2 + sigma_arr(i, j, k, 0) += mult_factor_pc * rho_cc; }); } From 6f29954d1098d23f9ad7c852acdfe60cd0bf22cf Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 5 Aug 2024 16:45:41 -0700 Subject: [PATCH 12/28] reset benchmark values --- .../Python_si_electrostatic.json | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json b/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json index 7edc638974d..c13b6d911fc 100644 --- a/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json +++ b/Regression/Checksum/benchmarks_json/Python_si_electrostatic.json @@ -1,15 +1,15 @@ { "lev=0": { - "Ex": 1166229.3217510823, + "Ex": 1428856.0323152007, "Ey": 0.0, - "Ez": 1236947.9176095286, - "T_electrons": 1702.6191883360068, - "T_ions": 858.7295550843751, - "jx": 154473.7772622879, - "jy": 470386.6132005394, - "jz": 153606.65643944603, - "phi": 500.8578828191274, - "rho_electrons": 12.0613993243032, - "rho_ions": 12.083188926525597 + "Ez": 1414591.6479792167, + "T_electrons": 1652.0361205183694, + "T_ions": 851.0244739005777, + "jx": 151920.26162133453, + "jy": 453186.62923648424, + "jz": 159070.6572202958, + "phi": 705.4775524885782, + "rho_electrons": 12.100065187070399, + "rho_ions": 12.149198603846399 } } From 818883935f93333160b88e42481fe3bcf8ee0f5c Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Tue, 6 Aug 2024 10:25:14 -0700 Subject: [PATCH 13/28] update documentation --- Docs/source/usage/parameters.rst | 10 ++++++++++ Python/pywarpx/picmi.py | 10 ++++++++++ Source/FieldSolver/ElectrostaticSolver.cpp | 1 + 3 files changed, 21 insertions(+) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 4a8c4be0b5f..9ef6fdf6d75 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -228,6 +228,16 @@ Overall simulation parameters \boldsymbol{\nabla}^2 \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi \\ \boldsymbol{\nabla}^2 \boldsymbol{A} = - \mu_0 \boldsymbol{j} \qquad \boldsymbol{B} = \boldsymbol{\nabla}\times\boldsymbol{A} + * ``labframe-semi-implicit``: Poisson's equation is solved with a modified dielectric function + to create a semi-implicit scheme in which grid heating effects due to unresolved + plasma modes are suppressed. If this option is used the additional parameter + ``warpx.semi_implicit_factor`` can also be specified to set the value of :math:`C_{SI}`. + The method is marginally stable for :math:`C_{SI} = 1`. Specifically, the code solves: + + .. math:: + + \boldsymbol{\nabla}\cdot\left(1+\frac{C_{SI}}{4}\sum_{s \in \text{species}}(\omega_{ps}\Delta t)^2 \right)\boldsymbol{\nabla} \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi + * ``relativistic``: Poisson's equation is solved **for each species** in their respective rest frame. The corresponding field is mapped back to the simulation frame and will produce both E and B diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index cca6de13e1d..499f881483b 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1545,6 +1545,16 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): warpx_self_fields_verbosity: integer, default=2 Level of verbosity for the lab frame solver + + warpx_magnetostatic: bool, default=False + Whether to use the magnetostatic solver + + warpx_semi_implicit: bool, default=False + Whether to use the semi-implicit Poisson solver + + warpx_semi_implicit_factor: float, default=4 + If the semi-implicit Poisson solver is used, this sets the value + of C_SI (the method is marginally stable at C_SI = 1) """ def init(self, kw): self.relativistic = kw.pop('warpx_relativistic', False) diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index f2bef2b37b8..31bb21b8791 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -450,6 +450,7 @@ WarpX::computePhi (const amrex::Vector >& rho, this->geom, this->dmap, this->grids, + WarpX::grid_type, this->m_poisson_boundary_handler, is_solver_igf_on_lev0, WarpX::do_single_precision_comms, From 75d4226b314e960f1f28e84168d6d3559a886e0f Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Thu, 8 Aug 2024 14:53:22 -0700 Subject: [PATCH 14/28] also reset fields for SI now that #5104 is merged --- Source/Evolve/WarpXEvolve.cpp | 2 +- Source/ablastr/fields/SemiImplicitPoissonSolver.H | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 24a165853b3..32f7a493916 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -238,7 +238,7 @@ WarpX::Evolve (int numsteps) // loop (i.e. immediately after a `Redistribute` and before particle // positions are next pushed) so that the particles do not deposit out of bounds // and so that the fields are at the correct time in the output. - bool const reset_fields = (electrostatic_solver_id != ElectrostaticSolverAlgo::LabFrameSemiImplicit); + bool const reset_fields = true; ComputeSpaceChargeField( reset_fields ); if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { // Call Magnetostatic Solver to solve for the vector potential A and compute the diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index a0c4ecaae9b..1d7b383b485 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -62,7 +62,7 @@ namespace ablastr::fields { * * More specifically, this solves the equation * \f[ - * ... \phi(r) = -\frac{r \rho(r)}{\epsilon(r)} + * \nabla \cdot \sigma \nabla \phi = - \rho/\epsilon_0 * \f] * * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler From 6e702a568e8c6966e44a15e7b4f9a3a3d07b3b62 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Wed, 14 Aug 2024 15:06:19 -0700 Subject: [PATCH 15/28] use `amrex::MLNodeLaplacian` if EB support is OFF --- Source/FieldSolver/ElectrostaticSolver.cpp | 31 ++----------------- .../fields/SemiImplicitPoissonSolver.H | 4 +++ 2 files changed, 6 insertions(+), 29 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 31bb21b8791..a12ad35a7dd 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -338,9 +338,9 @@ WarpX::computePhi (const amrex::Vector >& rho, sorted_phi.emplace_back(phi[lev].get()); } +#if defined(AMREX_USE_EB) std::optional post_phi_calculation; -#if defined(AMREX_USE_EB) // EB: use AMReX to directly calculate the electric field since with EB's the // simple finite difference scheme in WarpX::computeE sometimes fails if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || @@ -383,34 +383,7 @@ WarpX::computePhi (const amrex::Vector >& rho, } eb_farray_box_factory = factories; #else - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) - { - // TODO: maybe make this a helper function or pass Efield_fp directly - amrex::Vector< - amrex::Array - > e_field; - for (int lev = 0; lev <= finest_level; ++lev) { - e_field.push_back( -# if defined(WARPX_DIM_1D_Z) - amrex::Array{ - getFieldPointer(FieldType::Efield_fp, lev, 2) - } -# elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Array{ - getFieldPointer(FieldType::Efield_fp, lev, 0), - getFieldPointer(FieldType::Efield_fp, lev, 2) - } -# elif defined(WARPX_DIM_3D) - amrex::Array{ - getFieldPointer(FieldType::Efield_fp, lev, 0), - getFieldPointer(FieldType::Efield_fp, lev, 1), - getFieldPointer(FieldType::Efield_fp, lev, 2) - } -# endif - ); - } - post_phi_calculation = ElectrostaticSolver::EBCalcEfromPhiPerLevel(e_field); - } + const std::optional post_phi_calculation; const std::optional > eb_farray_box_factory; #endif diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index 1d7b383b485..0c903c33f2b 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -150,12 +150,16 @@ computeSemiImplicitPhi (amrex::Vector const & rho, using namespace ablastr::constant::SI; rho[lev]->mult(-1._rt/ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! +#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) // Use the EB enabled solver. amrex::MLEBNodeFDLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info #if defined(AMREX_USE_EB) , {eb_farray_box_factory.value()[lev]} #endif ); +#else + amrex::MLNodeLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info); +#endif #if defined(AMREX_USE_EB) // if the EB potential only depends on time, the potential can be passed From fd0bfbc8c4491b1201a3b6a97b39c0efdc5c053c Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Sun, 18 Aug 2024 11:35:40 -0700 Subject: [PATCH 16/28] update license statement Signed-off-by: roelof-groenewald --- Source/ablastr/fields/SemiImplicitPoissonSolver.H | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index 0c903c33f2b..b8a0596ef6a 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -1,9 +1,14 @@ -/* Copyright 2024 Roelof Groenewald (TAE Technologies) +/* Copyright 2024 The WarpX Community * * This file is part of WarpX. * + * Authors: Roelof Groenewald (TAE Technologies) + * * License: BSD-3-Clause-LBNL */ +/* + * This file was copied and edited from PoissonSolver.H in the same directory. + */ #ifndef ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H #define ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H From e0b48384d2daa0ce241546c477728294fa78f11c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 21 Aug 2024 01:03:09 +0000 Subject: [PATCH 17/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../PICMI_inputs_2d.py | 107 ++++++++++-------- 1 file changed, 60 insertions(+), 47 deletions(-) diff --git a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py index 46ebcc2201b..effff279f12 100644 --- a/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py +++ b/Examples/Tests/semi_implicit_electrostatic/PICMI_inputs_2d.py @@ -15,38 +15,36 @@ comm = mpi.COMM_WORLD -simulation = picmi.Simulation( - warpx_serialize_initial_conditions=True, - verbose=0 -) +simulation = picmi.Simulation(warpx_serialize_initial_conditions=True, verbose=0) class PlasmaExpansionSimulation(object): - '''Simulation setup for an expanding plasma.''' + """Simulation setup for an expanding plasma.""" + # Plasma parameters - m_ion = 250.0 # Ion mass (electron masses) - n_plasma = 1e18 # Plasma density (m^-3) - T_e = 2.0 # Electron temperature (eV) - T_i = 1.0 # Ion temperature (eV) + m_ion = 250.0 # Ion mass (electron masses) + n_plasma = 1e18 # Plasma density (m^-3) + T_e = 2.0 # Electron temperature (eV) + T_i = 1.0 # Ion temperature (eV) - scale_fac = 4.0 # scaling factor used to move between explicit and SI + scale_fac = 4.0 # scaling factor used to move between explicit and SI # Plasma size - plasma_radius = 20 # Plasma column radius (Debye lengths) + plasma_radius = 20 # Plasma column radius (Debye lengths) # Spatial domain - NZ = 256 / scale_fac # Number of cells in x and z directions + NZ = 256 / scale_fac # Number of cells in x and z directions # Temporal domain (if not run as a CI test) - LT = 0.1 # Simulation temporal length in ion crossing times + LT = 0.1 # Simulation temporal length in ion crossing times # Numerical parameters - NPPC = 750 # Seed number of particles per cell - DZ = 1.0 * scale_fac # Cell size (Debye lengths) - DT = 0.75 # Time step (electron CFL) + NPPC = 750 # Seed number of particles per cell + DZ = 1.0 * scale_fac # Cell size (Debye lengths) + DT = 0.75 # Time step (electron CFL) # Solver parameter - C_SI = 2 # Semi-implicit factor - marginal stability threshold = 1 + C_SI = 2 # Semi-implicit factor - marginal stability threshold = 1 def __init__(self, test, verbose): """Get input parameters for the specific case desired.""" @@ -118,17 +116,17 @@ def setup_run(self): self.grid = picmi.Cartesian2DGrid( number_of_cells=[self.NZ, self.NZ], - lower_bound=[-self.Lz/2.0, -self.Lz/2.0], - upper_bound=[self.Lz/2.0, self.Lz/2.0], - lower_boundary_conditions=['dirichlet']*2, - upper_boundary_conditions=['dirichlet']*2, - lower_boundary_conditions_particles=['absorbing']*2, - upper_boundary_conditions_particles=['absorbing']*2, - warpx_max_grid_size=self.NZ + lower_bound=[-self.Lz / 2.0, -self.Lz / 2.0], + upper_bound=[self.Lz / 2.0, self.Lz / 2.0], + lower_boundary_conditions=["dirichlet"] * 2, + upper_boundary_conditions=["dirichlet"] * 2, + lower_boundary_conditions_particles=["absorbing"] * 2, + upper_boundary_conditions_particles=["absorbing"] * 2, + warpx_max_grid_size=self.NZ, ) simulation.time_step_size = self.dt simulation.max_steps = self.total_steps - simulation.current_deposition_algo = 'direct' + simulation.current_deposition_algo = "direct" simulation.particle_shape = 1 simulation.verbose = self.verbose @@ -147,9 +145,11 @@ def setup_run(self): ####################################################################### solver = picmi.ElectrostaticSolver( - grid=self.grid, method='Multigrid', - warpx_semi_implicit=True, warpx_semi_implicit_factor=self.C_SI, - warpx_self_fields_verbosity=0 + grid=self.grid, + method="Multigrid", + warpx_semi_implicit=True, + warpx_semi_implicit_factor=self.C_SI, + warpx_self_fields_verbosity=0, ) simulation.solver = solver @@ -159,33 +159,36 @@ def setup_run(self): density_expr = f"if(x**2+z**2<{self.plasma_radius**2},{self.n_plasma},0)" self.electrons = picmi.Species( - name='electrons', particle_type='electron', + name="electrons", + particle_type="electron", initial_distribution=picmi.AnalyticDistribution( density_expression=density_expr, - warpx_momentum_spread_expressions=[self.v_te]*3, - warpx_density_min=self.n_plasma/10.0, - ) + warpx_momentum_spread_expressions=[self.v_te] * 3, + warpx_density_min=self.n_plasma / 10.0, + ), ) simulation.add_species( self.electrons, layout=picmi.PseudoRandomLayout( grid=self.grid, n_macroparticles_per_cell=self.NPPC - ) + ), ) self.ions = picmi.Species( - name='ions', charge='q_e', mass=self.M, + name="ions", + charge="q_e", + mass=self.M, initial_distribution=picmi.AnalyticDistribution( density_expression=density_expr, - warpx_momentum_spread_expressions=[self.v_ti]*3, - warpx_density_min=self.n_plasma/10.0, - ) + warpx_momentum_spread_expressions=[self.v_ti] * 3, + warpx_density_min=self.n_plasma / 10.0, + ), ) simulation.add_species( self.ions, layout=picmi.PseudoRandomLayout( grid=self.grid, n_macroparticles_per_cell=self.NPPC - ) + ), ) ####################################################################### @@ -193,15 +196,20 @@ def setup_run(self): ####################################################################### field_diag = picmi.FieldDiagnostic( - name='field_diag', + name="field_diag", grid=self.grid, period=self.diag_steps, - write_dir='.', + write_dir=".", data_list=[ - 'E', 'J', 'T_electrons', 'T_ions', 'phi', - 'rho_electrons', 'rho_ions' + "E", + "J", + "T_electrons", + "T_ions", + "phi", + "rho_electrons", + "rho_ions", ], - warpx_file_prefix='Python_si_electrostatic_plt' + warpx_file_prefix="Python_si_electrostatic_plt", ) simulation.add_diagnostic(field_diag) @@ -219,14 +227,19 @@ def setup_run(self): parser = argparse.ArgumentParser() parser.add_argument( - '-t', '--test', help='toggle whether this script is run as a short CI test', - action='store_true', + "-t", + "--test", + help="toggle whether this script is run as a short CI test", + action="store_true", ) parser.add_argument( - '-v', '--verbose', help='Verbose output', action='store_true', + "-v", + "--verbose", + help="Verbose output", + action="store_true", ) args, left = parser.parse_known_args() -sys.argv = sys.argv[:1]+left +sys.argv = sys.argv[:1] + left run = PlasmaExpansionSimulation(test=args.test, verbose=args.verbose) simulation.step() From cb6a5cad3cdde5a377ec7af2cdd32f2137c12996 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Sat, 28 Sep 2024 17:44:12 -0700 Subject: [PATCH 18/28] minor cleanup Signed-off-by: roelof-groenewald --- .../FieldSolver/ElectrostaticSolvers/SemiImplicitES.H | 11 +---------- Source/Python/WarpX.cpp | 1 - Source/ablastr/fields/SemiImplicitPoissonSolver.H | 3 ++- 3 files changed, 3 insertions(+), 12 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H index acf6b433d49..143209c2b17 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H +++ b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H @@ -15,15 +15,6 @@ #include #include -// #include "Particles/MultiParticleContainer.H" -// #include "Utils/WarpXConst.H" -// #include "Utils/WarpXProfilerWrapper.H" -// #include "WarpX.H" - -// #include -// #include -// #include - class SemiImplicitES final : public ElectrostaticSolver { public: @@ -55,7 +46,7 @@ public: * \f] * \param[out] phi The potential to be computed by this function * \param[in] rho The total charge density - * \param[in] sigma Represents the modified + * \param[in] sigma Represents the modified dielectric * \param[in] required_precision The relative convergence threshold for the MLMG solver * \param[in] absolute_tolerance The absolute convergence threshold for the MLMG solver * \param[in] max_iters The maximum number of iterations allowed for the MLMG solver diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp index 77dfb634d90..0aab95f78f8 100644 --- a/Source/Python/WarpX.cpp +++ b/Source/Python/WarpX.cpp @@ -11,7 +11,6 @@ #include #include #include -#include "FieldSolver/ElectrostaticSolvers/SemiImplicitES.H" #include #include #include diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index 9dd00fd0d7e..a63a715ecf6 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -174,8 +174,9 @@ computeSemiImplicitPhi ( if (boundary_handler.phi_EB_only_t) { linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value())); } - else + else { linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); + } #endif linop.setSigma(lev, sigma); From b23d3c61dc24617270ba6c3ff04e734da08af863 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 30 Sep 2024 11:34:19 -0700 Subject: [PATCH 19/28] overload compute phi function in semi-implicit Poisson solver Signed-off-by: roelof-groenewald --- .../ElectrostaticSolvers/SemiImplicitES.H | 6 +++++- .../ElectrostaticSolvers/SemiImplicitES.cpp | 16 ++++++++++++---- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H index 143209c2b17..48105d4de9d 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H +++ b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H @@ -35,7 +35,7 @@ public: MultiFluidContainer* mfl, int max_level) override; - void ComputeSigma ( amrex::MultiFab& sigma ); + void ComputeSigma ( amrex::MultiFab& sigma ) const; /** * Compute the potential `phi` by solving the semi-implicit Poisson equation @@ -52,6 +52,10 @@ public: * \param[in] max_iters The maximum number of iterations allowed for the MLMG solver * \param[in] verbosity The verbosity setting for the MLMG solver */ + void computePhi ( + ablastr::fields::MultiLevelScalarField const& rho, + ablastr::fields::MultiLevelScalarField const& phi + ) const; void computePhi ( ablastr::fields::MultiLevelScalarField const& rho, ablastr::fields::MultiLevelScalarField const& phi, diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp index f06d0d45d70..cf159abfe85 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp @@ -61,20 +61,28 @@ void SemiImplicitES::ComputeSpaceChargeField ( // set the boundary potentials appropriately setPhiBC(phi_fp, warpx.gett_new(0)); + // perform phi calculation + computePhi(rho_fp, phi_fp); +} + +void SemiImplicitES::computePhi ( + ablastr::fields::MultiLevelScalarField const& rho, + ablastr::fields::MultiLevelScalarField const& phi) const +{ // Calculate the mass enhancement factor: // The "sigma" multifab stores the dressing of the Poisson equation. It // is a cell-centered multifab. - auto const& ba = convert(rho_fp[0]->boxArray(), IntVect(AMREX_D_DECL(0,0,0))); - MultiFab sigma(ba, rho_fp[0]->DistributionMap(), 1, rho_fp[0]->nGrowVect()); + auto const& ba = convert(rho[0]->boxArray(), IntVect(AMREX_D_DECL(0,0,0))); + MultiFab sigma(ba, rho[0]->DistributionMap(), 1, rho[0]->nGrowVect()); ComputeSigma(sigma); // Use the AMREX MLMG solver - computePhi(rho_fp, phi_fp, sigma, self_fields_required_precision, + computePhi(rho, phi, sigma, self_fields_required_precision, self_fields_absolute_tolerance, self_fields_max_iters, self_fields_verbosity); } -void SemiImplicitES::ComputeSigma (MultiFab& sigma) +void SemiImplicitES::ComputeSigma (MultiFab& sigma) const { // Reset sigma to 1 sigma.setVal(1.0_rt); From 3063ad8fb66cd9f473ccad991b1826bf8c115818 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Tue, 22 Oct 2024 21:48:00 -0700 Subject: [PATCH 20/28] update documentation based on review comments --- Docs/source/refs.bib | 12 ++++++++++++ Docs/source/usage/parameters.rst | 13 +++++++++---- Source/Diagnostics/Diagnostics.cpp | 2 +- 3 files changed, 22 insertions(+), 5 deletions(-) diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib index 130e0ce4da7..d9cdfc7b0e9 100644 --- a/Docs/source/refs.bib +++ b/Docs/source/refs.bib @@ -444,3 +444,15 @@ @article{Vranic2015 issn = {0010-4655}, doi = {https://doi.org/10.1016/j.cpc.2015.01.020}, } + +@article{Barnes2021, + title = {Improved C1 shape functions for simplex meshes}, + author = {D.C. Barnes}, + journal = {Journal of Computational Physics}, + volume = {424}, + pages = {109852}, + year = {2021}, + issn = {0021-9991}, + doi = {https://doi.org/10.1016/j.jcp.2020.109852}, + url = {https://www.sciencedirect.com/science/article/pii/S0021999120306264}, +} diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 86ae66daf77..a9fe038822d 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -229,10 +229,15 @@ Overall simulation parameters \boldsymbol{\nabla}^2 \boldsymbol{A} = - \mu_0 \boldsymbol{j} \qquad \boldsymbol{B} = \boldsymbol{\nabla}\times\boldsymbol{A} * ``labframe-semi-implicit``: Poisson's equation is solved with a modified dielectric function - to create a semi-implicit scheme in which grid heating effects due to unresolved - plasma modes are suppressed. If this option is used the additional parameter - ``warpx.semi_implicit_factor`` can also be specified to set the value of :math:`C_{SI}`. - The method is marginally stable for :math:`C_{SI} = 1`. Specifically, the code solves: + to create a semi-implicit scheme which is robust to the numerical instability seen + in explicit electrostatic PIC when :math:`\Delta t \omega_{pe} > 2`. + If this option is used the additional parameter ``warpx.semi_implicit_factor`` can also be + specified to set the value of :math:`C_{SI}` (default 4). The method is stable for :math:`C_{SI} \geq 1` + regardless of :math:`\Delta t`, however, the larger :math:`C_{SI}` is set, the lower the numerical plasma + frequency will be and therefore care must be taken to not set :math:`C_{SI}` so high that the plasma mode + hybridizes with other modes of interest. + Details of the method can be found in Appendix A of :cite:t:`param-Barnes2021`. + In short, the code solves: .. math:: diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index c88b4f55fee..939a058214b 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -79,7 +79,7 @@ Diagnostics::BaseReadParameters () warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrame || warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameSemiImplicit, - "plot phi only works if do_electrostatic = labframe or do_electrostatic = labframe-electromagnetostatic"); + "plot phi only works if do_electrostatic = labframe, do_electrostatic = labframe-electromagnetostatic or do_electrostatic = labframe-semi-implicit"); } // Sanity check if user requests to plot A From 7f5f3666923214240bf0c98bedcf6dfc43054a13 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Wed, 23 Oct 2024 10:03:58 -0700 Subject: [PATCH 21/28] update `SemiImplicitPoissonSolver.H` to not rely on compile time settings for EB support --- .../ElectrostaticSolvers/SemiImplicitES.cpp | 2 +- Source/ablastr/fields/PoissonSolver.H | 4 +- .../fields/SemiImplicitPoissonSolver.H | 185 ++++++++++-------- 3 files changed, 105 insertions(+), 86 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp index cf159abfe85..77023f82fec 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp @@ -240,12 +240,12 @@ void SemiImplicitES::computePhi ( warpx.DistributionMap(), warpx.boxArray(), WarpX::grid_type, - *m_poisson_boundary_handler, false, EB::enabled(), WarpX::do_single_precision_comms, warpx.refRatio(), post_phi_calculation, + *m_poisson_boundary_handler, warpx.gett_new(0), eb_farray_box_factory ); diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H index d7eeecead1b..aa9288fe950 100644 --- a/Source/ablastr/fields/PoissonSolver.H +++ b/Source/ablastr/fields/PoissonSolver.H @@ -178,12 +178,12 @@ inline void interpolatePhiBetweenLevels ( * \param[in] dmap the distribution mapping per level (e.g., from AmrMesh) * \param[in] grids the grids per level (e.g., from AmrMesh) * \param[in] grid_type Integer that corresponds to the type of grid used in the simulation (collocated, staggered, hybrid) - * \param[in] boundary_handler a handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler * \param[in] is_solver_igf_on_lev0 boolean to select the Poisson solver: 1 for FFT on level 0 & Multigrid on other levels, 0 for Multigrid on all levels * \param[in] eb_enabled solve with embedded boundaries * \param[in] do_single_precision_comms perform communications in single precision * \param[in] rel_ref_ratio mesh refinement ratio between levels (default: 1) * \param[in] post_phi_calculation perform a calculation per level directly after phi was calculated; required for embedded boundaries (default: none) + * \param[in] boundary_handler a handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler * \param[in] current_time the current time; required for embedded boundaries (default: none) * \param[in] eb_farray_box_factory a factory for field data, @see amrex::EBFArrayBoxFactory; required for embedded boundaries (default: none) */ @@ -210,7 +210,7 @@ computePhi ( bool do_single_precision_comms = false, std::optional > rel_ref_ratio = std::nullopt, [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation = std::nullopt, - [[maybe_unused]] T_BoundaryHandler const boundary_handler = std::nullopt, // only used for EB + [[maybe_unused]] T_BoundaryHandler const boundary_handler = std::nullopt, [[maybe_unused]] std::optional current_time = std::nullopt, // only used for EB [[maybe_unused]] std::optional > eb_farray_box_factory = std::nullopt // only used for EB ) diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index a63a715ecf6..bff1a206087 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -19,6 +19,7 @@ #include #include #include +#include "PoissonSolver.H" #if defined(WARPX_USE_FFT) && defined(WARPX_DIM_3D) #include @@ -70,8 +71,8 @@ namespace ablastr::fields { * \nabla \cdot \sigma \nabla \phi = - \rho/\epsilon_0 * \f] * - * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler * \tparam T_PostPhiCalculationFunctor a calculation per level directly after phi was calculated + * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler * \tparam T_FArrayBoxFactory usually nothing or an amrex::EBFArrayBoxFactory (EB ONLY) * \param[in] rho The charge density a given species * \param[out] phi The potential to be computed by this function @@ -83,17 +84,17 @@ namespace ablastr::fields { * \param[in] geom the geometry per level (e.g., from AmrMesh) * \param[in] dmap the distribution mapping per level (e.g., from AmrMesh) * \param[in] grids the grids per level (e.g., from AmrMesh) - * \param[in] boundary_handler a handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler * \param[in] is_solver_igf_on_lev0 boolean to select the Poisson solver: 1 for FFT on level 0 & Multigrid on other levels, 0 for Multigrid on all levels * \param[in] do_single_precision_comms perform communications in single precision * \param[in] rel_ref_ratio mesh refinement ratio between levels (default: 1) * \param[in] post_phi_calculation perform a calculation per level directly after phi was calculated; required for embedded boundaries (default: none) + * \param[in] boundary_handler a handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler * \param[in] current_time the current time; required for embedded boundaries (default: none) * \param[in] eb_farray_box_factory a factory for field data, @see amrex::EBFArrayBoxFactory; required for embedded boundaries (default: none) */ template< - typename T_BoundaryHandler, typename T_PostPhiCalculationFunctor = std::nullopt_t, + typename T_BoundaryHandler = std::nullopt_t, typename T_FArrayBoxFactory = void > void @@ -109,12 +110,12 @@ computeSemiImplicitPhi ( amrex::Vector const& dmap, amrex::Vector const& grids, [[maybe_unused]] utils::enums::GridType grid_type, - T_BoundaryHandler const boundary_handler, bool is_solver_igf_on_lev0, - [[maybe_unused]] bool eb_enabled = false, + bool eb_enabled = false, bool do_single_precision_comms = false, std::optional > rel_ref_ratio = std::nullopt, [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation = std::nullopt, + [[maybe_unused]] T_BoundaryHandler const boundary_handler = std::nullopt, [[maybe_unused]] std::optional current_time = std::nullopt, // only used for EB [[maybe_unused]] std::optional > eb_farray_box_factory = std::nullopt // only used for EB ) { @@ -127,27 +128,29 @@ computeSemiImplicitPhi ( "rel_ref_ratio must be set if mesh-refinement is used"); rel_ref_ratio = amrex::Vector{{amrex::IntVect(AMREX_D_DECL(1, 1, 1))}}; } + +#if !defined(AMREX_USE_EB) + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(!eb_enabled, + "Embedded boundary solve requested but not compiled in"); +#endif + if (eb_enabled && std::is_same_v) { + throw std::runtime_error("EB requested by eb_farray_box_factory not provided!"); + } + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, "FFT solver cannot be used with semi-implicit Poisson solve"); +#ifdef WARPX_DIM_RZ + constexpr bool is_rz = true; +#else + constexpr bool is_rz = false; +#endif + auto const finest_level = static_cast(rho.size() - 1); // determine if rho is zero everywhere - amrex::Real max_norm_b = 0.0; - for (int lev=0; lev<=finest_level; lev++) { - max_norm_b = amrex::max(max_norm_b, rho[lev]->norm0()); - } - amrex::ParallelDescriptor::ReduceRealMax(max_norm_b); - - const bool always_use_bnorm = (max_norm_b > 0); - if (!always_use_bnorm) { - if (absolute_tolerance == 0.0) { absolute_tolerance = amrex::Real(1e-6); } - ablastr::warn_manager::WMRecordWarning( - "ElectrostaticSolver", - "Max norm of rho is 0", - ablastr::warn_manager::WarnPriority::low - ); - } + const amrex::Real max_norm_b = getMaxNormRho( + amrex::GetVecOfConstPtrs(rho), finest_level, absolute_tolerance); const amrex::LPInfo info; @@ -155,40 +158,85 @@ computeSemiImplicitPhi ( // Use the Multigrid (MLMG) solver but first scale rho appropriately using namespace ablastr::constant::SI; - rho[lev]->mult(-1._rt/ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! + rho[lev]->mult(-1._rt/ep0); -#if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ) - // Use the EB enabled solver. - amrex::MLEBNodeFDLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info + std::unique_ptr linop; + // In the presence of EB or RZ the EB enabled linear solver is used + if (eb_enabled) + { + auto linop_nodelap = std::make_unique(); + linop_nodelap->define( + amrex::Vector{geom[lev]}, + amrex::Vector{grids[lev]}, + amrex::Vector{dmap[lev]}, + info, + amrex::Vector{eb_farray_box_factory.value()[lev]} + ); #if defined(AMREX_USE_EB) - , {eb_farray_box_factory.value()[lev]} -#endif - ); -#else - amrex::MLNodeLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info); + if constexpr (!std::is_same_v) { + // if the EB potential only depends on time, the potential can be passed + // as a float instead of a callable + if (boundary_handler.phi_EB_only_t) { + linop_nodelap->setEBDirichlet(boundary_handler.potential_eb_t(current_time.value())); + } else { + linop_nodelap->setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); + } + } #endif - -#if defined(AMREX_USE_EB) - // if the EB potential only depends on time, the potential can be passed - // as a float instead of a callable - if (boundary_handler.phi_EB_only_t) { - linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value())); + linop_nodelap->setSigma(lev, sigma); + linop = std::move(linop_nodelap); } - else { - linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); + else if (is_rz) + { + auto linop_nodelap = std::make_unique(); + linop_nodelap->define( + amrex::Vector{geom[lev]}, + amrex::Vector{grids[lev]}, + amrex::Vector{dmap[lev]}, + info + ); + linop_nodelap->setRZ(true); + linop_nodelap->setSigma(lev, sigma); + linop = std::move(linop_nodelap); + } + else + { + auto linop_nodelap = std::make_unique(); + linop_nodelap->define( + amrex::Vector{geom[lev]}, + amrex::Vector{grids[lev]}, + amrex::Vector{dmap[lev]}, + info + ); + linop_nodelap->setSigma(lev, sigma); + linop = std::move(linop_nodelap); + } + + // Set domain boundary conditions + if constexpr (std::is_same_v) { + amrex::Array const lobc = {AMREX_D_DECL( + amrex::LinOpBCType::Dirichlet, + amrex::LinOpBCType::Dirichlet, + amrex::LinOpBCType::Dirichlet + )}; + amrex::Array const hibc = lobc; + linop->setDomainBC(lobc, hibc); + } else { + linop->setDomainBC(boundary_handler.lobc, boundary_handler.hibc); } -#endif - linop.setSigma(lev, sigma); // Solve the Poisson equation - linop.setDomainBC( boundary_handler.lobc, boundary_handler.hibc ); -#ifdef WARPX_DIM_RZ - linop.setRZ(true); -#endif - amrex::MLMG mlmg(linop); // actual solver defined here + amrex::MLMG mlmg(*linop); // actual solver defined here mlmg.setVerbose(verbosity); mlmg.setMaxIter(max_iters); - mlmg.setAlwaysUseBNorm(always_use_bnorm); + mlmg.setAlwaysUseBNorm((max_norm_b > 0)); + + const int ng = int(grid_type == utils::enums::GridType::Collocated); // ghost cells + if (ng) { + // In this case, computeE needs to use ghost nodes data. So we + // ask MLMG to fill BC for us after it solves the problem. + mlmg.setFinalFillBC(true); + } // Solve Poisson equation at lev mlmg.solve( {phi[lev]}, {rho[lev]}, @@ -200,44 +248,15 @@ computeSemiImplicitPhi ( // Interpolation from phi[lev] to phi[lev+1] // (This provides both the boundary conditions and initial guess for phi[lev+1]) if (lev < finest_level) { - - // Allocate phi_cp for lev+1 - amrex::BoxArray ba = phi[lev+1]->boxArray(); const amrex::IntVect& refratio = rel_ref_ratio.value()[lev]; - ba.coarsen(refratio); - const int ncomp = linop.getNComp(); - amrex::MultiFab phi_cp(ba, phi[lev+1]->DistributionMap(), ncomp, 1); - - // Copy from phi[lev] to phi_cp (in parallel) - const amrex::IntVect& ng = amrex::IntVect::TheUnitVector(); - const amrex::Periodicity& crse_period = geom[lev].periodicity(); - - ablastr::utils::communication::ParallelCopy( - phi_cp, - *phi[lev], - 0, - 0, - 1, - ng, - ng, - do_single_precision_comms, - crse_period - ); - - // Local interpolation from phi_cp to phi[lev+1] -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (amrex::MFIter mfi(*phi[lev + 1], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - amrex::Array4 const phi_fp_arr = phi[lev + 1]->array(mfi); - amrex::Array4 const phi_cp_arr = phi_cp.array(mfi); - - details::PoissonInterpCPtoFP const interp(phi_fp_arr, phi_cp_arr, refratio); - - amrex::Box const b = mfi.tilebox(phi[lev + 1]->ixType().toIntVect()); - amrex::ParallelFor(b, interp); - } - + const int ncomp = linop->getNComp(); + interpolatePhiBetweenLevels(phi[lev], + phi[lev+1], + geom[lev], + do_single_precision_comms, + refratio, + ncomp, + ng); } // Run additional operations, such as calculation of the E field for embedded boundaries @@ -246,7 +265,7 @@ computeSemiImplicitPhi ( post_phi_calculation.value()(mlmg, lev); } } - + rho[lev]->mult(-ep0); // Multiply rho by epsilon again } // loop over lev(els) } From 7af865054887a4a62ae10987e727c80bf87735c0 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Wed, 23 Oct 2024 10:29:14 -0700 Subject: [PATCH 22/28] replace unavoidable dependence on macros --- Source/ablastr/fields/SemiImplicitPoissonSolver.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/SemiImplicitPoissonSolver.H index bff1a206087..fa7d1e7679b 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/SemiImplicitPoissonSolver.H @@ -164,6 +164,7 @@ computeSemiImplicitPhi ( // In the presence of EB or RZ the EB enabled linear solver is used if (eb_enabled) { +#if defined(AMREX_USE_EB) auto linop_nodelap = std::make_unique(); linop_nodelap->define( amrex::Vector{geom[lev]}, @@ -172,7 +173,6 @@ computeSemiImplicitPhi ( info, amrex::Vector{eb_farray_box_factory.value()[lev]} ); -#if defined(AMREX_USE_EB) if constexpr (!std::is_same_v) { // if the EB potential only depends on time, the potential can be passed // as a float instead of a callable @@ -182,9 +182,9 @@ computeSemiImplicitPhi ( linop_nodelap->setEBDirichlet(boundary_handler.getPhiEB(current_time.value())); } } -#endif linop_nodelap->setSigma(lev, sigma); linop = std::move(linop_nodelap); +#endif } else if (is_rz) { From ed8618ba2dbb3e4f4efc4d82d9892e80b961ea00 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Wed, 23 Oct 2024 12:03:14 -0700 Subject: [PATCH 23/28] set `AMReX_LINEAR_SOLVERS_INCFLO` to ON --- GNUmakefile | 2 +- cmake/dependencies/AMReX.cmake | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GNUmakefile b/GNUmakefile index 1cc78403c7b..68d7ee1870e 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -44,7 +44,7 @@ USE_RZ = FALSE USE_EB = FALSE USE_LINEAR_SOLVERS_EM = TRUE -USE_LINEAR_SOLVERS_INCFLO = FALSE +USE_LINEAR_SOLVERS_INCFLO = TRUE WARPX_HOME := . include $(WARPX_HOME)/Source/Make.WarpX diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 2c4976777e2..ab667cb5e92 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -93,7 +93,7 @@ macro(find_amrex) set(AMReX_PROBINIT OFF CACHE INTERNAL "") set(AMReX_TINY_PROFILE ON CACHE BOOL "") set(AMReX_LINEAR_SOLVERS_EM ON CACHE INTERNAL "") - set(AMReX_LINEAR_SOLVERS_INCFLO OFF CACHE INTERNAL "") + set(AMReX_LINEAR_SOLVERS_INCFLO ON CACHE INTERNAL "") if(WarpX_ASCENT OR WarpX_SENSEI) set(AMReX_GPU_RDC ON CACHE BOOL "") From e41234e586df6ce15a172c7082f771fa3aa7aa38 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Wed, 23 Oct 2024 22:45:37 -0700 Subject: [PATCH 24/28] update CI test to use adiabatic expansion benchmark from https://doi.org/10.1109/TPS.2021.3072353 Signed-off-by: roelof-groenewald --- .../CMakeLists.txt | 8 +- .../semi_implicit_electrostatic/analysis.py | 79 ++++++++- ...t_3d_semi_implicit_electrostatic_picmi.py} | 157 ++++++++++-------- ..._2d_semi_implicit_electrostatic_picmi.json | 15 -- ..._3d_semi_implicit_electrostatic_picmi.json | 15 ++ 5 files changed, 177 insertions(+), 97 deletions(-) rename Examples/Tests/semi_implicit_electrostatic/{inputs_test_2d_semi_implicit_electrostatic_picmi.py => inputs_test_3d_semi_implicit_electrostatic_picmi.py} (55%) delete mode 100644 Regression/Checksum/benchmarks_json/test_2d_semi_implicit_electrostatic_picmi.json create mode 100644 Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json diff --git a/Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt b/Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt index 4ddd7610273..d10c209dbe1 100644 --- a/Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt +++ b/Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt @@ -2,11 +2,11 @@ # add_warpx_test( - test_2d_semi_implicit_electrostatic_picmi # name - 2 # dims + test_3d_semi_implicit_electrostatic_picmi # name + 3 # dims 2 # nprocs - inputs_test_2d_semi_implicit_electrostatic_picmi.py # inputs + inputs_test_3d_semi_implicit_electrostatic_picmi.py # inputs analysis.py # analysis - diags/field_diag000190 # output + diags/field_diag/ # output OFF # dependency ) diff --git a/Examples/Tests/semi_implicit_electrostatic/analysis.py b/Examples/Tests/semi_implicit_electrostatic/analysis.py index f62931ada1c..09b97fce0c3 100755 --- a/Examples/Tests/semi_implicit_electrostatic/analysis.py +++ b/Examples/Tests/semi_implicit_electrostatic/analysis.py @@ -3,10 +3,79 @@ import os import sys -sys.path.insert(1, "../../../../warpx/Regression/Checksum/") -import checksumAPI +import dill +import matplotlib.pyplot as plt +import numpy as np +from openpmd_viewer import OpenPMDTimeSeries +from scipy.interpolate import RegularGridInterpolator -filename = sys.argv[1] +from pywarpx import picmi -test_name = os.path.split(os.getcwd())[1] -checksumAPI.evaluate_checksum(test_name, filename) +constants = picmi.constants + +# load simulation parameters +with open("sim_parameters.dpkl", "rb") as f: + sim = dill.load(f) + +# characteristic expansion time +tau = np.sqrt(constants.m_e * sim.sigma_0**2 / (constants.kb * (sim.T_e + sim.T_i))) + + +def get_analytic_density(r, t): + expansion_factor = 1.0 + t**2 / tau**2 + T = sim.T_e / expansion_factor + sigma = sim.sigma_0 * np.sqrt(expansion_factor) + return sim.n_plasma * (T / sim.T_e) ** 1.5 * np.exp(-(r**2) / (2.0 * sigma**2)) + + +def get_radial_function(field, info): + """Helper function to transform Cartesian data to polar data and average + over theta and phi.""" + r_grid = np.linspace(0, np.max(info.z) - 1e-3, 30) + theta_grid = np.linspace(0, 2 * np.pi, 40, endpoint=False) + phi_grid = np.linspace(0, np.pi, 20) + + r, t, p = np.meshgrid(r_grid, theta_grid, phi_grid, indexing="ij") + x_sp = r * np.sin(p) * np.cos(t) + y_sp = r * np.sin(p) * np.sin(t) + z_sp = r * np.cos(p) + + interp = RegularGridInterpolator((info.x, info.y, info.z), field) + field_sp = interp((x_sp, y_sp, z_sp)) + return r_grid, np.mean(field_sp, axis=(1, 2)) + + +diag_dir = "diags/field_diag" +ts = OpenPMDTimeSeries(diag_dir, check_all_files=True) + +rms_errors = np.zeros(len(ts.iterations)) + +for ii, it in enumerate(ts.iterations): + rho_e, info = ts.get_field(field="rho_electrons", iteration=it) + r_grid, n_e = get_radial_function(-rho_e / constants.q_e, info) + + n_e_analytic = get_analytic_density(r_grid, ts.t[ii]) + rms_errors[ii] = ( + np.sqrt(np.mean(np.sum((n_e - n_e_analytic) ** 2))) / n_e_analytic[0] + ) + + plt.plot(r_grid, n_e_analytic, "k--", alpha=0.6) + plt.plot(r_grid, n_e, label=f"t = {ts.t[ii]*1e6:.1f} $\mu$s") + +print("RMS error (%) in density: ", rms_errors) +assert np.all(rms_errors < 0.06) + +plt.ylabel("$n_e$ (m$^{-3}$)") +plt.xlabel("r (m)") +plt.grid() +plt.legend() +plt.show() + +if len(sys.argv) > 1: + sys.path.insert(1, "../../../../warpx/Regression/Checksum/") + import checksumAPI + + filename = sys.argv[1] + + test_name = os.path.split(os.getcwd())[1] + checksumAPI.evaluate_checksum(test_name, filename, output_format="openpmd") diff --git a/Examples/Tests/semi_implicit_electrostatic/inputs_test_2d_semi_implicit_electrostatic_picmi.py b/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py similarity index 55% rename from Examples/Tests/semi_implicit_electrostatic/inputs_test_2d_semi_implicit_electrostatic_picmi.py rename to Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py index 4bf5dc9feb8..3decca089b0 100644 --- a/Examples/Tests/semi_implicit_electrostatic/inputs_test_2d_semi_implicit_electrostatic_picmi.py +++ b/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py @@ -1,13 +1,26 @@ #!/usr/bin/env python3 # -# --- Test script for the semi-implicit Poisson solver in which a plasma column -# --- expands inside a conducting cylinder. +# --- Test script for the semi-implicit Poisson solver. This test is based on the +# --- adiabatic plasma expansion benchmark from Connor et al. (2021) +# --- doi.org/10.1109/TPS.2021.3072353. +# --- In the benchmark an expanding plasma ball with Gaussian density distribution +# --- in the radial direction is simulated and the time evolution of the +# --- density of the electron species is compared to an approximate analytic solution. +# --- The example is modified slightly in the following ways: +# --- 1) The original example used an electromagnetic solver with absorbing +# --- boundaries while the present case uses an electrostatic solver with +# --- Neumann boundaries. +# --- 2) The domain and plasma parameters for this case has been modified to +# --- set the cell-size and time step such that the explicit electrostatic +# --- solver is unstable. import argparse import sys +import dill import numpy as np from mpi4py import MPI as mpi +from scipy.special import erf from pywarpx import picmi @@ -19,59 +32,60 @@ class PlasmaExpansionSimulation(object): - """Simulation setup for an expanding plasma.""" + """Simulation setup for an expanding plasma ball in 3d.""" - # Plasma parameters - m_ion = 250.0 # Ion mass (electron masses) - n_plasma = 1e18 # Plasma density (m^-3) - T_e = 2.0 # Electron temperature (eV) - T_i = 1.0 # Ion temperature (eV) - - scale_fac = 4.0 # scaling factor used to move between explicit and SI + m_ion = 87.62 # Ion mass (proton masses) - # Plasma size - plasma_radius = 20 # Plasma column radius (Debye lengths) + # Plasma parameters + n_plasma = 5e12 # Plasma density (m^-3) + sigma_0 = 48 # Initial Gaussian distribution standard deviation (Debye lengths) + T_e = 100.0 # Electron temperature (K) + T_i = 1.0 # Ion temperature (K) # Spatial domain - NZ = 256 / scale_fac # Number of cells in x and z directions + R = 200 # Radius of metallic sphere (Debye lengths) + NZ = 72 # Number of cells in each direction # Temporal domain (if not run as a CI test) - LT = 0.1 # Simulation temporal length in ion crossing times + LT = 0.8e-6 # Simulation temporal length (s) # Numerical parameters - NPPC = 750 # Seed number of particles per cell - DZ = 1.0 * scale_fac # Cell size (Debye lengths) - DT = 0.75 # Time step (electron CFL) + NPARTS = 500000 # Seed number of particles + DT = 0.8 # Time step (electron streaming) # Solver parameter - C_SI = 2 # Semi-implicit factor - marginal stability threshold = 1 + C_SI = 2.0 # Semi-implicit factor - def __init__(self, test, verbose): + def __init__(self, verbose): """Get input parameters for the specific case desired.""" - self.test = test - self.verbose = verbose or self.test + self.verbose = verbose # calculate various plasma parameters based on the simulation input self.get_plasma_quantities() - self.dz = self.DZ * self.lambda_e - self.Lz = self.NZ * self.dz + self.R *= self.lambda_e + self.sigma_0 *= self.lambda_e + # self.dz = self.DZ * self.lambda_e + self.dz = 2.0 * self.R / (self.NZ - 4) + self.Lz = self.dz * self.NZ self.dt = self.DT * self.dz / self.v_te - self.plasma_radius *= self.lambda_e + self.total_steps = int(self.LT / self.dt) + self.diag_steps = self.total_steps // 3 - self.total_steps = int(self.LT * self.Lz / (self.v_ti * self.dt)) - self.diag_steps = max(50, self.total_steps // 10) + # dump all the current attributes to a dill pickle file + if comm.rank == 0: + with open("sim_parameters.dpkl", "wb") as f: + dill.dump(self, f) # print out plasma parameters if comm.rank == 0: print( f"Initializing simulation with input parameters:\n" - f"\tT_e = {self.T_e:.1f} eV\n" - f"\tT_i = {self.T_i:.1f} eV\n" + f"\tT_e = {self.T_e:.1f} K\n" + f"\tT_i = {self.T_i:.1f} K\n" f"\tn = {self.n_plasma:.1e} m^-3\n" - f"\tM/m = {self.m_ion:.0f}\n" ) print( f"Plasma parameters:\n" @@ -81,8 +95,8 @@ def __init__(self, test, verbose): ) print( f"Numerical parameters:\n" - f"\tdz = {self.dz:.1e} m\n" - f"\tdt*f_pe = {self.dt*self.f_pe:.2f}\n" + f"\tdz/lambda_e = {self.dz/self.lambda_e:.2f}\n" + f"\tdt*w_pe = {self.dt*self.f_pe*2.0*np.pi:.2f}\n" f"\tdiag steps = {self.diag_steps:d}\n" f"\ttotal steps = {self.total_steps:d}\n" ) @@ -91,7 +105,7 @@ def __init__(self, test, verbose): def get_plasma_quantities(self): """Calculate various plasma parameters based on the simulation input.""" # Ion mass (kg) - self.M = self.m_ion * constants.m_e + self.M = self.m_ion * constants.m_p # Electron plasma frequency (Hz) self.f_pe = np.sqrt( @@ -100,12 +114,12 @@ def get_plasma_quantities(self): # Debye length (m) self.lambda_e = np.sqrt( - constants.ep0 * self.T_e / (self.n_plasma * constants.q_e) + constants.ep0 * constants.kb * self.T_e / (self.n_plasma * constants.q_e**2) ) # Thermal velocities (m/s) from v_th = np.sqrt(kT / m) - self.v_ti = np.sqrt(self.T_i * constants.q_e / self.M) - self.v_te = np.sqrt(self.T_e * constants.q_e / constants.m_e) + self.v_ti = np.sqrt(self.T_i * constants.kb / self.M) + self.v_te = np.sqrt(self.T_e * constants.kb / constants.m_e) def setup_run(self): """Setup simulation components.""" @@ -114,15 +128,15 @@ def setup_run(self): # Set geometry and boundary conditions # ####################################################################### - self.grid = picmi.Cartesian2DGrid( - number_of_cells=[self.NZ, self.NZ], - lower_bound=[-self.Lz / 2.0, -self.Lz / 2.0], - upper_bound=[self.Lz / 2.0, self.Lz / 2.0], - lower_boundary_conditions=["dirichlet"] * 2, - upper_boundary_conditions=["dirichlet"] * 2, - lower_boundary_conditions_particles=["absorbing"] * 2, - upper_boundary_conditions_particles=["absorbing"] * 2, - warpx_max_grid_size=self.NZ, + self.grid = picmi.Cartesian3DGrid( + number_of_cells=[self.NZ] * 3, + lower_bound=[-self.Lz / 2.0] * 3, + upper_bound=[self.Lz / 2.0] * 3, + lower_boundary_conditions=["neumann"] * 3, + upper_boundary_conditions=["neumann"] * 3, + lower_boundary_conditions_particles=["absorbing"] * 3, + upper_boundary_conditions_particles=["absorbing"] * 3, + warpx_max_grid_size=self.NZ // 2, ) simulation.time_step_size = self.dt simulation.max_steps = self.total_steps @@ -130,16 +144,6 @@ def setup_run(self): simulation.particle_shape = 1 simulation.verbose = self.verbose - ####################################################################### - # Insert probe as embedded boundary # - ####################################################################### - - embedded_boundary = picmi.EmbeddedBoundary( - implicit_function=f"(x**2+z**2-{self.Lz/2.-8*self.lambda_e}**2)", - potential=0.0, - ) - simulation.embedded_boundary = embedded_boundary - ####################################################################### # Field solver and external field # ####################################################################### @@ -149,7 +153,7 @@ def setup_run(self): method="Multigrid", warpx_semi_implicit=True, warpx_semi_implicit_factor=self.C_SI, - warpx_self_fields_verbosity=0, + warpx_self_fields_verbosity=self.verbose, ) simulation.solver = solver @@ -157,20 +161,30 @@ def setup_run(self): # Particle types setup # ####################################################################### - density_expr = f"if(x**2+z**2<{self.plasma_radius**2},{self.n_plasma},0)" + total_parts = ( + self.n_plasma + * self.sigma_0**2 + * ( + (2.0 * np.pi) ** 1.5 + * self.sigma_0 + * erf(self.R / (np.sqrt(2) * self.sigma_0)) + + 4.0 * np.pi * self.R * np.exp(-(self.R**2) / (2.0 * self.sigma_0**2)) + ) + ) + self.electrons = picmi.Species( name="electrons", particle_type="electron", - initial_distribution=picmi.AnalyticDistribution( - density_expression=density_expr, - warpx_momentum_spread_expressions=[self.v_te] * 3, - warpx_density_min=self.n_plasma / 10.0, + initial_distribution=picmi.GaussianBunchDistribution( + n_physical_particles=total_parts, + rms_bunch_size=[self.sigma_0] * 3, + rms_velocity=[self.v_te] * 3, ), ) simulation.add_species( self.electrons, layout=picmi.PseudoRandomLayout( - grid=self.grid, n_macroparticles_per_cell=self.NPPC + grid=self.grid, n_macroparticles=self.NPARTS ), ) @@ -178,16 +192,16 @@ def setup_run(self): name="ions", charge="q_e", mass=self.M, - initial_distribution=picmi.AnalyticDistribution( - density_expression=density_expr, - warpx_momentum_spread_expressions=[self.v_ti] * 3, - warpx_density_min=self.n_plasma / 10.0, + initial_distribution=picmi.GaussianBunchDistribution( + n_physical_particles=total_parts, + rms_bunch_size=[self.sigma_0] * 3, + rms_velocity=[self.v_ti] * 3, ), ) simulation.add_species( self.ions, layout=picmi.PseudoRandomLayout( - grid=self.grid, n_macroparticles_per_cell=self.NPPC + grid=self.grid, n_macroparticles=self.NPARTS ), ) @@ -208,6 +222,9 @@ def setup_run(self): "rho_electrons", "rho_ions", ], + write_dir="diags", + warpx_format="openpmd", + warpx_openpmd_backend="h5", ) simulation.add_diagnostic(field_diag) @@ -224,12 +241,6 @@ def setup_run(self): ########################## parser = argparse.ArgumentParser() -parser.add_argument( - "-t", - "--test", - help="toggle whether this script is run as a short CI test", - action="store_true", -) parser.add_argument( "-v", "--verbose", @@ -239,5 +250,5 @@ def setup_run(self): args, left = parser.parse_known_args() sys.argv = sys.argv[:1] + left -run = PlasmaExpansionSimulation(test=args.test, verbose=args.verbose) +run = PlasmaExpansionSimulation(verbose=args.verbose) simulation.step() diff --git a/Regression/Checksum/benchmarks_json/test_2d_semi_implicit_electrostatic_picmi.json b/Regression/Checksum/benchmarks_json/test_2d_semi_implicit_electrostatic_picmi.json deleted file mode 100644 index c13b6d911fc..00000000000 --- a/Regression/Checksum/benchmarks_json/test_2d_semi_implicit_electrostatic_picmi.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - "lev=0": { - "Ex": 1428856.0323152007, - "Ey": 0.0, - "Ez": 1414591.6479792167, - "T_electrons": 1652.0361205183694, - "T_ions": 851.0244739005777, - "jx": 151920.26162133453, - "jy": 453186.62923648424, - "jz": 159070.6572202958, - "phi": 705.4775524885782, - "rho_electrons": 12.100065187070399, - "rho_ions": 12.149198603846399 - } -} diff --git a/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json new file mode 100644 index 00000000000..40731746864 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json @@ -0,0 +1,15 @@ +{ + "lev=0": { + "Ex": 0.0, + "Ey": 0.0, + "Ez": 0.0, + "T_electrons": 109.03294855383206, + "T_ions": 3.068948580858831, + "jx": 141.86426220249493, + "jy": 141.98958981377464, + "jz": 141.74419356241674, + "phi": 312055.82596290566, + "rho_electrons": 0.005758521979913464, + "rho_ions": 0.00685854139567593 + } +} From 03cf3ccb832e42536e24e33360ef313b79475ffb Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Thu, 24 Oct 2024 01:19:57 -0700 Subject: [PATCH 25/28] use EB in CI test --- .../semi_implicit_electrostatic/analysis.py | 6 ++-- ...st_3d_semi_implicit_electrostatic_picmi.py | 29 +++++++++++++------ ..._3d_semi_implicit_electrostatic_picmi.json | 22 +++++++------- .../ElectrostaticSolvers/SemiImplicitES.cpp | 5 ++++ 4 files changed, 39 insertions(+), 23 deletions(-) diff --git a/Examples/Tests/semi_implicit_electrostatic/analysis.py b/Examples/Tests/semi_implicit_electrostatic/analysis.py index 09b97fce0c3..56cad050e2d 100755 --- a/Examples/Tests/semi_implicit_electrostatic/analysis.py +++ b/Examples/Tests/semi_implicit_electrostatic/analysis.py @@ -18,7 +18,7 @@ sim = dill.load(f) # characteristic expansion time -tau = np.sqrt(constants.m_e * sim.sigma_0**2 / (constants.kb * (sim.T_e + sim.T_i))) +tau = sim.sigma_0 * np.sqrt(sim.M / (constants.kb * (sim.T_e + sim.T_i))) def get_analytic_density(r, t): @@ -60,10 +60,10 @@ def get_radial_function(field, info): ) plt.plot(r_grid, n_e_analytic, "k--", alpha=0.6) - plt.plot(r_grid, n_e, label=f"t = {ts.t[ii]*1e6:.1f} $\mu$s") + plt.plot(r_grid, n_e, label=f"t = {ts.t[ii]*1e6:.2f} $\mu$s") print("RMS error (%) in density: ", rms_errors) -assert np.all(rms_errors < 0.06) +assert np.all(rms_errors < 0.05) plt.ylabel("$n_e$ (m$^{-3}$)") plt.xlabel("r (m)") diff --git a/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py b/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py index 3decca089b0..5f17ddc2485 100644 --- a/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py +++ b/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py @@ -8,8 +8,8 @@ # --- density of the electron species is compared to an approximate analytic solution. # --- The example is modified slightly in the following ways: # --- 1) The original example used an electromagnetic solver with absorbing -# --- boundaries while the present case uses an electrostatic solver with -# --- Neumann boundaries. +# --- boundaries while the present case encloses the plasma in a conducting +# --- sphere. # --- 2) The domain and plasma parameters for this case has been modified to # --- set the cell-size and time step such that the explicit electrostatic # --- solver is unstable. @@ -34,27 +34,27 @@ class PlasmaExpansionSimulation(object): """Simulation setup for an expanding plasma ball in 3d.""" - m_ion = 87.62 # Ion mass (proton masses) + m_ion = 25 # Ion mass (electron masses) # Plasma parameters n_plasma = 5e12 # Plasma density (m^-3) - sigma_0 = 48 # Initial Gaussian distribution standard deviation (Debye lengths) + sigma_0 = 22 # Initial Gaussian distribution standard deviation (Debye lengths) T_e = 100.0 # Electron temperature (K) - T_i = 1.0 # Ion temperature (K) + T_i = 10.0 # Ion temperature (K) # Spatial domain - R = 200 # Radius of metallic sphere (Debye lengths) + R = 86 # Radius of metallic sphere (Debye lengths) NZ = 72 # Number of cells in each direction # Temporal domain (if not run as a CI test) - LT = 0.8e-6 # Simulation temporal length (s) + LT = 0.6e-6 # Simulation temporal length (s) # Numerical parameters NPARTS = 500000 # Seed number of particles DT = 0.8 # Time step (electron streaming) # Solver parameter - C_SI = 2.0 # Semi-implicit factor + C_SI = 1.0 # Semi-implicit factor def __init__(self, verbose): """Get input parameters for the specific case desired.""" @@ -73,6 +73,7 @@ def __init__(self, verbose): self.total_steps = int(self.LT / self.dt) self.diag_steps = self.total_steps // 3 + self.total_steps = self.diag_steps * 3 # dump all the current attributes to a dill pickle file if comm.rank == 0: @@ -105,7 +106,7 @@ def __init__(self, verbose): def get_plasma_quantities(self): """Calculate various plasma parameters based on the simulation input.""" # Ion mass (kg) - self.M = self.m_ion * constants.m_p + self.M = self.m_ion * constants.m_e # Electron plasma frequency (Hz) self.f_pe = np.sqrt( @@ -144,6 +145,16 @@ def setup_run(self): simulation.particle_shape = 1 simulation.verbose = self.verbose + ####################################################################### + # Insert spherical boundary as EB # + ####################################################################### + + embedded_boundary = picmi.EmbeddedBoundary( + implicit_function=f"(x**2+y**2+z**2-{self.R**2})", + potential=0.0, + ) + simulation.embedded_boundary = embedded_boundary + ####################################################################### # Field solver and external field # ####################################################################### diff --git a/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json index 40731746864..9199fc5817e 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json @@ -1,15 +1,15 @@ { "lev=0": { - "Ex": 0.0, - "Ey": 0.0, - "Ez": 0.0, - "T_electrons": 109.03294855383206, - "T_ions": 3.068948580858831, - "jx": 141.86426220249493, - "jy": 141.98958981377464, - "jz": 141.74419356241674, - "phi": 312055.82596290566, - "rho_electrons": 0.005758521979913464, - "rho_ions": 0.00685854139567593 + "Ex": 149468.0491281357, + "Ey": 149110.49061700853, + "Ez": 149842.00955840104, + "T_electrons": 279.5716926772583, + "T_ions": 32.15328681363522, + "jx": 31.391568586093264, + "jy": 31.230032200603826, + "jz": 31.519509409154683, + "phi": 2010.2020313518924, + "rho_electrons": 0.00790935978091306, + "rho_ions": 0.008266945197224491 } } diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp index 77023f82fec..1bd2e5fa880 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp @@ -63,6 +63,11 @@ void SemiImplicitES::ComputeSpaceChargeField ( // perform phi calculation computePhi(rho_fp, phi_fp); + + // Compute the electric field. Note that if an EB is used the electric + // field will be calculated in the computePhi call. + const std::array beta = {0._rt}; + if (!EB::enabled()) { computeE( Efield_fp, phi_fp, beta ); } } void SemiImplicitES::computePhi ( From 06d281ac897e235f7d232cca9667d600acd1bef1 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Thu, 24 Oct 2024 08:57:10 -0700 Subject: [PATCH 26/28] update CI checksum values to Azure values --- ..._3d_semi_implicit_electrostatic_picmi.json | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json index 9199fc5817e..b61da644114 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json @@ -1,15 +1,15 @@ { "lev=0": { - "Ex": 149468.0491281357, - "Ey": 149110.49061700853, - "Ez": 149842.00955840104, - "T_electrons": 279.5716926772583, - "T_ions": 32.15328681363522, - "jx": 31.391568586093264, - "jy": 31.230032200603826, - "jz": 31.519509409154683, - "phi": 2010.2020313518924, - "rho_electrons": 0.00790935978091306, - "rho_ions": 0.008266945197224491 + "Ex": 148992.08170563323, + "Ey": 148964.62980386653, + "Ez": 149085.76473996745, + "T_electrons": 279.2796849134268, + "T_ions": 32.09677271761641, + "jx": 31.077856013948246, + "jy": 31.425549493321245, + "jz": 31.424168300110658, + "phi": 2002.2518068289028, + "rho_electrons": 0.007909027251581852, + "rho_ions": 0.008266762306092332 } } From fbfcdaab2dccde3070fa1c57e5fb66dfb07910cf Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Mon, 4 Nov 2024 10:08:50 -0800 Subject: [PATCH 27/28] rename "semi-implicit ES" to "effective potential ES" Signed-off-by: roelof-groenewald --- Docs/source/usage/parameters.rst | 20 +++++++++------- Examples/Tests/CMakeLists.txt | 2 +- .../CMakeLists.txt | 4 ++-- .../analysis.py | 0 ...ffective_potential_electrostatic_picmi.py} | 8 +++---- Python/pywarpx/picmi.py | 24 +++++++++++-------- ...ective_potential_electrostatic_picmi.json} | 0 Source/Diagnostics/Diagnostics.cpp | 4 ++-- .../ElectrostaticSolvers/CMakeLists.txt | 2 +- ...emiImplicitES.H => EffectivePotentialES.H} | 8 +++---- ...mplicitES.cpp => EffectivePotentialES.cpp} | 20 ++++++++-------- .../ElectrostaticSolvers/Make.package | 2 +- Source/Initialization/WarpXInitData.cpp | 4 ++-- Source/Utils/WarpXAlgorithmSelection.H | 2 +- Source/WarpX.cpp | 12 +++++----- ...er.H => EffectivePotentialPoissonSolver.H} | 12 +++++----- 16 files changed, 65 insertions(+), 59 deletions(-) rename Examples/Tests/{semi_implicit_electrostatic => effective_potential_electrostatic}/CMakeLists.txt (61%) rename Examples/Tests/{semi_implicit_electrostatic => effective_potential_electrostatic}/analysis.py (100%) rename Examples/Tests/{semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py => effective_potential_electrostatic/inputs_test_3d_effective_potential_electrostatic_picmi.py} (97%) rename Regression/Checksum/benchmarks_json/{test_3d_semi_implicit_electrostatic_picmi.json => test_3d_effective_potential_electrostatic_picmi.json} (100%) rename Source/FieldSolver/ElectrostaticSolvers/{SemiImplicitES.H => EffectivePotentialES.H} (87%) rename Source/FieldSolver/ElectrostaticSolvers/{SemiImplicitES.cpp => EffectivePotentialES.cpp} (94%) rename Source/ablastr/fields/{SemiImplicitPoissonSolver.H => EffectivePotentialPoissonSolver.H} (97%) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 115a41e1044..22483f509c0 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -228,20 +228,22 @@ Overall simulation parameters \boldsymbol{\nabla}^2 \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi \\ \boldsymbol{\nabla}^2 \boldsymbol{A} = - \mu_0 \boldsymbol{j} \qquad \boldsymbol{B} = \boldsymbol{\nabla}\times\boldsymbol{A} - * ``labframe-semi-implicit``: Poisson's equation is solved with a modified dielectric function - to create a semi-implicit scheme which is robust to the numerical instability seen - in explicit electrostatic PIC when :math:`\Delta t \omega_{pe} > 2`. - If this option is used the additional parameter ``warpx.semi_implicit_factor`` can also be - specified to set the value of :math:`C_{SI}` (default 4). The method is stable for :math:`C_{SI} \geq 1` - regardless of :math:`\Delta t`, however, the larger :math:`C_{SI}` is set, the lower the numerical plasma - frequency will be and therefore care must be taken to not set :math:`C_{SI}` so high that the plasma mode + * ``labframe-effective-potential``: Poisson's equation is solved with a modified dielectric function + (resulting in an "effective potential") to create a semi-implicit scheme which is robust to the + numerical instability seen in explicit electrostatic PIC when :math:`\Delta t \omega_{pe} > 2`. + If this option is used the additional parameter ``warpx.effective_potential_factor`` can also be + specified to set the value of :math:`C_{EP}` (default 4). The method is stable for :math:`C_{EP} \geq 1` + regardless of :math:`\Delta t`, however, the larger :math:`C_{EP}` is set, the lower the numerical plasma + frequency will be and therefore care must be taken to not set it so high that the plasma mode hybridizes with other modes of interest. - Details of the method can be found in Appendix A of :cite:t:`param-Barnes2021`. + Details of the method can be found in Appendix A of :cite:t:`param-Barnes2021` (note that in that paper + the method is referred to as "semi-implicit electrostatic" but here it has been renamed to "effective potential" + to avoid confusion with the semi-implicit method of Chen et al.). In short, the code solves: .. math:: - \boldsymbol{\nabla}\cdot\left(1+\frac{C_{SI}}{4}\sum_{s \in \text{species}}(\omega_{ps}\Delta t)^2 \right)\boldsymbol{\nabla} \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi + \boldsymbol{\nabla}\cdot\left(1+\frac{C_{EP}}{4}\sum_{s \in \text{species}}(\omega_{ps}\Delta t)^2 \right)\boldsymbol{\nabla} \phi = - \rho/\epsilon_0 \qquad \boldsymbol{E} = - \boldsymbol{\nabla}\phi * ``relativistic``: Poisson's equation is solved **for each species** in their respective rest frame. The corresponding field diff --git a/Examples/Tests/CMakeLists.txt b/Examples/Tests/CMakeLists.txt index a3233b4f4a9..c4713123ae6 100644 --- a/Examples/Tests/CMakeLists.txt +++ b/Examples/Tests/CMakeLists.txt @@ -71,7 +71,7 @@ add_subdirectory(restart) add_subdirectory(restart_eb) add_subdirectory(rigid_injection) add_subdirectory(scraping) -add_subdirectory(semi_implicit_electrostatic) +add_subdirectory(effective_potential_electrostatic) add_subdirectory(silver_mueller) add_subdirectory(single_particle) add_subdirectory(space_charge_initialization) diff --git a/Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt b/Examples/Tests/effective_potential_electrostatic/CMakeLists.txt similarity index 61% rename from Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt rename to Examples/Tests/effective_potential_electrostatic/CMakeLists.txt index d10c209dbe1..a6545e8c5f3 100644 --- a/Examples/Tests/semi_implicit_electrostatic/CMakeLists.txt +++ b/Examples/Tests/effective_potential_electrostatic/CMakeLists.txt @@ -2,10 +2,10 @@ # add_warpx_test( - test_3d_semi_implicit_electrostatic_picmi # name + test_3d_effective_potential_electrostatic_picmi # name 3 # dims 2 # nprocs - inputs_test_3d_semi_implicit_electrostatic_picmi.py # inputs + inputs_test_3d_effective_potential_electrostatic_picmi.py # inputs analysis.py # analysis diags/field_diag/ # output OFF # dependency diff --git a/Examples/Tests/semi_implicit_electrostatic/analysis.py b/Examples/Tests/effective_potential_electrostatic/analysis.py similarity index 100% rename from Examples/Tests/semi_implicit_electrostatic/analysis.py rename to Examples/Tests/effective_potential_electrostatic/analysis.py diff --git a/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py b/Examples/Tests/effective_potential_electrostatic/inputs_test_3d_effective_potential_electrostatic_picmi.py similarity index 97% rename from Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py rename to Examples/Tests/effective_potential_electrostatic/inputs_test_3d_effective_potential_electrostatic_picmi.py index 5f17ddc2485..e941179562d 100644 --- a/Examples/Tests/semi_implicit_electrostatic/inputs_test_3d_semi_implicit_electrostatic_picmi.py +++ b/Examples/Tests/effective_potential_electrostatic/inputs_test_3d_effective_potential_electrostatic_picmi.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 # -# --- Test script for the semi-implicit Poisson solver. This test is based on the +# --- Test script for the effective potential Poisson solver. This test is based on the # --- adiabatic plasma expansion benchmark from Connor et al. (2021) # --- doi.org/10.1109/TPS.2021.3072353. # --- In the benchmark an expanding plasma ball with Gaussian density distribution @@ -54,7 +54,7 @@ class PlasmaExpansionSimulation(object): DT = 0.8 # Time step (electron streaming) # Solver parameter - C_SI = 1.0 # Semi-implicit factor + C_EP = 1.0 # Effective potential factor def __init__(self, verbose): """Get input parameters for the specific case desired.""" @@ -162,8 +162,8 @@ def setup_run(self): solver = picmi.ElectrostaticSolver( grid=self.grid, method="Multigrid", - warpx_semi_implicit=True, - warpx_semi_implicit_factor=self.C_SI, + warpx_effective_potential=True, + warpx_effective_potential_factor=self.C_EP, warpx_self_fields_verbosity=self.verbose, ) simulation.solver = solver diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 21c9a62d5d3..8a1b265a9e9 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1883,12 +1883,12 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): warpx_magnetostatic: bool, default=False Whether to use the magnetostatic solver - warpx_semi_implicit: bool, default=False - Whether to use the semi-implicit Poisson solver + warpx_effective_potential: bool, default=False + Whether to use the effective potential Poisson solver (EP-PIC) - warpx_semi_implicit_factor: float, default=4 - If the semi-implicit Poisson solver is used, this sets the value - of C_SI (the method is marginally stable at C_SI = 1) + warpx_effective_potential_factor: float, default=4 + If the effective potential Poisson solver is used, this sets the value + of C_EP (the method is marginally stable at C_EP = 1) warpx_dt_update_interval: string, optional (default = -1) How frequently the timestep is updated. Adaptive timestepping is disabled when this is <= 0. @@ -1906,8 +1906,10 @@ def init(self, kw): self.absolute_tolerance = kw.pop("warpx_absolute_tolerance", None) self.self_fields_verbosity = kw.pop("warpx_self_fields_verbosity", None) self.magnetostatic = kw.pop("warpx_magnetostatic", False) - self.semi_implicit = kw.pop("warpx_semi_implicit", False) - self.semi_implicit_factor = kw.pop("warpx_semi_implicit_factor", None) + self.effective_potential = kw.pop("warpx_effective_potential", False) + self.effective_potential_factor = kw.pop( + "warpx_effective_potential_factor", None + ) self.cfl = kw.pop("warpx_cfl", None) self.dt_update_interval = kw.pop("dt_update_interval", None) self.max_dt = kw.pop("warpx_max_dt", None) @@ -1928,9 +1930,11 @@ def solver_initialize_inputs(self): else: if self.magnetostatic: pywarpx.warpx.do_electrostatic = "labframe-electromagnetostatic" - elif self.semi_implicit: - pywarpx.warpx.do_electrostatic = "labframe-semi-implicit" - pywarpx.warpx.semi_implicit_factor = self.semi_implicit_factor + elif self.effective_potential: + pywarpx.warpx.do_electrostatic = "labframe-effective-potential" + pywarpx.warpx.effective_potential_factor = ( + self.effective_potential_factor + ) else: pywarpx.warpx.do_electrostatic = "labframe" pywarpx.warpx.self_fields_required_precision = self.required_precision diff --git a/Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_effective_potential_electrostatic_picmi.json similarity index 100% rename from Regression/Checksum/benchmarks_json/test_3d_semi_implicit_electrostatic_picmi.json rename to Regression/Checksum/benchmarks_json/test_3d_effective_potential_electrostatic_picmi.json diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index 939a058214b..8a687dcd948 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -78,8 +78,8 @@ Diagnostics::BaseReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE( warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrame || warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameSemiImplicit, - "plot phi only works if do_electrostatic = labframe, do_electrostatic = labframe-electromagnetostatic or do_electrostatic = labframe-semi-implicit"); + warpx.electrostatic_solver_id==ElectrostaticSolverAlgo::LabFrameEffectivePotential, + "plot phi only works if do_electrostatic = labframe, do_electrostatic = labframe-electromagnetostatic or do_electrostatic = labframe-effective-potential"); } // Sanity check if user requests to plot A diff --git a/Source/FieldSolver/ElectrostaticSolvers/CMakeLists.txt b/Source/FieldSolver/ElectrostaticSolvers/CMakeLists.txt index ad3b2cf93a0..db728f6aaba 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/CMakeLists.txt +++ b/Source/FieldSolver/ElectrostaticSolvers/CMakeLists.txt @@ -2,10 +2,10 @@ foreach(D IN LISTS WarpX_DIMS) warpx_set_suffix_dims(SD ${D}) target_sources(lib_${SD} PRIVATE + EffectivePotentialES.cpp ElectrostaticSolver.cpp LabFrameExplicitES.cpp PoissonBoundaryHandler.cpp RelativisticExplicitES.cpp - SemiImplicitES.cpp ) endforeach() diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H b/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H similarity index 87% rename from Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H rename to Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H index 48105d4de9d..be189ec4102 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.H +++ b/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H @@ -13,17 +13,17 @@ #include "ElectrostaticSolver.H" #include -#include +#include -class SemiImplicitES final : public ElectrostaticSolver +class EffectivePotentialES final : public ElectrostaticSolver { public: - SemiImplicitES (int nlevs_max) : ElectrostaticSolver (nlevs_max) { + EffectivePotentialES (int nlevs_max) : ElectrostaticSolver (nlevs_max) { ReadParameters(); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (nlevs_max == 1), - "Semi-implicit electrostatic solver only supports one level at present" + "Effective potential electrostatic solver only supports one level at present" ); } diff --git a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.cpp similarity index 94% rename from Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp rename to Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.cpp index 1bd2e5fa880..e97cb18b5dd 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/SemiImplicitES.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.cpp @@ -7,7 +7,7 @@ * License: BSD-3-Clause-LBNL */ -#include "SemiImplicitES.H" +#include "EffectivePotentialES.H" #include "Fluids/MultiFluidContainer_fwd.H" #include "EmbeddedBoundary/Enabled.H" #include "Fields.H" @@ -17,18 +17,18 @@ using namespace amrex; -void SemiImplicitES::InitData() { +void EffectivePotentialES::InitData() { auto & warpx = WarpX::GetInstance(); m_poisson_boundary_handler->DefinePhiBCs(warpx.Geom(0)); } -void SemiImplicitES::ComputeSpaceChargeField ( +void EffectivePotentialES::ComputeSpaceChargeField ( ablastr::fields::MultiFabRegister& fields, MultiParticleContainer& mpc, [[maybe_unused]] MultiFluidContainer* mfl, int max_level) { - WARPX_PROFILE("SemiImplicitES::ComputeSpaceChargeField"); + WARPX_PROFILE("EffectivePotentialES::ComputeSpaceChargeField"); using ablastr::fields::MultiLevelScalarField; using ablastr::fields::MultiLevelVectorField; @@ -70,7 +70,7 @@ void SemiImplicitES::ComputeSpaceChargeField ( if (!EB::enabled()) { computeE( Efield_fp, phi_fp, beta ); } } -void SemiImplicitES::computePhi ( +void EffectivePotentialES::computePhi ( ablastr::fields::MultiLevelScalarField const& rho, ablastr::fields::MultiLevelScalarField const& phi) const { @@ -87,7 +87,7 @@ void SemiImplicitES::computePhi ( self_fields_verbosity); } -void SemiImplicitES::ComputeSigma (MultiFab& sigma) const +void EffectivePotentialES::ComputeSigma (MultiFab& sigma) const { // Reset sigma to 1 sigma.setVal(1.0_rt); @@ -95,7 +95,7 @@ void SemiImplicitES::ComputeSigma (MultiFab& sigma) const // Get the user set value for C_SI (defaults to 4) amrex::Real C_SI = 4.0; const ParmParse pp_warpx("warpx"); - utils::parser::queryWithParser(pp_warpx, "semi_implicit_factor", C_SI); + utils::parser::queryWithParser(pp_warpx, "effective_potential_factor", C_SI); int const lev = 0; @@ -119,7 +119,7 @@ void SemiImplicitES::ComputeSigma (MultiFab& sigma) const auto& warpx = WarpX::GetInstance(); auto& mypc = warpx.GetPartContainer(); - // The semi-implicit dielectric function is given by + // The effective potential dielectric function is given by // \varepsilon_{SI} = \varepsilon * (1 + \sum_{i in species} C_{SI}*(w_pi * dt)^2/4) // Note the use of the plasma frequency in rad/s (not Hz) and the factor of 1/4, // these choices make it so that C_SI = 1 is the marginal stability threshold. @@ -162,7 +162,7 @@ void SemiImplicitES::ComputeSigma (MultiFab& sigma) const } -void SemiImplicitES::computePhi ( +void EffectivePotentialES::computePhi ( ablastr::fields::MultiLevelScalarField const& rho, ablastr::fields::MultiLevelScalarField const& phi, amrex::MultiFab const& sigma, @@ -233,7 +233,7 @@ void SemiImplicitES::computePhi ( #endif } - ablastr::fields::computeSemiImplicitPhi( + ablastr::fields::computeEffectivePotentialPhi( sorted_rho, sorted_phi, sigma, diff --git a/Source/FieldSolver/ElectrostaticSolvers/Make.package b/Source/FieldSolver/ElectrostaticSolvers/Make.package index 2524fce891a..673f1c8aa7a 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/Make.package +++ b/Source/FieldSolver/ElectrostaticSolvers/Make.package @@ -1,7 +1,7 @@ CEXE_sources += PoissonBoundaryHandler.cpp CEXE_sources += LabFrameExplicitES.cpp CEXE_sources += RelativisticExplicitES.cpp -CEXE_sources += SemiImplicitES.cpp +CEXE_sources += EffectivePotentialES.cpp CEXE_sources += ElectrostaticSolver.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/ElectrostaticSolvers diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 72732b62b97..daecfac8bed 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -278,9 +278,9 @@ WarpX::PrintMainPICparameters () amrex::Print() << "Operation mode: | Electrostatic" << "\n"; amrex::Print() << " | - laboratory frame, electrostatic + magnetostatic" << "\n"; } - else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit){ + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential){ amrex::Print() << "Operation mode: | Electrostatic" << "\n"; - amrex::Print() << " | - laboratory frame, semi-implicit scheme" << "\n"; + amrex::Print() << " | - laboratory frame, effective potential scheme" << "\n"; } else{ amrex::Print() << "Operation mode: | Electromagnetic" << "\n"; diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index a3d587d933b..381bebb1cb6 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -61,7 +61,7 @@ AMREX_ENUM(ElectrostaticSolverAlgo, Relativistic, LabFrameElectroMagnetostatic, LabFrame, - LabFrameSemiImplicit, + LabFrameEffectivePotential, Default = None); AMREX_ENUM(PoissonSolverAlgo, diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 922111f38b7..c79bc6f42b9 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -19,7 +19,7 @@ #include "FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H" #include "FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.H" #include "FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H" -#include "FieldSolver/ElectrostaticSolvers/SemiImplicitES.H" +#include "FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H" @@ -332,10 +332,10 @@ WarpX::WarpX () { m_electrostatic_solver = std::make_unique(nlevs_max); } - // Initialize the semi-implicit electrostatic solver if required - else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) + // Initialize the effective potential electrostatic solver if required + else if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential) { - m_electrostatic_solver = std::make_unique(nlevs_max); + m_electrostatic_solver = std::make_unique(nlevs_max); } else { @@ -2355,7 +2355,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm int rho_ncomps = 0; if( (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame) || (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) || - (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit) || + (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential) || (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) ) { rho_ncomps = ncomps; } @@ -2377,7 +2377,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic || - electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit ) + electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameEffectivePotential ) { const IntVect ngPhi = IntVect( AMREX_D_DECL(1,1,1) ); m_fields.alloc_init(FieldType::phi_fp, lev, amrex::convert(ba, phi_nodal_flag), dm, diff --git a/Source/ablastr/fields/SemiImplicitPoissonSolver.H b/Source/ablastr/fields/EffectivePotentialPoissonSolver.H similarity index 97% rename from Source/ablastr/fields/SemiImplicitPoissonSolver.H rename to Source/ablastr/fields/EffectivePotentialPoissonSolver.H index fa7d1e7679b..5633110e847 100644 --- a/Source/ablastr/fields/SemiImplicitPoissonSolver.H +++ b/Source/ablastr/fields/EffectivePotentialPoissonSolver.H @@ -9,8 +9,8 @@ /* * This file was copied and edited from PoissonSolver.H in the same directory. */ -#ifndef ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H -#define ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H +#ifndef ABLASTR_EFFECTIVE_POTENTIAL_POISSON_SOLVER_H +#define ABLASTR_EFFECTIVE_POTENTIAL_POISSON_SOLVER_H #include #include @@ -98,7 +98,7 @@ template< typename T_FArrayBoxFactory = void > void -computeSemiImplicitPhi ( +computeEffectivePotentialPhi ( ablastr::fields::MultiLevelScalarField const& rho, ablastr::fields::MultiLevelScalarField const& phi, amrex::MultiFab const & sigma, @@ -121,7 +121,7 @@ computeSemiImplicitPhi ( ) { using namespace amrex::literals; - ABLASTR_PROFILE("computeSemiImplicitPhi"); + ABLASTR_PROFILE("computeEffectivePotentialPhi"); if (!rel_ref_ratio.has_value()) { ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(rho.size() == 1u, @@ -138,7 +138,7 @@ computeSemiImplicitPhi ( } ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, - "FFT solver cannot be used with semi-implicit Poisson solve"); + "FFT solver cannot be used with effective potential Poisson solve"); #ifdef WARPX_DIM_RZ constexpr bool is_rz = true; @@ -271,4 +271,4 @@ computeSemiImplicitPhi ( } // namespace ablastr::fields -#endif // ABLASTR_SEMI_IMPLICIT_POISSON_SOLVER_H +#endif // ABLASTR_EFFECTIVE_POTENTIAL_POISSON_SOLVER_H From aed936081aca47e4ca0149329a8254247770527e Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Mon, 4 Nov 2024 12:23:00 -0800 Subject: [PATCH 28/28] Fix missed name updates Co-authored-by: Thomas Marks --- .../ElectrostaticSolvers/EffectivePotentialES.H | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H b/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H index be189ec4102..2ade923211c 100644 --- a/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H +++ b/Source/FieldSolver/ElectrostaticSolvers/EffectivePotentialES.H @@ -7,8 +7,8 @@ * License: BSD-3-Clause-LBNL */ -#ifndef WARPX_SEMIIMPLICITES_H_ -#define WARPX_SEMIIMPLICITES_H_ +#ifndef WARPX_EFFECTIVEPOTENTIALES_H_ +#define WARPX_EFFECTIVEPOTENTIALES_H_ #include "ElectrostaticSolver.H" @@ -38,7 +38,7 @@ public: void ComputeSigma ( amrex::MultiFab& sigma ) const; /** - * Compute the potential `phi` by solving the semi-implicit Poisson equation + * Compute the potential `phi` by solving the semi-implicit Poisson equation using the Effective Potential method * with `rho` as the source. * More specifically, this solves the equation * \f[ @@ -68,4 +68,4 @@ public: }; -#endif // WARPX_SEMIIMPLICITES_H_ +#endif // WARPX_EFFECTIVEPOTENTIALES_H_