From 9476692839160aed8e3f389006d06f1fa9272776 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 11 Nov 2024 08:38:58 -0800 Subject: [PATCH] Implicit add filtering (#5086) This adds filtering to the implicit solver, replacing PR #4600. It is a simple change. All that is needed is adding a call to filter the `Efield_fp` just before the particles are pushed. The current density is already filtered in `SyncCurrentAndRho`. The name of the routine `ApplyFilterJ` was changed to `ApplyFilterMF` since it now has a more general usage. --- Docs/source/developers/fields.rst | 2 +- Examples/Tests/implicit/CMakeLists.txt | 10 ++ ...test_2d_theta_implicit_jfnk_vandb_filtered | 115 ++++++++++++++++++ ...2d_theta_implicit_jfnk_vandb_filtered.json | 31 +++++ Source/Evolve/WarpXEvolve.cpp | 6 +- .../ImplicitSolvers/WarpXImplicitOps.cpp | 3 + Source/Parallelization/WarpXComm.cpp | 34 +++--- Source/WarpX.H | 8 +- 8 files changed, 184 insertions(+), 25 deletions(-) create mode 100644 Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb_filtered create mode 100644 Regression/Checksum/benchmarks_json/test_2d_theta_implicit_jfnk_vandb_filtered.json diff --git a/Docs/source/developers/fields.rst b/Docs/source/developers/fields.rst index 9d980119814..bd6a886ae2a 100644 --- a/Docs/source/developers/fields.rst +++ b/Docs/source/developers/fields.rst @@ -119,7 +119,7 @@ Bilinear filter The multi-pass bilinear filter (applied on the current density) is implemented in ``Source/Filter/``, and class ``WarpX`` holds an instance of this class in member variable ``WarpX::bilinear_filter``. For performance reasons (to avoid creating too many guard cells), this filter is directly applied in communication routines, see ``WarpX::AddCurrentFromFineLevelandSumBoundary`` above and -.. doxygenfunction:: WarpX::ApplyFilterJ(const amrex::Vector, 3>> ¤t, int lev, int idim) +.. doxygenfunction:: WarpX::ApplyFilterMF(const amrex::Vector, 3>> &mfvec, int lev, int idim) .. doxygenfunction:: WarpX::SumBoundaryJ(const amrex::Vector, 3>> ¤t, int lev, int idim, const amrex::Periodicity &period) diff --git a/Examples/Tests/implicit/CMakeLists.txt b/Examples/Tests/implicit/CMakeLists.txt index dabd4de66b8..bf378631e16 100644 --- a/Examples/Tests/implicit/CMakeLists.txt +++ b/Examples/Tests/implicit/CMakeLists.txt @@ -31,6 +31,16 @@ add_warpx_test( OFF # dependency ) +add_warpx_test( + test_2d_theta_implicit_jfnk_vandb_filtered # name + 2 # dims + 2 # nprocs + inputs_test_2d_theta_implicit_jfnk_vandb_filtered # inputs + analysis_vandb_jfnk_2d.py # analysis + diags/diag1000020 # output + OFF # dependency +) + add_warpx_test( test_2d_theta_implicit_jfnk_vandb_picmi # name 2 # dims diff --git a/Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb_filtered b/Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb_filtered new file mode 100644 index 00000000000..4849a5e30a3 --- /dev/null +++ b/Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb_filtered @@ -0,0 +1,115 @@ +################################# +########## CONSTANTS ############ +################################# + +my_constants.n0 = 1.e30 # m^-3 +my_constants.Ti = 100. # eV +my_constants.Te = 100. # eV +my_constants.wpe = q_e*sqrt(n0/(m_e*epsilon0)) +my_constants.de0 = clight/wpe +my_constants.nppcz = 10 # number of particles/cell in z +my_constants.dt = 0.1/wpe # s + +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 20 +amr.n_cell = 40 40 +amr.max_grid_size = 8 +amr.blocking_factor = 8 +amr.max_level = 0 +geometry.dims = 2 +geometry.prob_lo = 0.0 0.0 # physical domain +geometry.prob_hi = 10.0*de0 10.0*de0 + +################################# +####### Boundary condition ###### +################################# +boundary.field_lo = periodic periodic +boundary.field_hi = periodic periodic + +################################# +############ NUMERICS ########### +################################# +warpx.abort_on_warning_threshold = high +warpx.serialize_initial_conditions = 1 +warpx.verbose = 1 +warpx.const_dt = dt +#warpx.cfl = 0.5656 +warpx.use_filter = 1 + +algo.maxwell_solver = Yee +algo.evolve_scheme = "theta_implicit_em" +#algo.evolve_scheme = "semi_implicit_em" + +implicit_evolve.theta = 0.5 +implicit_evolve.max_particle_iterations = 21 +implicit_evolve.particle_tolerance = 1.0e-12 + +#implicit_evolve.nonlinear_solver = "picard" +#picard.verbose = true +#picard.max_iterations = 25 +#picard.relative_tolerance = 0.0 #1.0e-12 +#picard.absolute_tolerance = 0.0 #1.0e-24 +#picard.require_convergence = false + +implicit_evolve.nonlinear_solver = "newton" +newton.verbose = true +newton.max_iterations = 20 +newton.relative_tolerance = 1.0e-12 +newton.absolute_tolerance = 0.0 +newton.require_convergence = false + +gmres.verbose_int = 2 +gmres.max_iterations = 1000 +gmres.relative_tolerance = 1.0e-8 +gmres.absolute_tolerance = 0.0 + +algo.particle_pusher = "boris" +#algo.particle_pusher = "higuera" + +algo.particle_shape = 2 +#algo.current_deposition = "direct" +#algo.current_deposition = "esirkepov" +algo.current_deposition = "villasenor" + +################################# +############ PLASMA ############# +################################# +particles.species_names = electrons protons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = nppcz nppcz +electrons.profile = constant +electrons.density = 1.e30 # number per m^3 +electrons.momentum_distribution_type = "gaussian" +electrons.ux_th = sqrt(Te*q_e/m_e)/clight +electrons.uy_th = sqrt(Te*q_e/m_e)/clight +electrons.uz_th = sqrt(Te*q_e/m_e)/clight + +protons.charge = q_e +protons.mass = m_p +protons.injection_style = "NUniformPerCell" +protons.num_particles_per_cell_each_dim = nppcz nppcz +protons.profile = constant +protons.density = 1.e30 # number per m^3 +protons.momentum_distribution_type = "gaussian" +protons.ux_th = sqrt(Ti*q_e/m_p)/clight +protons.uy_th = sqrt(Ti*q_e/m_p)/clight +protons.uz_th = sqrt(Ti*q_e/m_p)/clight + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 20 +diag1.diag_type = Full +diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE +diag1.electrons.variables = x z w ux uy uz +diag1.protons.variables = x z w ux uy uz + +warpx.reduced_diags_names = particle_energy field_energy +particle_energy.type = ParticleEnergy +particle_energy.intervals = 1 +field_energy.type = FieldEnergy +field_energy.intervals = 1 diff --git a/Regression/Checksum/benchmarks_json/test_2d_theta_implicit_jfnk_vandb_filtered.json b/Regression/Checksum/benchmarks_json/test_2d_theta_implicit_jfnk_vandb_filtered.json new file mode 100644 index 00000000000..d342c49e2fd --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_2d_theta_implicit_jfnk_vandb_filtered.json @@ -0,0 +1,31 @@ +{ + "lev=0": { + "Bx": 65625.24877705125, + "By": 71913.65275407257, + "Bz": 59768.79247890749, + "Ex": 56341360261928.086, + "Ey": 13926508614721.855, + "Ez": 56508162715968.17, + "divE": 5.5816922509658905e+22, + "jx": 1.8114330881270456e+19, + "jy": 2.0727708668063334e+19, + "jz": 1.7843765469944717e+19, + "rho": 494213515033.04443 + }, + "electrons": { + "particle_momentum_x": 4.888781979240524e-19, + "particle_momentum_y": 4.879904653089102e-19, + "particle_momentum_z": 4.878388335258947e-19, + "particle_position_x": 0.0042514822919144084, + "particle_position_y": 0.0042515394083575886, + "particle_weight": 2823958719279159.5 + }, + "protons": { + "particle_momentum_x": 2.0873319751377048e-17, + "particle_momentum_y": 2.0858882863041667e-17, + "particle_momentum_z": 2.0877426824914187e-17, + "particle_position_x": 0.004251275869325256, + "particle_position_y": 0.0042512738905204584, + "particle_weight": 2823958719279159.5 + } +} diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index a685afd28e7..e9540be3da7 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -616,7 +616,7 @@ void WarpX::SyncCurrentAndRho () // TODO This works only without mesh refinement const int lev = 0; if (use_filter) { - ApplyFilterJ(m_fields.get_mr_levels_alldirs(FieldType::current_fp_vay, finest_level), lev); + ApplyFilterMF(m_fields.get_mr_levels_alldirs(FieldType::current_fp_vay, finest_level), lev); } } } @@ -875,7 +875,7 @@ WarpX::OneStep_sub1 (Real cur_time) m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level), fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); if (use_filter) { - ApplyFilterJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev); + ApplyFilterMF( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev); } SumBoundaryJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), @@ -953,7 +953,7 @@ WarpX::OneStep_sub1 (Real cur_time) m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level), fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); if (use_filter) { - ApplyFilterJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev); + ApplyFilterMF( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev); } SumBoundaryJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev, Geom(fine_lev).periodicity()); ApplyFilterandSumBoundaryRho( diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp b/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp index 806c3412990..3cf42f18456 100644 --- a/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp +++ b/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp @@ -52,8 +52,11 @@ WarpX::ImplicitPreRHSOp ( amrex::Real a_cur_time, bool a_from_jacobian ) { using namespace amrex::literals; + using warpx::fields::FieldType; amrex::ignore_unused( a_full_dt, a_nl_iter, a_from_jacobian ); + if (use_filter) { ApplyFilterMF(m_fields.get_mr_levels_alldirs(FieldType::Efield_fp, finest_level), 0); } + // Advance the particle positions by 1/2 dt, // particle velocities by dt, then take average of old and new v, // deposit currents, giving J at n+1/2 diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index d64632d964a..a0ae7ed67e9 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1195,7 +1195,7 @@ WarpX::SyncCurrent (const std::string& current_fp_string) ablastr::fields::MultiLevelVectorField const& J_cp = m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level); if (use_filter) { - ApplyFilterJ(J_cp, lev+1, idim); + ApplyFilterMF(J_cp, lev+1, idim); } SumBoundaryJ(J_cp, lev+1, idim, period); } @@ -1232,7 +1232,7 @@ WarpX::SyncCurrent (const std::string& current_fp_string) if (use_filter) { - ApplyFilterJ(J_fp, lev, idim); + ApplyFilterMF(J_fp, lev, idim); } SumBoundaryJ(J_fp, lev, idim, period); } @@ -1354,32 +1354,32 @@ void WarpX::RestrictCurrentFromFineToCoarsePatch ( ablastr::coarsen::average::Coarsen(*crse[2], *fine[2], refinement_ratio ); } -void WarpX::ApplyFilterJ ( - const ablastr::fields::MultiLevelVectorField& current, +void WarpX::ApplyFilterMF ( + const ablastr::fields::MultiLevelVectorField& mfvec, const int lev, const int idim) { using ablastr::fields::Direction; - amrex::MultiFab& J = *current[lev][Direction{idim}]; + amrex::MultiFab& mf = *mfvec[lev][Direction{idim}]; - const int ncomp = J.nComp(); - const amrex::IntVect ngrow = J.nGrowVect(); - amrex::MultiFab Jf(J.boxArray(), J.DistributionMap(), ncomp, ngrow); - bilinear_filter.ApplyStencil(Jf, J, lev); + const int ncomp = mf.nComp(); + const amrex::IntVect ngrow = mf.nGrowVect(); + amrex::MultiFab mf_filtered(mf.boxArray(), mf.DistributionMap(), ncomp, ngrow); + bilinear_filter.ApplyStencil(mf_filtered, mf, lev); const int srccomp = 0; const int dstcomp = 0; - amrex::MultiFab::Copy(J, Jf, srccomp, dstcomp, ncomp, ngrow); + amrex::MultiFab::Copy(mf, mf_filtered, srccomp, dstcomp, ncomp, ngrow); } -void WarpX::ApplyFilterJ ( - const ablastr::fields::MultiLevelVectorField& current, +void WarpX::ApplyFilterMF ( + const ablastr::fields::MultiLevelVectorField& mfvec, const int lev) { for (int idim=0; idim<3; ++idim) { - ApplyFilterJ(current, lev, idim); + ApplyFilterMF(mfvec, lev, idim); } } @@ -1457,7 +1457,7 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary ( if (use_filter) { - ApplyFilterJ(J_fp, lev); + ApplyFilterMF(J_fp, lev); } SumBoundaryJ(J_fp, lev, period); @@ -1476,8 +1476,8 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary ( if (use_filter && J_buffer[lev+1][idim]) { - ApplyFilterJ(J_cp, lev+1, idim); - ApplyFilterJ(J_buffer, lev+1, idim); + ApplyFilterMF(J_cp, lev+1, idim); + ApplyFilterMF(J_buffer, lev+1, idim); MultiFab::Add( *J_buffer[lev+1][idim], *J_cp[lev+1][idim], @@ -1491,7 +1491,7 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary ( } else if (use_filter) // but no buffer { - ApplyFilterJ(J_cp, lev+1, idim); + ApplyFilterMF(J_cp, lev+1, idim); ablastr::utils::communication::ParallelAdd( mf, *J_cp[lev+1][idim], 0, 0, diff --git a/Source/WarpX.H b/Source/WarpX.H index 4dc3ab8c8be..da1a4b5a269 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1216,12 +1216,12 @@ private: int lev); void StoreCurrent (int lev); void RestoreCurrent (int lev); - void ApplyFilterJ ( - const ablastr::fields::MultiLevelVectorField& current, + void ApplyFilterMF ( + const ablastr::fields::MultiLevelVectorField& mfvec, int lev, int idim); - void ApplyFilterJ ( - const ablastr::fields::MultiLevelVectorField& current, + void ApplyFilterMF ( + const ablastr::fields::MultiLevelVectorField& mfvec, int lev); void SumBoundaryJ ( const ablastr::fields::MultiLevelVectorField& current,