Modify stair-case approximation to the EB #5534

RemiLehe · 2025-01-03T00:57:31Z

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), this is fundamentally because the particle does not deposit any charge in the valid cells, at the time when it is removed/emitted.

(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
this PR

Input script: inputs.txt
Analysis script: 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

When updating E (and when initializing the field values for E and B with a parser), we used to check somewhat complicated conditions that relied on the MultiFabs edge_lengths and face_areas.

This PR now introduces separate arrays (m_eb_update_E and m_eb_update_B, which are iMultiFab instead of MultiFab and thus take up much less memory). These arrays contain flags that keep track of where to udpate the fields (i.e. the black crosses in the above figure). These arrays are initialized in separate function (MarkUpdateECellsECT and MarkUpdateBCellsECT for the ECT sover - which preserve the previous behavior of the embedded boundary for this solver ; MarkUpdateCellsStairCase for the other FDTD solvers - which introduce the above-mentioned new stair-case definition). The code for the field pusher and field initialization uses these arrays (instead of edge_lengths and face_areas). It is thus significantly simpler and avoids duplicating complicated if conditions.

The above changes have not yet been implemented for the hybrid solver. This will instead be done in a separate PR: #5558. Once this is complete, the MultiFab edge_lengths and face_areas will not be needed anymore (except for the ECT solver), and we will thus only allocate them for the ECT solver. This should save a significant amount of memory.

…to move_stair_case_approx

This reverts commit 4fd4dc7.

RemiLehe · 2025-01-12T16:31:05Z

Source/Parallelization/WarpXRegrid.cpp

@@ -174,6 +174,14 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi
    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

+    const auto RemakeMultiFab = [&](auto& mf){


This reintroduces a function that we removed when introducing the MultiFabRegister (#5230).

It could be removed once we introduce an iMultiFabRegister.

RemiLehe · 2025-01-12T16:53:32Z

Source/EmbeddedBoundary/WarpXInitEB.cpp


 void
-WarpX::MarkCells ()
+WarpX::MarkExtensionCells ()


This function was renamed to avoid confusion with the other Mark... function.

RemiLehe · 2025-01-12T21:57:23Z

Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp

-#endif
+
+                // Skip field push in the embedded boundaries
+                if (update_Ey_arr && update_Ey_arr(i, j, k) == 0) { return; }


The if condition here is now more simple, since the logic that determines whether a cell should be updated has been moved to the function that initializes the update_E flag.

Source/Initialization/WarpXInitData.cpp

Source/WarpX.cpp

roelof-groenewald

This looks good to me. I was wondering, though, did you by any chance try running the CI tests using the new m_eb_update_E functionality but before moving the staircase approximation? That should be a really clean way to check that there are no hard to see bugs introduced in this PR, since all the CI tests should then pass with their original checksum values.

roelof-groenewald · 2025-01-17T00:28:54Z

Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json

-    "Ez": 68412803445328.25
+    "Bx": 144495.08082507108,
+    "By": 144495.08082507114,
+    "Bz": 8481.958724628861,


This change is quite big (more than 2x). Do you understand why this component had such a dramatic change?
I guess the Bz value is small compared to the others so this is probably fine.

Hm, yes, it is somewhat large. My guess is that this is because this PR changes the actual position of the EB for this test.

roelof-groenewald · 2025-01-17T01:09:02Z

Source/EmbeddedBoundary/WarpXInitEB.cpp

+                    // Stair-case approximation: If neighboring cells are either partially
+                    // or fully covered: do not update field


If I'm understanding the logic here correctly, eb_flag is defined on the cell centers and returns false for isRegular() if any part of the cell is touched by the EB. The reason you have to check the cell below the current one when dealing with a nodal MF is because for the same index, nodal MFs are offset by half a cell downward from the cell centered MF. I think it might be worth adding a bit more explanation in the comment to make it easier for future readers to understand why it is only necessary to check the index offset by -1 (but never +1).

Yes, that's correct. I added more comments ; could you check if that's clearer now?

roelof-groenewald · 2025-01-17T01:14:27Z

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp

+            lev, PatchType::fine,
+            warpx.GetEBUpdateEFlag());


This is a slight problem in that there will be an inconsistency in the hybrid solver until the follow-up PR is merged that also updates the field update for the hybrid solver. I guess we should just be sure to have both PRs merge in quick succession.

Sounds good.
I finalized the hybrid PR #5558 and the tests pass on that PR, so it is ready to be merged soon after this one.

roelof-groenewald · 2025-01-17T01:18:23Z

Source/FieldSolver/WarpXPushFieldsEM.cpp

@@ -1155,7 +1157,7 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt, amre
        m_fields.get_alldirs(FieldType::Efield_fp, lev),
        m_fields.get_alldirs(FieldType::Bfield_fp, lev),
        m_fields.get_alldirs(FieldType::current_fp, lev),
-        m_fields.get_alldirs(FieldType::edge_lengths, lev),
+        m_eb_update_E[lev],


I'm assuming the plan is to replace these arguments once iMultiFabRegister is introduced?

Yes, I agree.

roelof-groenewald · 2025-01-17T01:19:12Z

Source/Initialization/WarpXInitData.cpp

@@ -976,23 +976,21 @@ WarpX::InitLevelData (int lev, Real /*time*/)
    // The default maxlevel_extEMfield_init value is the total number of levels in the simulation
    if ((m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function)
         && (lev > 0) && (lev <= maxlevel_extEMfield_init)) {
+
+        // TODO: raise error when EB is activated


Why is this comment here? Why should an error be raised when EB is activated?

This is because I realized that the code is probably incorrect when using MR with the EM solver (i.e. when we use the cp version of the fields).

This is because we should also have cp and fp versions of the arrays eb_update_E / eb_update_B.

That being said, the problem already existed before this PR. (i.e. we should have had cp and fp version of edge_length and face_area).

For now, I cleaned up this comment and will instead open an issue with more explanation.

roelof-groenewald · 2025-01-17T01:20:26Z

Source/Initialization/WarpXInitData.cpp

-#ifdef AMREX_USE_EB
-#ifdef WARPX_DIM_3D
-                if(lx && ((topology=='e' and lx(i, j, k)<=0) or (topology=='f' and Sx(i, j, k)<=0))) { return; }
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-                //In XZ and RZ Ex is associated with a x-edge, while Bx is associated with a z-edge
-                if(lx && ((topology=='e' and lx(i, j, k)<=0) or (topology=='f' and lz(i, j, k)<=0))) { return; }
-#endif
-#endif


This is a great simplification! 🎉

RemiLehe · 2025-01-19T15:11:00Z

This looks good to me. I was wondering, though, did you by any chance try running the CI tests using the new m_eb_update_E functionality but before moving the staircase approximation? That should be a really clean way to check that there are no hard to see bugs introduced in this PR, since all the CI tests should then pass with their original checksum values.

That's a good point! I have done this in this PR: #5574
That PR is based on #5534 and uses the new eb_update_E/eb_update_B framework, but:

it initializes them with the previous definition of the EB, by using MarkUpdateECellsECT/MarkUpdateBCellsECT for all solvers (not just the ECT solver):
3a3f36e
it does not modify the checksum (compared to the current development)

The checksum tests indeed pass for #5574.

Modify stair-case approximation

3fa996f

RemiLehe changed the title ~~[WIP] Modify stair-case approximation~~ [WIP] Modify stair-case approximation to the EB Jan 3, 2025

RemiLehe added 4 commits January 5, 2025 08:10

Fix compilation without EB

112cce1

Slightly modify boundaries

b0ccce4

Added new arrays with flags

76fa82e

Finalize initialization of the array

3624a72

RemiLehe force-pushed the move_stair_case_approx branch from 3a60ef0 to 3181d45 Compare January 8, 2025 19:38

RemiLehe added 19 commits January 8, 2025 11:39

Use guard cell information

3181d45

Remove edge_length in field update

38d967d

Fix const-ness

79ed114

Fix compilation warning

8516d1a

Updated init functions

3c0867e

Fix initialization criterion

28d43a1

Add flag for B

0214769

Fix RZ compilation error

ab87809

Update initialization of external fields

78bb8c1

Fix setting of eb_update

df9ae82

Implement ECT condition for B

59579b4

Do not use any guard cells

4fd4dc7

Fixed ECT solver

4232db7

Merge branch 'move_stair_case_approx' of github.com:RemiLehe/WarpX in…

4b8079f

…to move_stair_case_approx

Revert "Do not use any guard cells"

f45b257

This reverts commit 4fd4dc7.

Fix compilation errors

2003d85

Fix compilation error

d6e64af

Fix const-ness

e6d93f9

Activate load-balancing

eb43087

RemiLehe commented Jan 12, 2025

View reviewed changes

Update checksums

2eb3383

RemiLehe commented Jan 12, 2025

View reviewed changes

Add documentation

8981179

RemiLehe mentioned this pull request Jan 14, 2025

[WIP] Use new stair-case approximation in hybrid solver #5558

Open

RemiLehe commented Jan 14, 2025

View reviewed changes

Source/Initialization/WarpXInitData.cpp Outdated Show resolved Hide resolved

RemiLehe commented Jan 14, 2025

View reviewed changes

Source/WarpX.cpp Outdated Show resolved Hide resolved

Apply suggestions from code review

8182df7

RemiLehe requested review from roelof-groenewald and dpgrote January 14, 2025 20:54

RemiLehe assigned roelof-groenewald Jan 14, 2025

RemiLehe changed the title ~~[WIP] Modify stair-case approximation to the EB~~ Modify stair-case approximation to the EB Jan 14, 2025

RemiLehe added component: boundary PML, embedded boundaries, et al. component: FDTD FDTD solvers component: initialization Changes related to the initialization of the simulation labels Jan 14, 2025

RemiLehe mentioned this pull request Jan 15, 2025

[WIP] Add automated test for particle absorption on new stair-case approximation #5562

Open

roelof-groenewald reviewed Jan 17, 2025

View reviewed changes

RemiLehe added 2 commits January 18, 2025 08:34

Minor cleanup

0980516

Modify initialization in ECT solver

ace399b

RemiLehe force-pushed the move_stair_case_approx branch from bc10261 to 16ec9b2 Compare January 18, 2025 18:20

Update comments

16ec9b2

RemiLehe mentioned this pull request Jan 19, 2025

[Do not merge] Test checksums with new EB structures #5574

Open

RemiLehe requested a review from roelof-groenewald January 19, 2025 15:18

RemiLehe enabled auto-merge (squash) January 20, 2025 23:32

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Modify stair-case approximation to the EB #5534

Modify stair-case approximation to the EB #5534

RemiLehe commented Jan 3, 2025 •

edited

Loading

RemiLehe Jan 12, 2025

RemiLehe Jan 12, 2025

RemiLehe Jan 12, 2025

roelof-groenewald left a comment

roelof-groenewald Jan 17, 2025

RemiLehe Jan 19, 2025

roelof-groenewald Jan 17, 2025

RemiLehe Jan 19, 2025

roelof-groenewald Jan 17, 2025

RemiLehe Jan 19, 2025

roelof-groenewald Jan 17, 2025

RemiLehe Jan 19, 2025

roelof-groenewald Jan 17, 2025

RemiLehe Jan 19, 2025

roelof-groenewald Jan 17, 2025

RemiLehe commented Jan 19, 2025

		// Stair-case approximation: If neighboring cells are either partially
		// or fully covered: do not update field

Modify stair-case approximation to the EB #5534

Are you sure you want to change the base?

Modify stair-case approximation to the EB #5534

Conversation

RemiLehe commented Jan 3, 2025 • edited Loading

Overview

Implementation

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

roelof-groenewald left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

RemiLehe commented Jan 19, 2025

RemiLehe commented Jan 3, 2025 •

edited

Loading