Skip to content

Commit 1231aed

Browse files
RemiLeheRemiLehe
authored
Modify stair-case approximation to the EB (#5534)
# Overview This PR changes the definition of the stair-case approximation of the EB (used in the Yee solver), so that the actual EB boundary (where e.g. particles are removed from the simulation) is **inside** the stair-case-approximated boundary - as opposed to the definition used in the current `development` branch, for which the actual EB is **outside** the stair-case-approximated boundary. This ensures that the algorithm remains charge conserving, when charged particles are absorbed or emitted by the embedded boundary, **for particle shape of order 1**. (Higher-order particle shapes will be addressed in #5209.) As illustrated in the figure below (and as discussed in [this paper](https://www.researchgate.net/publication/318642364_Charge_Conserving_Emission_from_Conformal_Boundaries_in_Electromagnetic_PIC_simulations)), this is fundamentally because the particle does not deposit any charge in the valid cells, at the time when it is removed/emitted. <img width="608" alt="Screenshot 2025-01-12 at 8 56 17 AM" src="https://github.com/user-attachments/assets/bd803ef0-faf8-4ea5-973b-240c15b2ba4b" /> (The black crosses show the locations where the electric field is not updated, and thus usually remains equal to 0. The red dots show the locations where the particle deposits charge, for particle shape of order 1.) The better behavior with respect to charge-conservation can be observed in the following animations, which show two particles of opposing charge separating and going into the embedded boundary. (Note that a static error in `divE` remains at the position where the particle was absorbed, with the `development` branch. The propagating errors in `divE` are expected: they are due to electromagnetic waves reflecting on the EB.) - **development branch** ![movie](https://github.com/user-attachments/assets/d486663d-a182-4751-b1d1-709b1a74ea44) - **this PR** ![movie](https://github.com/user-attachments/assets/94a5dea3-2bb8-4548-b320-7615cac86fe7) Input script: [inputs.txt](https://github.com/user-attachments/files/18428873/inputs.txt) Analysis script: [openPMD-visualization.ipynb.txt](https://github.com/user-attachments/files/18428878/openPMD-visualization.ipynb.txt) (An automated tests using a similar configuration has been added in a separate, follow-up PR: #5562) Note that, as part of this PR, the above new definition has been adopted for all the finite-difference solvers, except for the ECT solver (which uses a cut-cell representation instead of a stair-case representation). # Implementation This PR uses the changes of #5574. It still uses `MarkUpdateECellsECT` and `MarkUpdateBCellsECT` for the ECT sover - which preserve the previous behavior of the embedded boundary for this solver, but now uses `MarkUpdateCellsStairCase` for the other FDTD solvers - which introduce the above-mentioned new stair-case definition.
1 parent d7b195f commit 1231aed

14 files changed

+75
-67
lines changed

‎Examples/Tests/embedded_boundary_cube/inputs_base_3d

+2-2
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,8 @@ boundary.field_lo = pec pec pec
1313
boundary.field_hi = pec pec pec
1414

1515
eb2.geom_type = box
16-
eb2.box_lo = -0.5 -0.5 -0.5 # Ensures that the stair-case EB is exactly at -0.5
17-
eb2.box_hi = 0.5 0.5 0.5 # Ensures that the stair-case EB is exactly at 0.5
16+
eb2.box_lo = -0.501 -0.501 -0.501 # Ensures that the stair-case EB is exactly at -0.5
17+
eb2.box_hi = 0.501 0.501 0.501 # Ensures that the stair-case EB is exactly at 0.5
1818
eb2.box_has_fluid_inside = true
1919
# Alternatively one could use parser to build EB
2020
# Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular.

‎Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube

+4-4
Original file line numberDiff line numberDiff line change
@@ -12,10 +12,10 @@ warpx.abort_on_warning_threshold = medium
1212
boundary.field_lo = pec pec
1313
boundary.field_hi = pec pec
1414

15-
my_constants.xmin = -0.5
16-
my_constants.zmin = -0.5
17-
my_constants.xmax = 0.5
18-
my_constants.zmax = 0.5
15+
my_constants.xmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5
16+
my_constants.zmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5
17+
my_constants.xmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5
18+
my_constants.zmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5
1919
# Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular.
2020
warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))"
2121

