Skip to content

Commit e044430

Browse files
authored
Add the possibility to disable Schwinger in part of the domain (#1524)
* Add the possibility to disable Schwinger in part of the domain * Update checksum benchmarks * Only query ymin and ymax in 3D
1 parent 89c3715 commit e044430

File tree

8 files changed

+165
-74
lines changed

8 files changed

+165
-74
lines changed

Docs/source/running_cpp/parameters.rst

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -332,14 +332,12 @@ Particle initialization
332332
The mass of one `physical` particle of this species.
333333
If ``species_type`` is specified, the mass will be set to the physical value and ``mass`` is optional.
334334

335-
* ``<species_name>.xmin,ymin,zmin`` (`float`) optional (default unlimited)
336-
When ``<species_name>.xmin`` and ``<species_name>.xmax`` (see below) are set, they delimit the region within which particles are injected.
335+
* ``<species_name>.xmin,ymin,zmin`` and ``<species_name>.xmax,ymax,zmax`` (`float`) optional (default unlimited)
336+
When ``<species_name>.xmin`` and ``<species_name>.xmax`` are set, they delimit the region within which particles are injected.
337337
The WarpXParser (see :ref:`running-cpp-parameters-parser`) is used for the right-hand-side, so expressions like ``<species_name>.xmin = "2.+1."`` and/or using user-defined constants are accepted.
338338
The same is applicable in the other directions.
339339
If periodic boundary conditions are used in direction ``i``, then the default (i.e. if the range is not specified) range will be the simulation box, ``[geometry.prob_hi[i], geometry.prob_lo[i]]``.
340340

341-
* ``<species_name>.xmax,ymax,zmax`` (`float`) optional (default unlimited)
342-
343341
* ``<species_name>.injection_style`` (`string`)
344342
Determines how the particles will be injected in the simulation.
345343
The options are:
@@ -1978,6 +1976,11 @@ Lookup tables store pre-computed values for functions used by the QED modules.
19781976
This value should correspond to the typical transverse extent for which the EM field has a very high value
19791977
(e.g. the beam waist for a focused laser beam).
19801978

1979+
* ``qed_schwinger.xmin,ymin,zmin`` and ``qed_schwinger.xmax,ymax,zmax`` (`float`) optional (default unlimited)
1980+
When ``qed_schwinger.xmin`` and ``qed_schwinger.xmax`` are set, they delimit the region within
1981+
which Schwinger pairs can be created.
1982+
The same is applicable in the other directions.
1983+
19811984
* ``qed_schwinger.threshold_poisson_gaussian`` (`integer`) optional (default `25`)
19821985
If the expected number of physical pairs created in a cell at a given timestep is smaller than this threshold,
19831986
a Poisson distribution is used to draw the actual number of physical pairs created.

Examples/Modules/qed/schwinger/analysis_schwinger.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@
5353
Bx_test = 1679288857.0516706
5454
By_test = 525665014.1557486
5555
Bz_test = 1836353079.9561853
56+
dV = dV/2. # Schwinger is only activated in part of the simulation domain
5657
elif test_number == '3':
5758
# Third Schwinger test with intermediate electric field such that average created pair per cell
5859
# is 1. A Poisson distribution is used to obtain the weights of the particles.
@@ -63,6 +64,7 @@
6364
# case.
6465
Ez_test = 2.5e+20
6566
By_test = 833910140000.
67+
dV = dV*(3./4.)**2. # Schwinger is only activated in part of the simulation domain
6668
else:
6769
assert(False)
6870

Lines changed: 26 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -1,35 +1,35 @@
11
{
22
"ele_schwinger": {
3-
"particle_cpu": 1024.0,
4-
"particle_id": 1573888.0,
5-
"particle_momentum_x": 1.2538560192562119145907523613836e-09,
6-
"particle_momentum_y": 7.98841234718825980904544944075e-11,
7-
"particle_momentum_z": 2.802314569448479707680714186416e-10,
8-
"particle_position_x": 0.0005120000000000001,
9-
"particle_position_y": 0.000512,
10-
"particle_position_z": 0.000512,
11-
"particle_weight": 8.3009747141536481280e+19
3+
"particle_cpu": 512.0,
4+
"particle_id": 393728.0,
5+
"particle_momentum_x": 0.0000000005389357130722444579163396297604,
6+
"particle_momentum_y": 0.0000000000380620585874242066235441537959,
7+
"particle_momentum_z": 0.0000000001303705903189427521990174631828,
8+
"particle_position_x": 0.0001414301708416780501906262479394627007,
9+
"particle_position_y": 0.0002559999999999999342054080031516605231,
10+
"particle_position_z": 0.0002559999999999999884155166274268822235,
11+
"particle_weight": 41504873569658028032.0
1212
},
1313
"lev=0": {
1414
"Bx": 3439183579241.8213,
15-
"By": 1076561948990.9731,
16-
"Bz": 3760851107750.2676,
17-
"Ex": 1.61382622190626294479192064e+26,
18-
"Ey": 1.0214830510117843745374208e+25,
19-
"Ez": 3.5684390178813696941228032e+25,
20-
"jx": 1.4539783886540389390252047859712e+31,
21-
"jy": 9.20294522889443324046591655936e+29,
22-
"jz": 3.214947991714962836023450533888e+30
15+
"By": 3467279530258676.5,
16+
"Bz": 1160277379606216.0,
17+
"Ex": 80690799093154680420696064.0,
18+
"Ey": 5505248329779462566576128.0,
19+
"Ez": 17842195088929589557198848.0,
20+
"jx": 7269891943075730389116072755200.0,
21+
"jy": 495989618234384662296063377408.0,
22+
"jz": 1607473995814483300569405456384.0
2323
},
2424
"pos_schwinger": {
25-
"particle_cpu": 1024.0,
26-
"particle_id": 2622464.0,
27-
"particle_momentum_x": 1.2705199397782006324793234588533e-09,
28-
"particle_momentum_y": 8.00263182163650910948918243636e-11,
29-
"particle_momentum_z": 2.809319203195289649657045801174e-10,
30-
"particle_position_x": 0.0005120000000000001,
31-
"particle_position_y": 0.000512,
32-
"particle_position_z": 0.000512,
33-
"particle_weight": 8.3009747141536481280e+19
25+
"particle_cpu": 512.0,
26+
"particle_id": 655872.0,
27+
"particle_momentum_x": 0.0000000005389349271464131667351435623151,
28+
"particle_momentum_y": 0.0000000000364586108023806536132049519030,
29+
"particle_momentum_z": 0.0000000001311702385821036629827073120337,
30+
"particle_position_x": 0.0001280000000000000213128126258510519619,
31+
"particle_position_y": 0.0002559999999999999884155166274268822235,
32+
"particle_position_z": 0.0002559999999999999884155166274268822235,
33+
"particle_weight": 41504873569658028032.0
3434
}
3535
}
Lines changed: 23 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,35 +1,35 @@
11
{
22
"ele_schwinger": {
3-
"particle_cpu": 1024.0,
4-
"particle_id": 1573888.0,
5-
"particle_momentum_x": 6.569424596732217e-12,
6-
"particle_momentum_y": 3.3588025031324e-27,
7-
"particle_momentum_z": 5.121822174810324e-12,
8-
"particle_position_x": 0.0005120000000000001,
9-
"particle_position_y": 0.0005120000000000001,
10-
"particle_position_z": 0.0005119999999999999,
11-
"particle_weight": 114054.5697167885082308202981948852539062500000
3+
"particle_cpu": 384.0,
4+
"particle_id": 516672.0,
5+
"particle_momentum_x": 0.0000000000036953013356618679404339996551,
6+
"particle_momentum_y": 0.0000000000000000000000000098544358005409,
7+
"particle_momentum_z": 0.0000000000028810249733308377960296374853,
8+
"particle_position_x": 0.0002880000000000000615063555642336723395,
9+
"particle_position_y": 0.0002159999999999999919196580489000325542,
10+
"particle_position_z": 0.0002540000004700523764458730546778042481,
11+
"particle_weight": 63968.7024717521271668374538421630859375000000
1212
},
1313
"lev=0": {
14-
"Bx": 1.4000605828203265268427912815241143107414,
14+
"Bx": 10.2823324587886126835201139329001307487488,
1515
"By": 1707847966720000.0,
16-
"Bz": 0.3684660435860281357811629732168512418866,
17-
"Ex": 1316425589.8291716575622558593750000000000000000000,
16+
"Bz": 0.5093466537801442095556581080018077045679,
17+
"Ex": 3460290775.8437714576721191406250000000000000000000,
1818
"Ey": 0.0,
1919
"Ez": 5.119999999998239e+23,
20-
"jx": 195136835283257.5000000000000000000000000000000000000000,
20+
"jx": 354749214244724.1250000000000000000000000000000000000000,
2121
"jy": 0.0,
22-
"jz": 1.5866778635486594e+16
22+
"jz": 8899049325589567.0000000000000000000000000000000000000000
2323
},
2424
"pos_schwinger": {
25-
"particle_cpu": 1024.0,
26-
"particle_id": 2622464.0,
27-
"particle_momentum_x": 6.569424596732212e-12,
28-
"particle_momentum_y": 3.4619972045013e-27,
29-
"particle_momentum_z": 5.1218221748102914e-12,
30-
"particle_position_x": 0.0005120000000000001,
31-
"particle_position_y": 0.0005120000000000001,
32-
"particle_position_z": 0.000512,
33-
"particle_weight": 114054.5697167885227827355265617370605468750000
25+
"particle_cpu": 384.0,
26+
"particle_id": 959040.0,
27+
"particle_momentum_x": 0.0000000000036953013356618671326404327088,
28+
"particle_momentum_y": 0.0000000000000000000000000103981113667251,
29+
"particle_momentum_z": 0.0000000000028810249733308252752293498174,
30+
"particle_position_x": 0.0002880000000000000072962469399584506391,
31+
"particle_position_y": 0.0002159999999999999919196580489000325542,
32+
"particle_position_z": 0.0002379999995299475479271222866373136640,
33+
"particle_weight": 63968.7024717521271668374538421630859375000000
3434
}
3535
}

Regression/WarpX-tests.ini

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1142,7 +1142,7 @@ tolerance = 1.e-14
11421142
[qed_schwinger2]
11431143
buildDir = .
11441144
inputFile = Examples/Modules/qed/schwinger/inputs_3d_schwinger
1145-
runtime_params = warpx.E_external_grid = 1.e18 0 0 warpx.B_external_grid = 1679288857.0516706 525665014.1557486 1836353079.9561853
1145+
runtime_params = warpx.E_external_grid = 1.e18 0 0 warpx.B_external_grid = 1679288857.0516706 525665014.1557486 1836353079.9561853 qed_schwinger.xmin = -2.5e-7 qed_schwinger.xmax = 2.49e-7
11461146
dim = 3
11471147
addToCompileString = QED=TRUE
11481148
restartTest = 0
@@ -1174,7 +1174,7 @@ tolerance = 1.e-14
11741174
[qed_schwinger4]
11751175
buildDir = .
11761176
inputFile = Examples/Modules/qed/schwinger/inputs_3d_schwinger
1177-
runtime_params = warpx.E_external_grid = 0 0 2.5e+20 warpx.B_external_grid = 0 833910140000. 0
1177+
runtime_params = warpx.E_external_grid = 0 0 2.5e+20 warpx.B_external_grid = 0 833910140000. 0 qed_schwinger.ymin = -2.5e-7 qed_schwinger.zmax = 2.49e-7
11781178
dim = 3
11791179
addToCompileString = QED=TRUE
11801180
restartTest = 0

Source/Particles/MultiParticleContainer.H

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,12 @@ public:
148148
* particles.
149149
*/
150150
void doQEDSchwinger ();
151+
152+
/** This function computes the box outside which Schwinger process is disabled. The box is
153+
* defined by m_qed_schwinger_xmin/xmax/ymin/ymax/zmin/zmax and the warpx level 0 geometry
154+
* object (to make the link between Real and int quatities).
155+
*/
156+
amrex::Box ComputeSchwingerGlobalBox () const;
151157
#endif
152158

153159
void Restart (const std::string& dir);
@@ -391,6 +397,15 @@ protected:
391397
* a Poisson distribution for the pair production rate calculations
392398
*/
393399
int m_qed_schwinger_threshold_poisson_gaussian = 25;
400+
/** The 6 following variables are spatial boundaries beyond which Schwinger process is
401+
* deactivated
402+
*/
403+
amrex::Real m_qed_schwinger_xmin = std::numeric_limits<amrex::Real>::lowest();
404+
amrex::Real m_qed_schwinger_xmax = std::numeric_limits<amrex::Real>::max();
405+
amrex::Real m_qed_schwinger_ymin = std::numeric_limits<amrex::Real>::lowest();
406+
amrex::Real m_qed_schwinger_ymax = std::numeric_limits<amrex::Real>::max();
407+
amrex::Real m_qed_schwinger_zmin = std::numeric_limits<amrex::Real>::lowest();
408+
amrex::Real m_qed_schwinger_zmax = std::numeric_limits<amrex::Real>::max();
394409

395410
#endif
396411

Source/Particles/MultiParticleContainer.cpp

Lines changed: 82 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -277,8 +277,15 @@ MultiParticleContainer::ReadParameters ()
277277
#if (AMREX_SPACEDIM == 2)
278278
ppq.get("y_size",m_qed_schwinger_y_size);
279279
#endif
280-
ppq.query("threshold_poisson_gaussian",
281-
m_qed_schwinger_threshold_poisson_gaussian);
280+
ppq.query("threshold_poisson_gaussian", m_qed_schwinger_threshold_poisson_gaussian);
281+
ppq.query("xmin", m_qed_schwinger_xmin);
282+
ppq.query("xmax", m_qed_schwinger_xmax);
283+
#if (AMREX_SPACEDIM == 3)
284+
ppq.query("ymin", m_qed_schwinger_ymin);
285+
ppq.query("ymax", m_qed_schwinger_ymax);
286+
#endif
287+
ppq.query("zmin", m_qed_schwinger_zmin);
288+
ppq.query("zmax", m_qed_schwinger_zmax);
282289
}
283290
#endif
284291
initialized = true;
@@ -1099,7 +1106,7 @@ MultiParticleContainer::doQEDSchwinger ()
10991106
"ERROR: Schwinger process only implemented for warpx.do_nodal = 1"
11001107
"or algo.field_gathering = momentum-conserving");
11011108

