Fix bug with ES solver and MR: E_aux=E_fp in UpdateAuxilaryData (#…

…4922)

* Removal of asserttion which prevented from usung the averaged PSATD algorithms with PML BC

* Clean-up

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Fixed to arr_aux(j,k,l) = fine when ES solve is used

* Removed temporary print statements

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Clean-up

* Added CI test ElectrostaticSphereEB_RZ_MR_lev_1 to check the fields on the level=1

* Updated becnmarks for ElectrostaticSphereLabFrame_MR_emass_10

* Imported regular expression (re) in the analysis script.

* Fixed typo

* United two CI tests for different levels of MR in one test.

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Updated CI test ElectrostaticSphereEB_RZ_MR and the corresponding analysis script with  smaller MR patch.

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Clean-up

* Do deepcopy for lev>0 and collocated grid

* Fix bugs to resolve failure of CI tests

* Preserve plotfile output, update benchmark file

* Working on  CI test

* Update benchmark of `ElectrostaticSphereEB_RZ_MR`

* Initialize with value all `aux`, `cax` fields

* Remove changes related to averaged Galilean PSATD with PML

* Remove style changes (e.g., changes to empty lines)

* Replace `MFIter`/`ParallelFor` loop with simple copy

* Apply suggestions from code review

* Revert part of the code to its previous, equivalent state

* Removed DeepCopy & no need for ghost cells in temp phi_cp.

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update Source/ablastr/fields/PoissonSolver.H

* Add inline comments

* Update analysis script

* Use `TilingIfNotGPU`, `growntilebox` in E loop

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <ezoni@lbl.gov>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+
+# Copyright 2024 Olga Shapoval, Edoardo Zoni
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# This script tests the embedded boundary in RZ.
+# A cylindrical surface (r=0.1) has a fixed potential 1 V.
+# The outer surface has 0 V fixed.
+# Thus the analytical solution has the form:
+# phi(r) = A+B*log(r), Er(r) = -B/r.
+
+import os
+import sys
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+tolerance = 0.004
+print(f'tolerance = {tolerance}')
+
+fn = sys.argv[1]
+
+def find_first_non_zero_from_bottom_left(matrix):
+    for i in range(matrix.shape[0]):
+        for j in range(matrix.shape[1]):
+            if (matrix[i][j] != 0) and (matrix[i][j] != np.nan):
+                return (i, j)
+    return i, j
+
+def find_first_non_zero_from_upper_right(matrix):
+    for i in range(matrix.shape[0]-1, -1, -1):
+        for j in range(matrix.shape[1]-1, -1, -1):
+            if (matrix[i][j] != 0) and (matrix[i][j] != np.nan):
+                return (i, j)
+    return i,j
+
+def get_fields(ts, level):
+    if level == 0:
+        Er, info = ts.get_field('E', 'r', iteration=0)
+        phi, info = ts.get_field('phi', iteration=0)
+    else:
+        Er, info = ts.get_field(f'E_lvl{level}', 'r', iteration=0)
+        phi, info = ts.get_field(f'phi_lvl{level}', iteration=0)
+    return Er, phi, info
+
+def get_error_per_lev(ts,level):
+    Er, phi, info = get_fields(ts, level)
+
+    nr_half = info.r.shape[0] // 2
+    dr = info.dr
+
+    Er_patch = Er[:,nr_half:]
+    phi_patch = phi[:,nr_half:]
+    r1 = info.r[nr_half:]
+    patch_left_lower_i, patch_left_lower_j = find_first_non_zero_from_bottom_left(Er_patch)
+    patch_right_upper_i, patch_right_upper_j = find_first_non_zero_from_upper_right(Er_patch)
+
+    # phi and Er field on the MR patch
+    phi_sim = phi_patch[patch_left_lower_i:patch_right_upper_i+1, patch_left_lower_j:patch_right_upper_j+1]
+    Er_sim =  Er_patch[patch_left_lower_i:patch_right_upper_i+1, patch_left_lower_j:patch_right_upper_j+1]
+    r =  r1[patch_left_lower_j:patch_right_upper_j+1]
+
+    B = 1.0/np.log(0.1/0.5)
+    A = -B*np.log(0.5)
+
+    # outside EB and last cutcell
+    rmin = np.min(np.argwhere(r >= (0.1+dr)))
+    rmax = -1
+    r = r[rmin:rmax]
+    phi_sim = phi_sim[:,rmin:rmax]
+    Er_sim = Er_sim[:,rmin:rmax]
+
+    phi_theory = A + B*np.log(r)
+    phi_theory = np.tile(phi_theory, (phi_sim.shape[0],1))
+    phi_error = np.max(np.abs(phi_theory-phi_sim) / np.abs(phi_theory))
+
+    Er_theory = -B/r
+    Er_theory = np.tile(Er_theory, (Er_sim.shape[0],1))
+    Er_error = np.max(np.abs(Er_theory-Er_sim) / np.abs(Er_theory))
+
+    print(f'max error of phi[lev={level}]: {phi_error}')
+    print(f'max error of  Er[lev={level}]: {Er_error}')
+    assert(phi_error < tolerance)
+    assert(Er_error < tolerance)
+
+ts = OpenPMDTimeSeries(fn)
+level_fields = [field for field in ts.avail_fields if 'lvl' in field]
+nlevels = 0 if level_fields == [] else int(level_fields[-1][-1])
+for level in range(nlevels+1):
+    get_error_per_lev(ts,level)
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn, output_format="openpmd")

Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr

@@ -30,3 +30,4 @@ diagnostics.diags_names = diag1
 diag1.intervals = 1
 diag1.diag_type = Full
 diag1.fields_to_plot = Er phi
+diag1.format = openpmd

Regression/Checksum/benchmarks_json/ElectrostaticSphereEB_RZ_MR.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Er": 8487.661571739109,
-    "phi": 2036.0428085225362
+    "Er": 16975.32314347822,
+    "phi": 4072.085617045073
   },
   "lev=1": {
-    "Er": 19519.172334977942,
-    "phi": 3291.0262856782897
+    "Er": 26818.189739547757,
+    "phi_lvl1": 8731.176548788893
   }
 }

Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json

@@ -6,15 +6,15 @@
     "rho": 0.0
   },
   "lev=1": {
-    "Ex": 14.281015560380963,
-    "Ey": 14.281015560380965,
-    "Ez": 14.281015560380965,
+    "Ex": 7.170105936287823,
+    "Ey": 7.17010593628782,
+    "Ez": 7.170105936287821,
     "rho": 2.6092568008333786e-10
   },
   "electron": {
-    "particle_momentum_x": 1.80842228672388e-24,
-    "particle_momentum_y": 1.8084222867238806e-24,
-    "particle_momentum_z": 1.7598771525647628e-24,
+    "particle_momentum_x": 9.257577597262615e-25,
+    "particle_momentum_y": 9.257577597262618e-25,
+    "particle_momentum_z": 9.257577597262624e-25,
     "particle_position_x": 327.46875,
     "particle_position_y": 327.46875,
     "particle_position_z": 327.46875,

Regression/WarpX-tests.ini

@@ -483,7 +483,7 @@ analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz.py
 [ElectrostaticSphereEB_RZ_MR]
 buildDir = .
 inputFile = Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr
-runtime_params =  warpx.abort_on_warning_threshold = medium
+runtime_params =  warpx.abort_on_warning_threshold = medium amr.ref_ratio_vect = 2 2 2
 dim = 2
 addToCompileString = USE_EB=TRUE USE_RZ=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_EB=ON
@@ -492,7 +492,8 @@ useMPI = 1
 numprocs = 2
 useOMP = 1
 numthreads = 1
-analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz.py
+outputFile = ElectrostaticSphereEB_RZ_MR_plt
+analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py
 
 [ElectrostaticSphereLabFrame]
 buildDir = .