‎Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json

+2-2
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,8 @@
33
"Bx": 9.263694545408503e-05,
44
"By": 0.00031905198933489145,
55
"Bz": 7.328424783762594e-05,
6-
"Ex": 8553.906698053046,
6+
"Ex": 8553.90669811286,
77
"Ey": 60867.04830538045,
8-
"Ez": 8.439422682267567e-07
8+
"Ez": 4.223902107031194e-06
99
}
1010
}

‎Regression/Checksum/benchmarks_json/test_2d_field_probe.json

+4-4
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
33
"Bx": 0.0,
4-
"By": 126826.78487921853,
4+
"By": 123510.69657444415,
55
"Bz": 0.0,
6-
"Ex": 32517064310550.266,
6+
"Ex": 31206368949280.34,
77
"Ey": 0.0,
8-
"Ez": 17321323003697.61
8+
"Ez": 16921005306450.537
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json

+7-7
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Bx": 148673.005859208,
4-
"By": 148673.00585920806,
5-
"Bz": 3371.758117878558,
6-
"Ex": 55378581103426.71,
7-
"Ey": 55378581103426.72,
8-
"Ez": 68412803445328.25
3+
"Bx": 144495.08082507108,
4+
"By": 144495.08082507114,
5+
"Bz": 8481.958724628861,
6+
"Ex": 54500496182517.92,
7+
"Ey": 54500496182517.91,
8+
"Ez": 70231240245509.39
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json

+5-5
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Bx": 4.060477854092961e-18,
3+
"Bx": 4.166971025838921e-18,
44
"By": 0.006628374119786834,
55
"Bz": 0.006628374119786834,
6-
"Ex": 5102618.4711524295,
7-
"Ey": 6.323754160591239e-05,
8-
"Ez": 6.323754160591239e-05
6+
"Ex": 5102618.471153786,
7+
"Ey": 1.4283859321773714e-05,
8+
"Ez": 1.4283859321773714e-05
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json

+5-5
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Bx": 4.20930075273562e-18,
3+
"Bx": 4.228863291892693e-18,
44
"By": 0.005101824310293573,
55
"Bz": 0.005101824310293573,
6-
"Ex": 4414725.184731115,
7-
"Ey": 6.32375413967707e-05,
8-
"Ez": 6.32375413967707e-05
6+
"Ex": 4414725.184732471,
7+
"Ey": 1.4283895626502055e-05,
8+
"Ez": 1.4283895626502055e-05
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json

+1-2
Original file line numberDiff line numberDiff line change
@@ -7,5 +7,4 @@
77
"Ey": 1.042254197269831e+04,
88
"Ez": 1.040011664019071e+04
99
}
10-
}
11-
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json

+7-7
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Bx": 202106.71291347666,
4-
"By": 202106.71291347663,
5-
"Bz": 3371.897999274175,
6-
"Ex": 38304043178806.11,
7-
"Ey": 38304043178806.11,
8-
"Ez": 83057027925874.84
3+
"Bx": 198610.0530604908,
4+
"By": 198610.0530604909,
5+
"Bz": 8482.656173586969,
6+
"Ex": 37232105734622.53,
7+
"Ey": 37232105734622.54,
8+
"Ez": 85094015810307.19
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json

+7-7
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Bx": 148673.005859208,
4-
"By": 148673.00585920803,
5-
"Bz": 3371.758117878557,
6-
"Ex": 55378581103426.695,
7-
"Ey": 55378581103426.7,
8-
"Ez": 68412803445328.25
3+
"Bx": 144495.08082507108,
4+
"By": 144495.0808250711,
5+
"Bz": 8481.95872462886,
6+
"Ex": 54500496182517.914,
7+
"Ey": 54500496182517.914,
8+
"Ez": 70231240245509.4
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json