1102-
const int level_0 = 0;
1109+
constexpr int level_0 = 0;
11031110

11041111
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(warpx.maxLevel() == level_0,
11051112
"ERROR: Schwinger process not implemented with mesh refinement");
@@ -1143,7 +1150,15 @@ MultiParticleContainer::doQEDSchwinger ()
11431150
for (MFIter mfi(Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi )
11441151
{
11451152
// Make the box cell centered to avoid creating particles twice on the tile edges
1146-
const Box& box = enclosedCells(mfi.nodaltilebox());
1153+
amrex::Box box = enclosedCells(mfi.nodaltilebox());
1154+
1155+
// Get the box representing global Schwinger boundaries
1156+
const amrex::Box global_schwinger_box = ComputeSchwingerGlobalBox();
1157+
1158+
// If Schwinger process is not activated anywhere in the current box, we move to the next
1159+
// one. Otherwise we use the intersection of current box with global Schwinger box.
1160+
if (!box.intersects(global_schwinger_box)) {continue;}
1161+
box &= global_schwinger_box;
11471162

11481163
const auto& arrEx = Ex[mfi].array();
11491164
const auto& arrEy = Ey[mfi].array();
@@ -1184,6 +1199,69 @@ MultiParticleContainer::doQEDSchwinger ()
11841199
}
11851200
}
11861201

