Implement stair-case Yee solver with EB in RZ geometry (ECP-WarpX#2707)

* Allow compilation with RZ EB * Do not push cells for RZ Yee solver, when covered with EB * Fix compilation errors * Fix additional compilation errors * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix additional compilation errors * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Add automated test * Add automated test * Fix path in tests * Enable parser in RZ * Update example script * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Clean-up PR * Initialize EB quantities * Modified EM field initialization in 2D with EB * Typo fix * Typo fix * Ignoring unused variables correctly * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Correct condition for updating E * Correct update of B * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update B push * Update input script * Revert "Update input script" This reverts commit 5087485. * Update initialization * Updated test * Move test to a different folder * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Add test for WarpX-test.ini * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix path for tests * Update test description * Update test metadata * Add checksum file * Revert changes * Revert changes * Change lx to lr * Revert "Change lx to lr" This reverts commit be3039a. * Change lx to lr --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: lgiacome <[email protected]>
clarkse · Mar 9, 2024 · e38acc4 · e38acc4
1 parent ebbe634
commit e38acc4
Show file tree

Hide file tree

Showing 9 changed files with 144 additions and 22 deletions.
diff --git a/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py b/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+"""
+This test checks the implementation of the embedded boundary in cylindrical geometry,
+by checking the diffraction of a laser by an embedded boundary here.
+We then check that the first minimum of the diffracted intensity pattern
+occurs along the angle given by the theoretical Airy pattern, i.e.
+theta_diffraction = 1.22 * lambda / d
+"""
+import os
+import sys
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.ndimage import gaussian_filter1d
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+ts = OpenPMDTimeSeries('./EmbeddedBoundaryDiffraction_plt/')
+
+# Extract the intensity as a function of r and z
+Ex, info = ts.get_field('E', 'x', iteration=300)
+I = gaussian_filter1d(Ex**2, sigma=5, axis=0) # Extract intensity by averaging E^2 over wavelength
+irmax = np.argmax( I, axis=-1)
+
+# Find the radius of the first minimum, as a function of z
+def r_first_minimum(iz):
+    ir = len(info.r)//2
+    while I[iz, ir+1] < I[iz, ir]:
+        ir += 1
+    return info.r[ir]
+r = np.array([ r_first_minimum(iz) for iz in range(len(info.z)) ])
+
+# Check that this corresponds to the prediction from the Airy pattern
+theta_diffraction = np.arcsin(1.22*0.1/0.4)/2
+assert np.all( abs(r[50:] - theta_diffraction*info.z[50:]) < 0.03 )
+
+# Open the right plot file
+filename = sys.argv[1]
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename, output_format='openpmd')
diff --git a/Examples/Tests/embedded_boundary_diffraction/inputs_rz b/Examples/Tests/embedded_boundary_diffraction/inputs_rz
@@ -0,0 +1,42 @@
+# This script tests the diffraction of a laser by
+# a cylindrical object, represented by an embedded boundary here
+
+max_step = 300
+amr.n_cell = 128 256
+amr.max_grid_size = 256
+amr.max_level = 0
+
+geometry.dims = RZ
+geometry.prob_lo     =  0. -0.2
+geometry.prob_hi     =  2   1.4
+warpx.cfl = 0.99
+algo.particle_shape = 1
+
+boundary.field_lo = none absorbing_silver_mueller
+boundary.field_hi = pec absorbing_silver_mueller
+
+# Make the cylindrical object that the laser will diffract on
+my_constants.aperture_l = 0.01
+my_constants.aperture_r = 0.2
+warpx.eb_implicit_function = "if( (abs(z)<0.5*aperture_l) and (x<aperture_r), 1, -1 )"
+
+warpx.n_rz_azimuthal_modes = 2
+
+# Laser
+lasers.names        = laser1
+laser1.profile      = Gaussian
+laser1.position     = 0. 0. -0.1
+laser1.direction    = 0. 0. 1.
+laser1.polarization = 1. 0. 0.
+laser1.profile_waist = 1.
+laser1.profile_duration = 100
+laser1.profile_t_peak = 0
+laser1.profile_focal_distance = 0
+laser1.e_max        = 1.
+laser1.wavelength   = 0.1
+
+diagnostics.diags_names = diag1
+diag1.intervals = 100
+diag1.diag_type = Full
+diag1.fields_to_plot = Er Et Ez Br Bt Bz
+diag1.format = openpmd
diff --git a/Regression/Checksum/benchmarks_json/EmbeddedBoundaryDiffraction.json b/Regression/Checksum/benchmarks_json/EmbeddedBoundaryDiffraction.json
@@ -0,0 +1,10 @@
+{
+  "lev=0": {
+    "Br": 6.821267675779345e-19,
+    "Bt": 5.564905732478707e-05,
+    "Bz": 2.368259586613272e-19,
+    "Er": 16503.98082446463,
+    "Et": 1.5299584682447838e-10,
+    "Ez": 1466.854467399728
+  }
+}
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
@@ -499,6 +499,24 @@ doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/analysis_default_regression.py
 
