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

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.

…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/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.

Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json

Source/EmbeddedBoundary/WarpXInitEB.cpp

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.

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 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 functions which keep the same EB definition. (PR #5534 will then change this 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. --------- Co-authored-by: Roelof Groenewald <[email protected]>

Source/Initialization/WarpXInitData.cpp

roelof-groenewald

Looks good. Thanks!

Merge #5534 first. This extends the changes of #5534 and #5574 to the hybrid solver. Note that this effectively changes the definition of the stair-cased embedded boundary for the hybrid solver. --------- Co-authored-by: Roelof Groenewald <[email protected]> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

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

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

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

New embedded boundary data structures in preparation for #5534 #5574

Merged

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

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

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

Merge branch 'development' into move_stair_case_approx

9b8282a

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

Move the boundary

ac82004

RemiLehe disabled auto-merge January 21, 2025 17:35

roelof-groenewald reviewed Jan 21, 2025

View reviewed changes

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

Correct code

6147b95

roelof-groenewald approved these changes Jan 21, 2025

View reviewed changes

roelof-groenewald merged commit 1231aed into ECP-WarpX:development Jan 21, 2025
37 checks passed

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 commented Jan 19, 2025

roelof-groenewald left a comment

Modify stair-case approximation to the EB #5534

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

RemiLehe commented Jan 19, 2025

roelof-groenewald left a comment

Choose a reason for hiding this comment

RemiLehe commented Jan 3, 2025 •

edited

Loading