1202+
amrex::Box
1203+
MultiParticleContainer::ComputeSchwingerGlobalBox () const
1204+
{
1205+
auto & warpx = WarpX::GetInstance();
1206+
constexpr int level_0 = 0;
1207+
amrex::Geometry const & geom = warpx.Geom(level_0);
1208+
1209+
#if (AMREX_SPACEDIM == 3)
1210+
const amrex::Array<amrex::Real,3> schwinger_min{m_qed_schwinger_xmin,
1211+
m_qed_schwinger_ymin,
1212+
m_qed_schwinger_zmin};
1213+
const amrex::Array<amrex::Real,3> schwinger_max{m_qed_schwinger_xmax,
1214+
m_qed_schwinger_ymax,
1215+
m_qed_schwinger_zmax};
1216+
#else
1217+
const amrex::Array<amrex::Real,2> schwinger_min{m_qed_schwinger_xmin,
1218+
m_qed_schwinger_zmin};
1219+
const amrex::Array<amrex::Real,2> schwinger_max{m_qed_schwinger_xmax,
1220+
m_qed_schwinger_zmax};
1221+
#endif
1222+
1223+
// Box inside which Schwinger is activated
1224+
amrex::Box schwinger_global_box;
1225+
1226+
for (int dir=0; dir<AMREX_SPACEDIM; dir++)
1227+
{
1228+
// Dealing with these corner cases should ensure that we don't overflow on the integers
1229+
if (schwinger_min[dir] < geom.ProbLo(dir))
1230+
{
1231+
schwinger_global_box.setSmall(dir, std::numeric_limits<int>::lowest());
1232+
}
1233+
else if (schwinger_min[dir] > geom.ProbHi(dir))
1234+
{
1235+
schwinger_global_box.setSmall(dir, std::numeric_limits<int>::max());
1236+
}
1237+
else
1238+
{
1239+
// Schwinger pairs are currently created on the lower nodes of a cell. Using ceil here
1240+
// excludes all cells whose lower node is strictly lower than schwinger_min[dir].
1241+
schwinger_global_box.setSmall(dir, static_cast<int>(std::ceil(
1242+
(schwinger_min[dir] - geom.ProbLo(dir)) / geom.CellSize(dir))));
1243+
}
1244+
1245+
if (schwinger_max[dir] < geom.ProbLo(dir))
1246+
{
1247+
schwinger_global_box.setBig(dir, std::numeric_limits<int>::lowest());
1248+
}
1249+
else if (schwinger_max[dir] > geom.ProbHi(dir))
1250+
{
1251+
schwinger_global_box.setBig(dir, std::numeric_limits<int>::max());
1252+
}
1253+
else
1254+
{
1255+
// Schwinger pairs are currently created on the lower nodes of a cell. Using floor here
1256+
// excludes all cells whose lower node is strictly higher than schwinger_max[dir].
1257+
schwinger_global_box.setBig(dir, static_cast<int>(std::floor(
1258+
(schwinger_max[dir] - geom.ProbLo(dir)) / geom.CellSize(dir))));
1259+
}
1260+
}
1261+
1262+
return schwinger_global_box;
1263+
}
1264+
11871265
void MultiParticleContainer::doQedEvents (int lev,
11881266
const MultiFab& Ex,
11891267
const MultiFab& Ey,

Source/Particles/ParticleCreation/FilterCreateTransformFromFAB.H

Lines changed: 8 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -58,19 +58,14 @@ Index filterCreateTransformFromFAB (DstTile& dst1, DstTile& dst2, const amrex::B
5858
const int level_0 = 0;
5959
Geometry const & geom = warpx.Geom(level_0);
6060

61-
// Note that we use here x and y for the 2 directions in 2D, which is not consistent with
62-
// the rest of the code (x and z are usually chosen)
61+
constexpr int spacedim = AMREX_SPACEDIM;
62+
6363
const Real xlo_global = geom.ProbLo(0);
64-
const Real ylo_global = geom.ProbLo(1);
6564
const Real dx = geom.CellSize(0);
66-
const Real dy = geom.CellSize(1);
67-
#if (AMREX_SPACEDIM == 3)
68-
const Real zlo_global = geom.ProbLo(2);
69-
const Real dz = geom.CellSize(2);
70-
#else
71-
const Real zlo_global = 0.0;
72-
const Real dz = 0.0;
73-
#endif
65+
const Real ylo_global = (spacedim == 3) ? geom.ProbLo(1) : amrex::Real(0.);
66+
const Real dy = (spacedim == 3) ? geom.CellSize(1) : amrex::Real(0.);
67+
const Real zlo_global = (spacedim == 3) ? geom.ProbLo(2) : geom.ProbLo(1);
68+
const Real dz = (spacedim == 3) ? geom.CellSize(2) : geom.CellSize(1);
7469

7570
const auto arrNumPartCreation = src_FAB->array();
7671
Gpu::DeviceVector<Index> offsets(ncells);
@@ -84,8 +79,6 @@ Index filterCreateTransformFromFAB (DstTile& dst1, DstTile& dst2, const amrex::B
8479
const auto dst1_data = dst1.getParticleTileData();
8580
const auto dst2_data = dst2.getParticleTileData();
8681

87-
constexpr int spacedim = AMREX_SPACEDIM;
88-
8982
// For loop over all cells in the box. If mask is true in the given cell,
9083
// we create the particles in the cell and apply a transform function to the
9184
// created particles.
@@ -102,8 +95,8 @@ Index filterCreateTransformFromFAB (DstTile& dst1, DstTile& dst2, const amrex::B
10295
// Currently all particles are created on nodes. This makes it useless
10396
// to use N>1 (for now).
10497
const Real x = xlo_global + j*dx;
105-
const Real y = (spacedim == 3) ? ylo_global + k*dy : 0.0_rt;
106-
const Real z = (spacedim == 3) ? zlo_global + l*dz : ylo_global + k*dy;
98+
const Real y = ylo_global + k*dy;
99+
const Real z = (spacedim == 3) ? zlo_global + l*dz : zlo_global + k*dz;
107100

108101
for (int n = 0; n < N; ++n)
109102
{

0 commit comments

Comments
 (0)