Skip to content

Commit 9476692

Browse files
dpgrotedpgrote
authored
Implicit add filtering (BLAST-WarpX#5086)
This adds filtering to the implicit solver, replacing PR BLAST-WarpX#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.
1 parent c36d6aa commit 9476692

File tree

8 files changed

+184
-25
lines changed

8 files changed

+184
-25
lines changed

‎Docs/source/developers/fields.rst‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ Bilinear filter
119119

120120
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
121121

122-
.. doxygenfunction:: WarpX::ApplyFilterJ(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &current, int lev, int idim)
122+
.. doxygenfunction:: WarpX::ApplyFilterMF(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &mfvec, int lev, int idim)
123123

124124
.. doxygenfunction:: WarpX::SumBoundaryJ(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3>> &current, int lev, int idim, const amrex::Periodicity &period)
125125

‎Examples/Tests/implicit/CMakeLists.txt‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,16 @@ add_warpx_test(
3131
OFF # dependency
3232
)
3333

34+
add_warpx_test(
35+
test_2d_theta_implicit_jfnk_vandb_filtered # name
36+
2 # dims
37+
2 # nprocs
38+
inputs_test_2d_theta_implicit_jfnk_vandb_filtered # inputs
39+
analysis_vandb_jfnk_2d.py # analysis
40+
diags/diag1000020 # output
41+
OFF # dependency
42+
)
43+
3444
add_warpx_test(
3545
test_2d_theta_implicit_jfnk_vandb_picmi # name
3646
2 # dims

‎Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb_filtered‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
#################################
2+
########## CONSTANTS ############
3+
#################################
4+
5+
my_constants.n0 = 1.e30 # m^-3
6+
my_constants.Ti = 100. # eV
7+
my_constants.Te = 100. # eV
8+
my_constants.wpe = q_e*sqrt(n0/(m_e*epsilon0))
9+
my_constants.de0 = clight/wpe
10+
my_constants.nppcz = 10 # number of particles/cell in z
11+
my_constants.dt = 0.1/wpe # s
12+
13+
#################################
14+
####### GENERAL PARAMETERS ######
15+
#################################
16+
max_step = 20
17+
amr.n_cell = 40 40
18+
amr.max_grid_size = 8
19+
amr.blocking_factor = 8
20+
amr.max_level = 0
21+
geometry.dims = 2
22+
geometry.prob_lo = 0.0 0.0 # physical domain
23+
geometry.prob_hi = 10.0*de0 10.0*de0
24+
25+
#################################
26+
####### Boundary condition ######
27+
#################################
28+
boundary.field_lo = periodic periodic
29+
boundary.field_hi = periodic periodic
30+
31+
#################################
32+
############ NUMERICS ###########
33+
#################################
34+
warpx.abort_on_warning_threshold = high
35+
warpx.serialize_initial_conditions = 1
36+
warpx.verbose = 1
37+
warpx.const_dt = dt
38+
#warpx.cfl = 0.5656
39+
warpx.use_filter = 1
40+
41+
algo.maxwell_solver = Yee
42+
algo.evolve_scheme = "theta_implicit_em"
43+
#algo.evolve_scheme = "semi_implicit_em"
44+
45+
implicit_evolve.theta = 0.5
46+
implicit_evolve.max_particle_iterations = 21
47+
implicit_evolve.particle_tolerance = 1.0e-12
48+
49+
#implicit_evolve.nonlinear_solver = "picard"
50+
#picard.verbose = true
51+
#picard.max_iterations = 25
52+
#picard.relative_tolerance = 0.0 #1.0e-12
53+
#picard.absolute_tolerance = 0.0 #1.0e-24
54+
#picard.require_convergence = false
55+
56+
implicit_evolve.nonlinear_solver = "newton"
57+
newton.verbose = true
58+
newton.max_iterations = 20
59+
newton.relative_tolerance = 1.0e-12
60+
newton.absolute_tolerance = 0.0
61+
newton.require_convergence = false
62+
63+
gmres.verbose_int = 2
64+
gmres.max_iterations = 1000
65+
gmres.relative_tolerance = 1.0e-8
66+
gmres.absolute_tolerance = 0.0
67+
68+
algo.particle_pusher = "boris"
69+
#algo.particle_pusher = "higuera"
70+
71+
algo.particle_shape = 2
72+
#algo.current_deposition = "direct"
73+
#algo.current_deposition = "esirkepov"
74+
algo.current_deposition = "villasenor"
75+
76+
#################################
77+
############ PLASMA #############
78+
#################################
79+
particles.species_names = electrons protons
80+
81+
electrons.charge = -q_e
82+
electrons.mass = m_e
83+
electrons.injection_style = "NUniformPerCell"
84+
electrons.num_particles_per_cell_each_dim = nppcz nppcz
85+
electrons.profile = constant
86+
electrons.density = 1.e30 # number per m^3
87+
electrons.momentum_distribution_type = "gaussian"
88+
electrons.ux_th = sqrt(Te*q_e/m_e)/clight
89+
electrons.uy_th = sqrt(Te*q_e/m_e)/clight
90+
electrons.uz_th = sqrt(Te*q_e/m_e)/clight
91+
92+
protons.charge = q_e
93+
protons.mass = m_p
94+
protons.injection_style = "NUniformPerCell"
95+
protons.num_particles_per_cell_each_dim = nppcz nppcz
96+
protons.profile = constant
97+
protons.density = 1.e30 # number per m^3
98+
protons.momentum_distribution_type = "gaussian"
99+
protons.ux_th = sqrt(Ti*q_e/m_p)/clight
100+
protons.uy_th = sqrt(Ti*q_e/m_p)/clight
101+
protons.uz_th = sqrt(Ti*q_e/m_p)/clight
102+
103+
# Diagnostics
104+
diagnostics.diags_names = diag1
105+
diag1.intervals = 20
106+
diag1.diag_type = Full
107+
diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
108+
diag1.electrons.variables = x z w ux uy uz
109+
diag1.protons.variables = x z w ux uy uz
110+
111+
warpx.reduced_diags_names = particle_energy field_energy
112+
particle_energy.type = ParticleEnergy
113+
particle_energy.intervals = 1
114+
field_energy.type = FieldEnergy
115+
field_energy.intervals = 1

‎Regression/Checksum/benchmarks_json/test_2d_theta_implicit_jfnk_vandb_filtered.json‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
{
2+
"lev=0": {
3+
"Bx": 65625.24877705125,
4+
"By": 71913.65275407257,
5+
"Bz": 59768.79247890749,
6+
"Ex": 56341360261928.086,
7+
"Ey": 13926508614721.855,
8+
"Ez": 56508162715968.17,
9+
"divE": 5.5816922509658905e+22,
10+
"jx": 1.8114330881270456e+19,
11+
"jy": 2.0727708668063334e+19,
12+
"jz": 1.7843765469944717e+19,
13+
"rho": 494213515033.04443
14+
},
15+
"electrons": {
16+
"particle_momentum_x": 4.888781979240524e-19,
17+
"particle_momentum_y": 4.879904653089102e-19,
18+
"particle_momentum_z": 4.878388335258947e-19,
19+
"particle_position_x": 0.0042514822919144084,
20+
"particle_position_y": 0.0042515394083575886,
21+
"particle_weight": 2823958719279159.5
22+
},
23+
"protons": {
24+
"particle_momentum_x": 2.0873319751377048e-17,
25+
"particle_momentum_y": 2.0858882863041667e-17,
26+
"particle_momentum_z": 2.0877426824914187e-17,
27+
"particle_position_x": 0.004251275869325256,
28+
"particle_position_y": 0.0042512738905204584,
29+
"particle_weight": 2823958719279159.5
30+
}
31+
}

‎Source/Evolve/WarpXEvolve.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -616,7 +616,7 @@ void WarpX::SyncCurrentAndRho ()
616616
// TODO This works only without mesh refinement
617617
const int lev = 0;
618618
if (use_filter) {
619-
ApplyFilterJ(m_fields.get_mr_levels_alldirs(FieldType::current_fp_vay, finest_level), lev);
619+
ApplyFilterMF(m_fields.get_mr_levels_alldirs(FieldType::current_fp_vay, finest_level), lev);
620620
}
621621
}
622622
}
@@ -875,7 +875,7 @@ WarpX::OneStep_sub1 (Real cur_time)
875875
m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level), fine_lev);
876876
RestrictRhoFromFineToCoarsePatch(fine_lev);
877877
if (use_filter) {
878-
ApplyFilterJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev);
878+
ApplyFilterMF( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev);
879879
}
880880
SumBoundaryJ(
881881
m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level),
@@ -953,7 +953,7 @@ WarpX::OneStep_sub1 (Real cur_time)
953953
m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level), fine_lev);
954954
RestrictRhoFromFineToCoarsePatch(fine_lev);
955955
if (use_filter) {
956-
ApplyFilterJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev);
956+
ApplyFilterMF( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev);
957957
}
958958
SumBoundaryJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev, Geom(fine_lev).periodicity());
959959
ApplyFilterandSumBoundaryRho(

‎Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,8 +52,11 @@ WarpX::ImplicitPreRHSOp ( amrex::Real a_cur_time,
5252
bool a_from_jacobian )
5353
{
5454
using namespace amrex::literals;
55+
using warpx::fields::FieldType;
5556
amrex::ignore_unused( a_full_dt, a_nl_iter, a_from_jacobian );
5657

58+
if (use_filter) { ApplyFilterMF(m_fields.get_mr_levels_alldirs(FieldType::Efield_fp, finest_level), 0); }
59+
5760
// Advance the particle positions by 1/2 dt,
5861
// particle velocities by dt, then take average of old and new v,
5962
// deposit currents, giving J at n+1/2

‎Source/Parallelization/WarpXComm.cpp‎

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1195,7 +1195,7 @@ WarpX::SyncCurrent (const std::string& current_fp_string)
11951195
ablastr::fields::MultiLevelVectorField const& J_cp = m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level);
11961196
if (use_filter)
11971197
{
1198-
ApplyFilterJ(J_cp, lev+1, idim);
1198+
ApplyFilterMF(J_cp, lev+1, idim);
11991199
}
12001200
SumBoundaryJ(J_cp, lev+1, idim, period);
12011201
}
@@ -1232,7 +1232,7 @@ WarpX::SyncCurrent (const std::string& current_fp_string)
12321232