+[EmbeddedBoundaryDiffraction]
+buildDir = .
+inputFile = Examples/Tests/embedded_boundary_diffraction/inputs_rz
+runtime_params =
+dim = 2
+addToCompileString = USE_EB=TRUE USE_RZ=TRUE
+cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_EB=ON
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+outputFile = EmbeddedBoundaryDiffraction_plt
+analysisRoutine = Examples/Tests/embedded_boundary_diffraction/analysis_fields.py
+
 [ElectrostaticSphereEB_RZ]
 buildDir = .
 inputFile = Examples/Tests/electrostatic_sphere_eb/inputs_rz

diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp
@@ -180,11 +180,10 @@ WarpX::ComputeEdgeLengths (std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ed
 void
 WarpX::ComputeFaceAreas (std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face_areas,
                          const amrex::EBFArrayBoxFactory& eb_fact) {
-#ifndef WARPX_DIM_RZ
     BL_PROFILE("ComputeFaceAreas");
 
     auto const &flags = eb_fact.getMultiEBCellFlagFab();
-#ifdef WARPX_DIM_XZ
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
     //In 2D the volume frac is actually the area frac.
     auto const &area_frac = eb_fact.getVolFrac();
 #elif defined(WARPX_DIM_3D)
@@ -194,12 +193,12 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face
         "ComputeFaceAreas: Only implemented in 2D3V and 3D3V");
 #endif
 
-#ifdef WARPX_DIM_XZ
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
     face_areas[0]->setVal(0.);
     face_areas[2]->setVal(0.);
 #endif
     for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
-#ifdef WARPX_DIM_XZ
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
         // In 2D we change the extrema of the for loop so that we only have the case idim=1
         for (int idim = 1; idim < AMREX_SPACEDIM; ++idim) {
 #elif defined(WARPX_DIM_3D)
@@ -223,7 +222,7 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face
                     face_areas_dim(i, j, k) = amrex::Real(0.);
                 });
             } else {
-#ifdef WARPX_DIM_XZ
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
                 auto const &face = area_frac.const_array(mfi);
 #elif defined(WARPX_DIM_3D)
                 auto const &face = area_frac[idim]->const_array(mfi);
@@ -237,7 +236,6 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face
             }
         }
     }
-#endif
 }
 
 
