aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Evolve/WarpXEvolve.cpp2
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp155
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt1
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H6
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/Make.package1
-rw-r--r--Source/FieldSolver/WarpXPushFieldsEM.cpp7
-rw-r--r--Source/WarpX.H2
-rw-r--r--Source/WarpX.cpp4
8 files changed, 178 insertions, 0 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 8e67af998..d2392f21b 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -358,7 +358,9 @@ WarpX::OneStep_nosub (Real cur_time)
FillBoundaryF(guard_cells.ng_FieldSolverF);
EvolveB(0.5_rt * dt[0]); // We now have B^{n+1/2}
+ if (do_silver_mueller) ApplySilverMuellerBoundary( dt[0] );
FillBoundaryB(guard_cells.ng_FieldSolver);
+
if (WarpX::em_solver_medium == MediumForEM::Vacuum) {
// vacuum medium
EvolveE(dt[0]); // We now have E^{n+1}
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp
new file mode 100644
index 000000000..abceead1f
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp
@@ -0,0 +1,155 @@
+/* Copyright 2021
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "FiniteDifferenceSolver.H"
+#include "Utils/WarpXConst.H"
+#include <AMReX_Gpu.H>
+
+using namespace amrex;
+
+/**
+ * \brief Update the B field at the boundary, using the Silver-Mueller condition
+ */
+void FiniteDifferenceSolver::ApplySilverMuellerBoundary (
+ std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
+ std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+ amrex::Box domain_box,
+ amrex::Real const dt ) {
+
+#ifdef WARPX_DIM_RZ
+ amrex::Abort("The Silver-Mueller boundary conditions cannot be used in RZ geometry.");
+#else
+
+ // Ensure that we are using the Yee solver
+ if (m_fdtd_algo != MaxwellSolverAlgo::Yee) {
+ amrex::Abort("The Silver-Mueller boundary conditions can only be used with the Yee solver.");
+ }
+
+ // Ensure that we are using the cells the domain
+ domain_box.enclosedCells();
+
+ // Calculate relevant coefficients
+ amrex::Real const cdt_over_dx = PhysConst::c*dt*m_stencil_coefs_x[0];
+ amrex::Real const coef1_x = (1._rt - cdt_over_dx)/(1._rt + cdt_over_dx);
+ amrex::Real const coef2_x = 2._rt*cdt_over_dx/(1._rt + cdt_over_dx) / PhysConst::c;
+#ifdef WARPX_DIM_3D
+ amrex::Real const cdt_over_dy = PhysConst::c*dt*m_stencil_coefs_y[0];
+ amrex::Real const coef1_y = (1._rt - cdt_over_dy)/(1._rt + cdt_over_dx);
+ amrex::Real const coef2_y = 2._rt*cdt_over_dy/(1._rt + cdt_over_dy) / PhysConst::c;
+#endif
+ amrex::Real const cdt_over_dz = PhysConst::c*dt*m_stencil_coefs_z[0];
+ amrex::Real const coef1_z = (1._rt - cdt_over_dz)/(1._rt + cdt_over_dx);
+ amrex::Real const coef2_z = 2._rt*cdt_over_dz/(1._rt + cdt_over_dz) / PhysConst::c;
+
+ // Loop through the grids, and over the tiles within each grid
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+ // tiling is usually set by TilingIfNotGPU()
+ // but here, we set it to false because of potential race condition,
+ // since we grow the tiles by one guard cell after creating them.
+ for ( MFIter mfi(*Efield[0], false); mfi.isValid(); ++mfi ) {
+
+ // Extract field data for this grid/tile
+ Array4<Real> const& Ex = Efield[0]->array(mfi);
+ Array4<Real> const& Ey = Efield[1]->array(mfi);
+ Array4<Real> const& Ez = Efield[2]->array(mfi);
+ Array4<Real> const& Bx = Bfield[0]->array(mfi);
+ Array4<Real> const& By = Bfield[1]->array(mfi);
+ Array4<Real> const& Bz = Bfield[2]->array(mfi);
+
+ // Extract the tileboxes for which to loop
+ Box tbx = mfi.tilebox(Bfield[0]->ixType().toIntVect());
+ Box tby = mfi.tilebox(Bfield[1]->ixType().toIntVect());
+ Box tbz = mfi.tilebox(Bfield[2]->ixType().toIntVect());
+
+ // We will modify the first (i.e. innermost) guard cell
+ // (if it is outside of the physical domain)
+ // Thus, the tileboxes here are grown by 1 guard cell
+ tbx.grow(1);
+ tby.grow(1);
+ tbz.grow(1);
+
+ // Loop over cells
+ amrex::ParallelFor(tbx, tby, tbz,
+
+ // Apply Boundary condition to Bx
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+#ifdef WARPX_DIM_3D
+ // At the +y boundary (innermost guard cell)
+ if ( j==domain_box.bigEnd(1)+1 )
+ Bx(i,j,k) = coef1_y * Bx(i,j,k) + coef2_y * Ez(i,j,k);
+ // At the -y boundary (innermost guard cell)
+ if ( j==domain_box.smallEnd(1)-1 )
+ Bx(i,j,k) = coef1_y * Bx(i,j,k) - coef2_y * Ez(i,j+1,k);
+ // At the +z boundary (innermost guard cell)
+ if ( k==domain_box.bigEnd(2)+1 )
+ Bx(i,j,k) = coef1_z * Bx(i,j,k) - coef2_z * Ey(i,j,k);
+ // At the -z boundary (innermost guard cell)
+ if ( k==domain_box.smallEnd(2)-1 )
+ Bx(i,j,k) = coef1_z * Bx(i,j,k) + coef2_z * Ey(i,j,k+1);
+#elif WARPX_DIM_XZ
+ // At the +z boundary (innermost guard cell)
+ if ( j==domain_box.bigEnd(1)+1 )
+ Bx(i,j,k) = coef1_z * Bx(i,j,k) - coef2_z * Ey(i,j,k);
+ // At the -z boundary (innermost guard cell)
+ if ( j==domain_box.smallEnd(1)-1 )
+ Bx(i,j,k) = coef1_z * Bx(i,j,k) + coef2_z * Ey(i,j+1,k);
+#endif
+ },
+
+ // Apply Boundary condition to By
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+ // At the +x boundary (innermost guard cell)
+ if ( i==domain_box.bigEnd(0)+1 )
+ By(i,j,k) = coef1_x * By(i,j,k) - coef2_x * Ez(i,j,k);
+ // At the -x boundary (innermost guard cell)
+ if ( i==domain_box.smallEnd(0)-1 )
+ By(i,j,k) = coef1_x * By(i,j,k) + coef2_x * Ez(i+1,j,k);
+#ifdef WARPX_DIM_3D
+ // At the +z boundary (innermost guard cell)
+ if ( k==domain_box.bigEnd(2)+1 )
+ By(i,j,k) = coef1_z * By(i,j,k) + coef2_z * Ex(i,j,k);
+ // At the -z boundary (innermost guard cell)
+ if ( k==domain_box.smallEnd(2)-1 )
+ By(i,j,k) = coef1_z * By(i,j,k) - coef2_z * Ex(i,j,k+1);
+#elif WARPX_DIM_XZ
+ // At the +z boundary (innermost guard cell)
+ if ( j==domain_box.bigEnd(1)+1 )
+ By(i,j,k) = coef1_z * By(i,j,k) + coef2_z * Ex(i,j,k);
+ // At the -z boundary (innermost guard cell)
+ if ( j==domain_box.smallEnd(1)-1 )
+ By(i,j,k) = coef1_z * By(i,j,k) - coef2_z * Ex(i,j+1,k);
+#endif
+ },
+
+ // Apply Boundary condition to Bz
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+ // At the +x boundary (innermost guard cell)
+ if ( i==domain_box.bigEnd(0)+1 )
+ Bz(i,j,k) = coef1_x * Bz(i,j,k) + coef2_x * Ey(i,j,k);
+ // At the -x boundary (innermost guard cell)
+ if ( i==domain_box.smallEnd(0)-1 )
+ Bz(i,j,k) = coef1_x * Bz(i,j,k) - coef2_x * Ey(i+1,j,k);
+#ifdef WARPX_DIM_3D
+ // At the +y boundary (innermost guard cell)
+ if ( j==domain_box.bigEnd(1)+1 )
+ Bz(i,j,k) = coef1_y * Bz(i,j,k) - coef2_y * Ex(i,j,k);
+ // At the -y boundary (innermost guard cell)
+ if ( j==domain_box.smallEnd(1)-1 )
+ Bz(i,j,k) = coef1_y * Bz(i,j,k) + coef2_y * Ex(i,j+1,k);
+#endif
+ }
+ );
+
+ }
+#endif // WARPX_DIM_RZ
+}
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt
index 93dbd7938..4b1388f42 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt
+++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt
@@ -9,6 +9,7 @@ target_sources(WarpX
EvolveFPML.cpp
FiniteDifferenceSolver.cpp
MacroscopicEvolveE.cpp
+ ApplySilverMuellerBoundary.cpp
)
add_subdirectory(MacroscopicProperties)
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index 72fabb54e..c5041a53f 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -52,6 +52,12 @@ class FiniteDifferenceSolver
int const rhocomp,
amrex::Real const dt );
+ void ApplySilverMuellerBoundary(
+ std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
+ std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+ amrex::Box domain_box,
+ amrex::Real const dt );
+
void ComputeDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
amrex::MultiFab& divE );
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
index 31bccd1e4..ea94f2c66 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package
+++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
@@ -7,6 +7,7 @@ CEXE_sources += MacroscopicEvolveE.cpp
CEXE_sources += EvolveBPML.cpp
CEXE_sources += EvolveEPML.cpp
CEXE_sources += EvolveFPML.cpp
+CEXE_sources += ApplySilverMuellerBoundary.cpp
include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/Make.package
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 1e138b4c1..120dc8b5c 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -211,6 +211,13 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
}
void
+WarpX::ApplySilverMuellerBoundary (amrex::Real a_dt) {
+ // Only apply to level 0
+ m_fdtd_solver_fp[0]->ApplySilverMuellerBoundary(
+ Efield_fp[0], Bfield_fp[0], Geom(0).Domain(), a_dt );
+}
+
+void
WarpX::EvolveE (amrex::Real a_dt)
{
for (int lev = 0; lev <= finest_level; ++lev)
diff --git a/Source/WarpX.H b/Source/WarpX.H
index a03f0e717..9cfbc3b72 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -327,6 +327,7 @@ public:
void EvolveB (int lev, PatchType patch_type, amrex::Real dt);
void EvolveE (int lev, PatchType patch_type, amrex::Real dt);
void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type);
+ void ApplySilverMuellerBoundary (amrex::Real dt);
void MacroscopicEvolveE ( amrex::Real dt);
void MacroscopicEvolveE (int lev, amrex::Real dt);
@@ -820,6 +821,7 @@ private:
// PML
int do_pml = 1;
+ int do_silver_mueller = 0;
int pml_ncell = 10;
int pml_delta = 10;
int pml_has_particles = 0;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 2fc6fabc6..b4528b1bc 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -554,6 +554,10 @@ WarpX::ReadParameters ()
}
pp.query("do_pml", do_pml);
+ pp.query("do_silver_mueller", do_silver_mueller);
+ if ( (do_pml==1)&&(do_silver_mueller==1) ) {
+ amrex::Abort("PML and Silver-Mueller boundary conditions cannot be activated at the same time.");
+ }
pp.query("pml_ncell", pml_ncell);
pp.query("pml_delta", pml_delta);
pp.query("pml_has_particles", pml_has_particles);
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-14 02:53:51 +0000