12331233
if (use_filter)
12341234
{
1235-
ApplyFilterJ(J_fp, lev, idim);
1235+
ApplyFilterMF(J_fp, lev, idim);
12361236
}
12371237
SumBoundaryJ(J_fp, lev, idim, period);
12381238
}
@@ -1354,32 +1354,32 @@ void WarpX::RestrictCurrentFromFineToCoarsePatch (
13541354
ablastr::coarsen::average::Coarsen(*crse[2], *fine[2], refinement_ratio );
13551355
}
13561356

1357-
void WarpX::ApplyFilterJ (
1358-
const ablastr::fields::MultiLevelVectorField& current,
1357+
void WarpX::ApplyFilterMF (
1358+
const ablastr::fields::MultiLevelVectorField& mfvec,
13591359
const int lev,
13601360
const int idim)
13611361
{
13621362
using ablastr::fields::Direction;
13631363

1364-
amrex::MultiFab& J = *current[lev][Direction{idim}];
1364+
amrex::MultiFab& mf = *mfvec[lev][Direction{idim}];
13651365

1366-
const int ncomp = J.nComp();
1367-
const amrex::IntVect ngrow = J.nGrowVect();
1368-
amrex::MultiFab Jf(J.boxArray(), J.DistributionMap(), ncomp, ngrow);
1369-
bilinear_filter.ApplyStencil(Jf, J, lev);
1366+
const int ncomp = mf.nComp();
1367+
const amrex::IntVect ngrow = mf.nGrowVect();
1368+
amrex::MultiFab mf_filtered(mf.boxArray(), mf.DistributionMap(), ncomp, ngrow);
1369+
bilinear_filter.ApplyStencil(mf_filtered, mf, lev);
13701370

13711371
const int srccomp = 0;
13721372
const int dstcomp = 0;
1373-
amrex::MultiFab::Copy(J, Jf, srccomp, dstcomp, ncomp, ngrow);
1373+
amrex::MultiFab::Copy(mf, mf_filtered, srccomp, dstcomp, ncomp, ngrow);
13741374
}
13751375

1376-
void WarpX::ApplyFilterJ (
1377-
const ablastr::fields::MultiLevelVectorField& current,
1376+
void WarpX::ApplyFilterMF (
1377+
const ablastr::fields::MultiLevelVectorField& mfvec,
13781378
const int lev)
13791379
{
13801380
for (int idim=0; idim<3; ++idim)
13811381
{
1382-
ApplyFilterJ(current, lev, idim);
1382+
ApplyFilterMF(mfvec, lev, idim);
13831383
}
13841384
}
13851385

@@ -1457,7 +1457,7 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary (
14571457

14581458
if (use_filter)
14591459
{
1460-
ApplyFilterJ(J_fp, lev);
1460+
ApplyFilterMF(J_fp, lev);
14611461
}
14621462
SumBoundaryJ(J_fp, lev, period);
14631463

@@ -1476,8 +1476,8 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary (
14761476

14771477
if (use_filter && J_buffer[lev+1][idim])
14781478
{
1479-
ApplyFilterJ(J_cp, lev+1, idim);
1480-
ApplyFilterJ(J_buffer, lev+1, idim);
1479+
ApplyFilterMF(J_cp, lev+1, idim);
1480+
ApplyFilterMF(J_buffer, lev+1, idim);
14811481

14821482
MultiFab::Add(
14831483
*J_buffer[lev+1][idim], *J_cp[lev+1][idim],
@@ -1491,7 +1491,7 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary (
14911491
}
14921492
else if (use_filter) // but no buffer
14931493
{
1494-
ApplyFilterJ(J_cp, lev+1, idim);
1494+
ApplyFilterMF(J_cp, lev+1, idim);
14951495

14961496
ablastr::utils::communication::ParallelAdd(
14971497
mf, *J_cp[lev+1][idim], 0, 0,

‎Source/WarpX.H‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1216,12 +1216,12 @@ private:
12161216
int lev);
12171217
void StoreCurrent (int lev);
12181218
void RestoreCurrent (int lev);
1219-
void ApplyFilterJ (
1220-
const ablastr::fields::MultiLevelVectorField& current,
1219+
void ApplyFilterMF (
1220+
const ablastr::fields::MultiLevelVectorField& mfvec,
12211221
int lev,
12221222
int idim);
1223-
void ApplyFilterJ (
1224-
const ablastr::fields::MultiLevelVectorField& current,
1223+
void ApplyFilterMF (
1224+
const ablastr::fields::MultiLevelVectorField& mfvec,
12251225
int lev);
12261226
void SumBoundaryJ (
12271227
const ablastr::fields::MultiLevelVectorField& current,

0 commit comments

Comments
 (0)