@@ -269,13 +267,12 @@ WarpX::ScaleEdges (std::array< std::unique_ptr<amrex::MultiFab>, 3 >& edge_lengt
 void
 WarpX::ScaleAreas(std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face_areas,
                   const std::array<amrex::Real,3>& cell_size) {
-#ifndef WARPX_DIM_RZ
     BL_PROFILE("ScaleAreas");
 
     amrex::Real full_area;
 
     for (amrex::MFIter mfi(*face_areas[0]); mfi.isValid(); ++mfi) {
-#ifdef WARPX_DIM_XZ
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
         // In 2D we change the extrema of the for loop so that we only have the case idim=1
         for (int idim = 1; idim < AMREX_SPACEDIM; ++idim) {
 #elif defined(WARPX_DIM_3D)
@@ -286,7 +283,7 @@ WarpX::ScaleAreas(std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face_areas,
 #endif
             const amrex::Box& box = mfi.tilebox(face_areas[idim]->ixType().toIntVect(),
                                                 face_areas[idim]->nGrowVect() );
-#ifdef WARPX_DIM_XZ
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
             full_area = cell_size[0]*cell_size[2];
 #elif defined(WARPX_DIM_3D)
             if (idim == 0) {
@@ -308,7 +305,6 @@ WarpX::ScaleAreas(std::array< std::unique_ptr<amrex::MultiFab>, 3 >& face_areas,
 
         }
     }
-#endif
 }
 
 

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -59,8 +59,8 @@ void FiniteDifferenceSolver::EvolveB (
     std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing,
     int lev, amrex::Real const dt ) {
 
-#ifndef AMREX_USE_EB
-    amrex::ignore_unused(area_mod, ECTRhofield, Venl, flag_info_cell, borrowing);
+#if defined(WARPX_DIM_RZ) || !defined(AMREX_USE_EB)
+  amrex::ignore_unused(area_mod, ECTRhofield, Venl, flag_info_cell, borrowing);
 #endif
 
     // Select algorithm (The choice of algorithm is a runtime option,

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
@@ -67,8 +67,7 @@ void FiniteDifferenceSolver::EvolveE (
     // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
     if (m_fdtd_algo == ElectromagneticSolverAlgo::Yee){
-        ignore_unused(edge_lengths);
-        EvolveECylindrical <CylindricalYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
+        EvolveECylindrical <CylindricalYeeAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
 #else
     if (m_grid_type == GridType::Collocated) {
 
@@ -181,7 +180,6 @@ void FiniteDifferenceSolver::EvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-
 #ifdef AMREX_USE_EB
                 // Skip field push if this cell is fully covered by embedded boundaries
                 if (lz(i,j,k) <= 0) return;
@@ -235,9 +233,14 @@ void FiniteDifferenceSolver::EvolveECylindrical (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
     int lev, amrex::Real const dt ) {
 
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(edge_lengths);
+#endif
+
     amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
     // Loop through the grids, and over the tiles within each grid
@@ -262,6 +265,11 @@ void FiniteDifferenceSolver::EvolveECylindrical (
         Array4<Real> const& jt = Jfield[1]->array(mfi);
         Array4<Real> const& jz = Jfield[2]->array(mfi);
 
+#ifdef AMREX_USE_EB
+        amrex::Array4<amrex::Real> const& lr = edge_lengths[0]->array(mfi);
+        amrex::Array4<amrex::Real> const& lz = edge_lengths[2]->array(mfi);
+#endif
+
         // Extract stencil coefficients
         Real const * const AMREX_RESTRICT coefs_r = m_stencil_coefs_r.dataPtr();
         auto const n_coefs_r = static_cast<int>(m_stencil_coefs_r.size());
@@ -284,6 +292,10 @@ void FiniteDifferenceSolver::EvolveECylindrical (
         amrex::ParallelFor(ter, tet, tez,
 
             [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){
+#ifdef AMREX_USE_EB
+                // Skip field push if this cell is fully covered by embedded boundaries
+                if (lr(i, j, 0) <= 0) return;
+#endif
                 Real const r = rmin + (i + 0.5_rt)*dr; // r on cell-centered point (Er is cell-centered in r)
                 Er(i, j, 0, 0) +=  c2 * dt*(
                     - T_Algo::DownwardDz(Bt, coefs_z, n_coefs_z, i, j, 0, 0)
@@ -301,6 +313,11 @@ void FiniteDifferenceSolver::EvolveECylindrical (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){
+#ifdef AMREX_USE_EB
+                // Skip field push if this cell is fully covered by embedded boundaries
+                // The Et field is at a node, so we need to check if the node is covered
+                if (lr(i, j, 0)<=0 || lr(i-1, j, 0)<=0 || lz(i, j-1, 0)<=0 || lz(i, j, 0)<=0) return;
+#endif
                 Real const r = rmin + i*dr; // r on a nodal grid (Et is nodal in r)
                 if (r != 0) { // Off-axis, regular Maxwell equations
                     Et(i, j, 0, 0) += c2 * dt*(
@@ -342,6 +359,10 @@ void FiniteDifferenceSolver::EvolveECylindrical (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){
+#ifdef AMREX_USE_EB
+                // Skip field push if this cell is fully covered by embedded boundaries
+                if (lz(i, j, 0) <= 0) return;
+#endif
                 Real const r = rmin + i*dr; // r on a nodal grid (Ez is nodal in r)
                 if (r != 0) { // Off-axis, regular Maxwell equations
                     Ez(i, j, 0, 0) += c2 * dt*(

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -213,6 +213,7 @@ class FiniteDifferenceSolver
             std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
             std::unique_ptr<amrex::MultiFab> const& Ffield,
             int lev,
             amrex::Real dt );

diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
@@ -750,13 +750,6 @@ WarpX::ReadParameters ()
             electromagnetic_solver_id = ElectromagneticSolverAlgo::None;
         }
 
-#if defined(AMREX_USE_EB) && defined(WARPX_DIM_RZ)
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-        electromagnetic_solver_id==ElectromagneticSolverAlgo::None
-        || electromagnetic_solver_id==ElectromagneticSolverAlgo::HybridPIC,
-        "Currently, the embedded boundary in RZ only works for electrostatic solvers, the Ohm's law solver or with no solver installed.");
-#endif
-
         if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame ||
             electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic)
         {