+7-7
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Bx": 148673.005859208,
4-
"By": 148673.00585920803,
5-
"Bz": 3371.758117878557,
6-
"Ex": 55378581103426.695,
7-
"Ey": 55378581103426.7,
8-
"Ez": 68412803445328.25
3+
"Bx": 144495.08082507108,
4+
"By": 144495.08082507114,
5+
"Bz": 8481.958724628861,
6+
"Ex": 54500496182517.92,
7+
"Ey": 54500496182517.91,
8+
"Ez": 70231240245509.39
99
}
10-
}
10+
}

‎Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json

+6-6
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
{
22
"lev=0": {
3-
"Br": 6.821267675779345e-19,
4-
"Bt": 5.564905732478707e-05,
5-
"Bz": 2.368259586613272e-19,
6-
"Er": 16503.98082446463,
7-
"Et": 1.5299584682447838e-10,
8-
"Ez": 1466.854467399728
3+
"Br": 6.7914286131989935e-19,
4+
"Bt": 5.4557350206853276e-05,
5+
"Bz": 2.357229221622199e-19,
6+
"Er": 16481.39008058988,
7+
"Et": 1.5258937379236053e-10,
8+
"Ez": 1508.1064116028576
99
}
1010
}

‎Source/EmbeddedBoundary/WarpXInitEB.cpp

+3-4
Original file line numberDiff line numberDiff line change
@@ -320,14 +320,14 @@ WarpX::MarkUpdateCellsStairCase (
320320

321321
if (fab_type == amrex::FabType::regular) { // All cells in the box are regular
322322

323-
// Every cell in box is all regular: update field in every cell
323+
// Every cell in box is regular: update field in every cell
324324
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
325325
eb_update_arr(i, j, k) = 1;
326326
});
327327

328328
} else if (fab_type == amrex::FabType::covered) { // All cells in the box are covered
329329

330-
// Every cell in box is all covered: do not update field
330+
// Every cell in box is fully covered: do not update field
331331
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
332332
eb_update_arr(i, j, k) = 0;
333333
});
@@ -369,8 +369,7 @@ WarpX::MarkUpdateCellsStairCase (
369369
for (int k_cell = k_start; k_cell <= k; ++k_cell) {
370370
// If one of the neighboring is either partially or fully covered
371371
// (i.e. if they are not regular cells), do not update field
372-
// (Note that `flag` is a cell-centered object, and `isRegular`
373-
// returns `false` if the cell is either partially or fully covered.)
372+
// (`isRegular` returns `false` if the cell is either partially or fully covered.)
374373
if ( !flag(i_cell, j_cell, k_cell).isRegular() ) {
375374
eb_update = 0;
376375
}

‎Source/Initialization/WarpXInitData.cpp

+15-5
Original file line numberDiff line numberDiff line change
@@ -1253,13 +1253,23 @@ void WarpX::InitializeEBGridData (int lev)
12531253
// Compute additional quantities required for the ECT solver
12541254
MarkExtensionCells();
12551255
ComputeFaceExtensions();
1256+
// Mark on which grid points E should be updated
1257+
MarkUpdateECellsECT( m_eb_update_E[lev], edge_lengths_lev );
1258+
// Mark on which grid points B should be updated
1259+
MarkUpdateBCellsECT( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev);
1260+
} else {
1261+
// Mark on which grid points E should be updated (stair-case approximation)
1262+
MarkUpdateCellsStairCase(
1263+
m_eb_update_E[lev],
1264+
m_fields.get_alldirs(FieldType::Efield_fp, lev),
1265+
eb_fact );
1266+
// Mark on which grid points B should be updated (stair-case approximation)
1267+
MarkUpdateCellsStairCase(
1268+
m_eb_update_B[lev],
1269+
m_fields.get_alldirs(FieldType::Bfield_fp, lev),
1270+
eb_fact );
12561271
}
12571272

1258-
// Mark on which grid points E should be updated
1259-
MarkUpdateECellsECT( m_eb_update_E[lev], edge_lengths_lev );
1260-
// Mark on which grid points B should be updated
1261-
MarkUpdateBCellsECT( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev);
1262-
12631273
}
12641274

12651275
ComputeDistanceToEB();

0 commit comments

Comments
 (0)