Implicit add filtering

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<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &current, int lev, int idim)
+.. doxygenfunction:: WarpX::ApplyFilterMF(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &mfvec, int lev, int idim)
 
 .. doxygenfunction:: WarpX::SumBoundaryJ(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &current, int lev, int idim, const amrex::Periodicity &period)
 

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

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

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
+  }
+}

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(

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

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,

Source/WarpX.H

@@ -1205,12 +1205,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,