From e38acc4da155087b3cdccb502fc808639c0c5c30 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 8 Mar 2024 18:12:04 -0800 Subject: [PATCH] Implement stair-case Yee solver with EB in RZ geometry (#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 5087485db606f77806c08849c51af9be635ca622. * 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 be3039af229d5dee509c8de0a712addf099a2b82. * Change lx to lr --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: lgiacome --- .../analysis_fields.py | 41 ++++++++++++++++++ .../embedded_boundary_diffraction/inputs_rz | 42 +++++++++++++++++++ .../EmbeddedBoundaryDiffraction.json | 10 +++++ Regression/WarpX-tests.ini | 18 ++++++++ Source/EmbeddedBoundary/WarpXInitEB.cpp | 16 +++---- .../FiniteDifferenceSolver/EvolveB.cpp | 4 +- .../FiniteDifferenceSolver/EvolveE.cpp | 27 ++++++++++-- .../FiniteDifferenceSolver.H | 1 + Source/WarpX.cpp | 7 ---- 9 files changed, 144 insertions(+), 22 deletions(-) create mode 100755 Examples/Tests/embedded_boundary_diffraction/analysis_fields.py create mode 100644 Examples/Tests/embedded_boundary_diffraction/inputs_rz create mode 100644 Regression/Checksum/benchmarks_json/EmbeddedBoundaryDiffraction.json diff --git a/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py b/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py new file mode 100755 index 00000000000..da344f332a1 --- /dev/null +++ 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 new file mode 100644 index 00000000000..705c9641dfb --- /dev/null +++ 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, 3 >& ed void WarpX::ComputeFaceAreas (std::array< std::unique_ptr, 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, 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, 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, 3 >& face } } } -#endif } @@ -269,13 +267,12 @@ WarpX::ScaleEdges (std::array< std::unique_ptr, 3 >& edge_lengt void WarpX::ScaleAreas(std::array< std::unique_ptr, 3 >& face_areas, const std::array& 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, 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, 3 >& face_areas, } } -#endif } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index 4425e655917..fbc1397b413 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -59,8 +59,8 @@ void FiniteDifferenceSolver::EvolveB ( std::array< std::unique_ptr >, 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 index e52cca846bb..5f707fbc927 100644 --- 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 ( Efield, Bfield, Jfield, Ffield, lev, dt ); + EvolveECylindrical ( 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, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real const dt ) { +#ifndef AMREX_USE_EB + amrex::ignore_unused(edge_lengths); +#endif + amrex::LayoutData* cost = WarpX::getCosts(lev); // Loop through the grids, and over the tiles within each grid @@ -262,6 +265,11 @@ void FiniteDifferenceSolver::EvolveECylindrical ( Array4 const& jt = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); +#ifdef AMREX_USE_EB + amrex::Array4 const& lr = edge_lengths[0]->array(mfi); + amrex::Array4 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(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 index 861b2648c1e..d045a30dd44 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -213,6 +213,7 @@ class FiniteDifferenceSolver std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real dt ); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d453fa6f941..6c3186c5c8d 100644 --- 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) {