From a4d563147edf63b11834f2a8900b28b32dd80d10 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 5 Nov 2024 13:25:11 -0800 Subject: [PATCH] Add `PECInsulator` boundary condition (#4943) This PR adds a mixed PEC and insulator boundary condition. This allows an insulator to be placed on a portion of the boundary. The rest of that boundary will be PEC. Within the insulator portion, the tangential fields can be specified on the boundary (as functions of space and time). The normal fields and fields not specified are extrapolated to the guard cells from the valid cells. The fields are specified in pairs, the two tangential electric fields, and the two tangential magnetic fields. In each pair, if one is set, the other will be zeroed if not set. A use case is the simulation of a dynamic pinch, driven by an external current, represented as a time dependent B field on the boundary. [PECinsulatorBC_warpX_summaryOnly.pdf](https://github.com/user-attachments/files/17637695/PECinsulatorBC_warpX_summaryOnly.pdf) --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> --- Docs/source/usage/parameters.rst | 31 + Examples/Tests/pec/CMakeLists.txt | 10 + .../pec/inputs_test_2d_pec_field_insulator | 34 ++ .../test_2d_pec_field_insulator.json | 13 + Source/BoundaryConditions/CMakeLists.txt | 1 + Source/BoundaryConditions/Make.package | 1 + Source/BoundaryConditions/PEC_Insulator.H | 180 ++++++ Source/BoundaryConditions/PEC_Insulator.cpp | 561 ++++++++++++++++++ Source/BoundaryConditions/PEC_Insulator_fwd.H | 13 + .../WarpXFieldBoundaries.cpp | 63 ++ Source/Python/WarpX.cpp | 1 + Source/Utils/Parser/ParserUtils.H | 14 + Source/Utils/Parser/ParserUtils.cpp | 13 + Source/Utils/WarpXAlgorithmSelection.H | 1 + Source/WarpX.H | 4 + Source/WarpX.cpp | 4 + 16 files changed, 944 insertions(+) create mode 100644 Examples/Tests/pec/inputs_test_2d_pec_field_insulator create mode 100644 Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json create mode 100644 Source/BoundaryConditions/PEC_Insulator.H create mode 100644 Source/BoundaryConditions/PEC_Insulator.cpp create mode 100644 Source/BoundaryConditions/PEC_Insulator_fwd.H diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index af559aa1fba..1b5c3e7b186 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -493,6 +493,37 @@ Domain Boundary Conditions * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. Please see the :ref:`PEC theory section ` for more details. Note that PEC boundary is invalid at `r=0` for the RZ solver. Please use ``none`` option. This boundary condition does not work with the spectral solver. + * ``pec_insulator``: This option specifies a mixed perfect electric conductor and insulator boundary, where some part of the + boundary is PEC and some is insulator. In the insulator portion, the normal fields are extrapolated and the tangential fields + are either set to the specified value or extrapolated. The region that is insulator is specified using a spatially dependent expression with the insulator being in the area where the value of the expression is greater than zero. + The expressions are given for the low and high boundary on each axis, as listed below. The tangential fields are specified as + expressions that can depend on the location and time. The tangential fields are in two pairs, the electric fields and the + magnetic fields. In each pair, if one is specified, the other will be set to zero if not also specified. + + * ``insulator.area_x_lo(y,z)``: For the lower x (or r) boundary, expression specifying the insulator location + + * ``insulator.area_x_hi(y,z)``: For the upper x (or r) boundary, expression specifying the insulator location + + * ``insulator.area_y_lo(x,z)``: For the lower y boundary, expression specifying the insulator location + + * ``insulator.area_y_hi(x,z)``: For the upper y boundary, expression specifying the insulator location + + * ``insulator.area_z_lo(x,y)``: For the lower z boundary, expression specifying the insulator location + + * ``insulator.area_z_hi(x,y)``: For the upper z boundary, expression specifying the insulator location + + * ``insulator.Ey_x_lo(y,z,t)``, ``insulator.Ez_x_lo(y,z,t)``, ``insulator.By_x_lo(y,z,t)``, ``insulator.Bz_x_lo(y,z,t)``: expressions of the tangential field values for the lower x (or r) boundary + + * ``insulator.Ey_x_hi(y,z,t)``, ``insulator.Ez_x_hi(y,z,t)``, ``insulator.By_x_hi(y,z,t)``, ``insulator.Bz_x_hi(y,z,t)``: expressions of the tangential field values for the upper x (or r) boundary + + * ``insulator.Ex_y_lo(x,z,t)``, ``insulator.Ez_y_lo(x,z,t)``, ``insulator.Bx_y_lo(x,z,t)``, ``insulator.Bz_y_lo(x,z,t)``: expressions of the tangential field values for the lower y boundary + + * ``insulator.Ex_y_hi(x,z,t)``, ``insulator.Ez_y_hi(x,z,t)``, ``insulator.Bx_y_hi(x,z,t)``, ``insulator.Bz_y_hi(x,z,t)``: expressions of the tangential field values for the upper y boundary + + * ``insulator.Ex_z_lo(x,y,t)``, ``insulator.Ey_z_lo(x,y,t)``, ``insulator.Bx_z_lo(x,y,t)``, ``insulator.By_z_lo(x,y,t)``: expressions of the tangential field values for the lower z boundary + + * ``insulator.Ex_z_hi(x,y,t)``, ``insulator.Ey_z_hi(x,y,t)``, ``insulator.Bx_z_hi(x,y,t)``, ``insulator.By_z_hi(x,y,t)``: expressions of the tangential field values for the upper z boundary + * ``none``: No boundary condition is applied to the fields with the electromagnetic solver. This option must be used for the RZ-solver at `r=0`. * ``neumann``: For the electrostatic multigrid solver, a Neumann boundary condition (with gradient of the potential equal to 0) will be applied on the specified boundary. diff --git a/Examples/Tests/pec/CMakeLists.txt b/Examples/Tests/pec/CMakeLists.txt index ec710f7d919..e0bab40d058 100644 --- a/Examples/Tests/pec/CMakeLists.txt +++ b/Examples/Tests/pec/CMakeLists.txt @@ -30,3 +30,13 @@ add_warpx_test( diags/diag1000020 # output OFF # dependency ) + +add_warpx_test( + test_2d_pec_field_insulator # name + 2 # dims + 2 # nprocs + inputs_test_2d_pec_field_insulator # inputs + analysis_default_regression.py # analysis + diags/diag1000010 # output + OFF # dependency +) diff --git a/Examples/Tests/pec/inputs_test_2d_pec_field_insulator b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator new file mode 100644 index 00000000000..68a8df1b600 --- /dev/null +++ b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator @@ -0,0 +1,34 @@ +# Maximum number of time steps +max_step = 10 + +# number of grid points +amr.n_cell = 32 32 +amr.blocking_factor = 16 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +# Geometry +geometry.dims = 2 +geometry.prob_lo = 0. 2.e-2 # physical domain +geometry.prob_hi = 1.e-2 3.e-2 + +# Boundary condition +boundary.field_lo = neumann periodic +boundary.field_hi = PECInsulator periodic + +warpx.serialize_initial_conditions = 1 + +# Verbosity +warpx.verbose = 1 + +# CFL +warpx.cfl = 1.0 + +insulator.area_x_hi(y,z) = (2.25e-2 <= z and z <= 2.75e-2) +insulator.By_x_hi(y,z,t) = min(t/1.0e-12,1)*1.e1*3.3e-4 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 10 +diag1.diag_type = Full diff --git a/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json b/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json new file mode 100644 index 00000000000..ca6f38977ae --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator.json @@ -0,0 +1,13 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.34938851065132936, + "Bz": 0.0, + "Ex": 31871402.236828588, + "Ey": 0.0, + "Ez": 104908439.18998256, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + } +} \ No newline at end of file diff --git a/Source/BoundaryConditions/CMakeLists.txt b/Source/BoundaryConditions/CMakeLists.txt index 751e52abdd9..c560d121385 100644 --- a/Source/BoundaryConditions/CMakeLists.txt +++ b/Source/BoundaryConditions/CMakeLists.txt @@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS) warpx_set_suffix_dims(SD ${D}) target_sources(lib_${SD} PRIVATE + PEC_Insulator.cpp PML.cpp WarpXEvolvePML.cpp WarpXFieldBoundaries.cpp diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package index 43d18425ffc..452c9c18b7e 100644 --- a/Source/BoundaryConditions/Make.package +++ b/Source/BoundaryConditions/Make.package @@ -1,3 +1,4 @@ +CEXE_sources += PEC_Insulator.cpp CEXE_sources += PML.cpp WarpXEvolvePML.cpp CEXE_sources += WarpXFieldBoundaries.cpp WarpX_PEC.cpp diff --git a/Source/BoundaryConditions/PEC_Insulator.H b/Source/BoundaryConditions/PEC_Insulator.H new file mode 100644 index 00000000000..5cfdf6488f0 --- /dev/null +++ b/Source/BoundaryConditions/PEC_Insulator.H @@ -0,0 +1,180 @@ +#ifndef PEC_INSULATOR_H_ +#define PEC_INSULATOR_H_ + +#include "Utils/WarpXAlgorithmSelection.H" + +#include +#include +#include + +#include + +#include +#include + +class PEC_Insulator +{ +public: + + PEC_Insulator(); + + /** + * \brief Apply either the PEC or insulator boundary condition on the boundary and in the + * guard cells. + * In the PEC, the nodal fields (in a Yee mesh) are made even relative to the boundary, + * the non-nodal fields are made odd. + * In the insulator, the tangential fields are set to the value if specified, otherwise unchanged, + * and the normal fields extrapolated from the valid cells. + * + * \param[in,out] Efield + * \param[in] field_boundary_lo lower field boundary conditions + * \param[in] field_boundary_hi upper field boundary conditions + * \param[in] ng_fieldgather number of guard cells used by field gather + * \param[in] geom geometry object of level "lev" + * \param[in] lev level of the Multifab + * \param[in] patch_type coarse or fine + * \param[in] ref_ratios vector containing the refinement ratios of the refinement levels + * \param[in] time current time of the simulation + * \param[in] split_pml_field whether pml the multifab is the regular Efield or + * split pml field + */ + void ApplyPEC_InsulatortoEfield (std::array Efield, + amrex::Array const & field_boundary_lo, + amrex::Array const & field_boundary_hi, + amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom, + int lev, PatchType patch_type, amrex::Vector const & ref_ratios, + amrex::Real time, + bool split_pml_field = false); + /** + * \brief Apply either the PEC or insulator boundary condition on the boundary and in the + * guard cells. + * In the PEC, the nodal fields (in a Yee mesh) are made even relative to the boundary, + * the non-nodal fields are made odd. + * In the insulator, the tangential fields are set to the value if specified, otherwise unchanged, + * and the normal fields extrapolated from the valid cells. + * + * \param[in,out] Bfield + * \param[in] field_boundary_lo lower field boundary conditions + * \param[in] field_boundary_hi upper field boundary conditions + * \param[in] ng_fieldgather number of guard cells used by field gather + * \param[in] geom geometry object of level "lev" + * \param[in] lev level of the Multifab + * \param[in] patch_type coarse or fine + * \param[in] ref_ratios vector containing the refinement ratios of the refinement levels + * \param[in] time current time of the simulation + */ + void ApplyPEC_InsulatortoBfield (std::array Bfield, + amrex::Array const & field_boundary_lo, + amrex::Array const & field_boundary_hi, + amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom, + int lev, PatchType patch_type, amrex::Vector const & ref_ratios, + amrex::Real time); + + /** + * \brief The work routine applying the boundary condition + * + * \param[in,out] field + * \param[in] field_boundary_lo lower field boundary conditions + * \param[in] field_boundary_hi upper field boundary conditions + * \param[in] ng_fieldgather number of guard cells used by field gather + * \param[in] geom geometry object of level "lev" + * \param[in] lev level of the Multifab + * \param[in] patch_type coarse or fine + * \param[in] ref_ratios vector containing the refinement ratios of the refinement levels + * \param[in] time current time of the simulation + * \param[in] split_pml_field whether pml the multifab is the regular Efield or + * split pml field + * \param[in] E_like whether the field is E like or B like + * \param[in] set_F_x_lo whether the tangential field at the boundary was specified + * \param[in] set_F_x_hi whether the tangential field at the boundary was specified + * \param[in] a_Fy_x_lo the parser for the tangential field at the boundary + * \param[in] a_Fz_x_lo the parser for the tangential field at the boundary + * \param[in] a_Fy_x_hi the parser for the tangential field at the boundary + * \param[in] a_Fz_x_hi the parser for the tangential field at the boundary + * \param[in] set_F_y_lo whether the tangential field at the boundary was specified + * \param[in] set_F_y_hi whether the tangential field at the boundary was specified + * \param[in] a_Fx_y_lo the parser for the tangential field at the boundary + * \param[in] a_Fz_y_lo the parser for the tangential field at the boundary + * \param[in] a_Fx_y_hi the parser for the tangential field at the boundary + * \param[in] a_Fz_y_hi the parser for the tangential field at the boundary + * \param[in] set_F_z_lo whether the tangential field at the boundary was specified + * \param[in] set_F_z_hi whether the tangential field at the boundary was specified + * \param[in] a_Fx_z_lo the parser for the tangential field at the boundary + * \param[in] a_Fy_z_lo the parser for the tangential field at the boundary + * \param[in] a_Fx_z_hi the parser for the tangential field at the boundary + * \param[in] a_Fy_z_hi the parser for the tangential field at the boundary + */ + void + ApplyPEC_InsulatortoField (std::array field, + amrex::Array const & field_boundary_lo, + amrex::Array const & field_boundary_hi, + amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom, + int lev, PatchType patch_type, amrex::Vector const & ref_ratios, + amrex::Real time, + bool split_pml_field, + bool E_like, +#if (AMREX_SPACEDIM > 1) + bool set_F_x_lo, bool set_F_x_hi, + std::unique_ptr const & a_Fy_x_lo, std::unique_ptr const & a_Fz_x_lo, + std::unique_ptr const & a_Fy_x_hi, std::unique_ptr const & a_Fz_x_hi, +#endif +#if defined(WARPX_DIM_3D) + bool set_F_y_lo, bool set_F_y_hi, + std::unique_ptr const & a_Fx_y_lo, std::unique_ptr const & a_Fz_y_lo, + std::unique_ptr const & a_Fx_y_hi, std::unique_ptr const & a_Fz_y_hi, +#endif + bool set_F_z_lo, bool set_F_z_hi, + std::unique_ptr const & a_Fx_z_lo, std::unique_ptr const & a_Fy_z_lo, + std::unique_ptr const & a_Fx_z_hi, std::unique_ptr const & a_Fy_z_hi); + +private: + + /* \brief Reads in the parsers for the tangential fields, returning whether + * the input parameter was specified. + * \param[in] pp_insulator ParmParse instance + * \param[out] parser the parser generated from the input + * \param[in] input_name the name of the input parameter + * \param[in] coord1 the first coordinate in the plane + * \param[in] coord2 the second coordinate in the plane + */ + bool ReadTangentialFieldParser (amrex::ParmParse const & pp_insulator, + std::unique_ptr & parser, + std::string const & input_name, + std::string const & coord1, + std::string const & coord2); + + std::vector> m_insulator_area_lo; + std::vector> m_insulator_area_hi; + +#if (AMREX_SPACEDIM > 1) + bool m_set_B_x_lo = false, m_set_B_x_hi = false; + std::unique_ptr m_By_x_lo, m_Bz_x_lo; + std::unique_ptr m_By_x_hi, m_Bz_x_hi; +#endif +#if defined(WARPX_DIM_3D) + bool m_set_B_y_lo = false, m_set_B_y_hi = false; + std::unique_ptr m_Bx_y_lo, m_Bz_y_lo; + std::unique_ptr m_Bx_y_hi, m_Bz_y_hi; +#endif + bool m_set_B_z_lo = false, m_set_B_z_hi = false; + std::unique_ptr m_Bx_z_lo, m_By_z_lo; + std::unique_ptr m_Bx_z_hi, m_By_z_hi; + + +#if (AMREX_SPACEDIM > 1) + bool m_set_E_x_lo = false, m_set_E_x_hi = false; + std::unique_ptr m_Ey_x_lo, m_Ez_x_lo; + std::unique_ptr m_Ey_x_hi, m_Ez_x_hi; +#endif +#if defined(WARPX_DIM_3D) + bool m_set_E_y_lo = false, m_set_E_y_hi = false; + std::unique_ptr m_Ex_y_lo, m_Ez_y_lo; + std::unique_ptr m_Ex_y_hi, m_Ez_y_hi; +#endif + bool m_set_E_z_lo = false, m_set_E_z_hi = false; + std::unique_ptr m_Ex_z_lo, m_Ey_z_lo; + std::unique_ptr m_Ex_z_hi, m_Ey_z_hi; + + +}; +#endif // PEC_INSULATOR_H_ diff --git a/Source/BoundaryConditions/PEC_Insulator.cpp b/Source/BoundaryConditions/PEC_Insulator.cpp new file mode 100644 index 00000000000..df411f8e908 --- /dev/null +++ b/Source/BoundaryConditions/PEC_Insulator.cpp @@ -0,0 +1,561 @@ +#include "BoundaryConditions/PEC_Insulator.H" +#include "Utils/Parser/ParserUtils.H" +#include "WarpX.H" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + /** + * \brief At the specified grid location, apply either the PEC or insulator boundary condition if + * the cell is on the boundary or in the guard cells. + * + * \param[in] icomp component of the field being updated + * (0=x, 1=y, 2=z in Cartesian) + * (0=r, 1=theta, 2=z in RZ) + * \param[in] dom_lo index value of the lower domain boundary (cell-centered) + * \param[in] dom_hi index value of the higher domain boundary (cell-centered) + * \param[in] ijk_vec indices along the x(i), y(j), z(k) of field Array4 + * \param[in] n index of the MultiFab component being updated + * \param[in] field field data to be updated if (ijk) is at the boundary + * or a guard cell + * \param[in] E_like whether the field behaves like E field or B field + * \param[in] is_nodal staggering of the field data being updated. + * \param[in] is_insulator_lo Specifies whether lower boundaries are insulators + * \param[in] is_insulator_hi Specifies whether upper boundaries are insulators + * \param[in] field_lo the values of the field for the lower insulator boundary cell + * \param[in] field_hi the values of the field for the upper insulator boundary cell + * \param[in] set_field_lo whether to set the field for the direction on the lower boundary + * \param[in] set_field_hi whether to set the field for the direction on the upper boundary + * \param[in] fbndry_lo specified values of the field at the lower boundaries in the insulator + * \param[in] fbndry_hi specified values of the field at the upper boundaries in the insulator + */ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void SetFieldOnPEC_Insulator (int icomp, + amrex::IntVect const & dom_lo, + amrex::IntVect const & dom_hi, + amrex::IntVect const & ijk_vec, int n, + amrex::Array4 const & field, + bool const E_like, + amrex::IntVect const & is_nodal, + amrex::IntVect const & is_insulator_lo, + amrex::IntVect const & is_insulator_hi, + amrex::RealVect const & field_lo, + amrex::RealVect const & field_hi, + amrex::IntVect const & set_field_lo, + amrex::IntVect const & set_field_hi, + amrex::GpuArray const fbndry_lo, + amrex::GpuArray const fbndry_hi) + { + using namespace amrex::literals; + amrex::IntVect ijk_mirror = ijk_vec; + amrex::IntVect ijk_mirrorp1 = ijk_vec; + bool OnBoundary = false; + bool GuardCell = false; + bool isInsulatorBoundary = false; + amrex::Real sign = +1._rt; + bool is_normal_to_boundary; + amrex::Real field_value = 0._rt; + bool set_field = false; + // Loop over all dimensions + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // Loop over sides, iside = -1 (lo), iside = +1 (hi) + for (int iside = -1; iside <= +1; iside += 2) { + bool const isPEC_InsulatorBoundary = ( (iside == -1) + ? fbndry_lo[idim] == FieldBoundaryType::PECInsulator + : fbndry_hi[idim] == FieldBoundaryType::PECInsulator ); + if (isPEC_InsulatorBoundary) { + isInsulatorBoundary = ( (iside == -1) + ? is_insulator_lo[idim] == 1 + : is_insulator_hi[idim] == 1 ); + } + if (isPEC_InsulatorBoundary) { + // Calculates the number of grid points ijk_vec is beyond the + // domain boundary i.e. a value of +1 means the current cell is + // outside of the simulation domain by 1 cell. Note that the high + // side domain boundary is between cell dom_hi and dom_hi+1 for cell + // centered grids and on cell dom_hi+1 for nodal grid. This is why + // (dom_hi[idim] + is_nodal[idim]) is used. + int const ig = ((iside == -1) ? (dom_lo[idim] - ijk_vec[idim]) + : (ijk_vec[idim] - (dom_hi[idim] + is_nodal[idim]))); + +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + // For 2D : for icomp==1, (Fy in XZ, Ftheta in RZ), + // icomp=1 is not normal to x or z boundary + // The logic below ensures that the flags are set right for 2D + is_normal_to_boundary = (icomp == (2*idim)); +#elif (defined WARPX_DIM_1D_Z) + // For 1D : icomp=0 and icomp=1 (Fx and Fy are not normal to the z boundary) + // The logic below ensures that the flags are set right for 1D + is_normal_to_boundary = (icomp == 2); +#else + is_normal_to_boundary = (icomp == idim); +#endif + + if (ig == 0) { + // Check if field is on the boundary + if (is_nodal[idim] == 1) { + OnBoundary = true; + } + } else if (ig > 0) { + GuardCell = true; + + // Mirror location inside the domain by "ig" number of cells + ijk_mirror[idim] = ( (iside == -1) + ? (dom_lo[idim] + ig - (1 - is_nodal[idim])) + : (dom_hi[idim] + 1 - ig)); + // Location twice as far in, for extrapolation + ijk_mirrorp1[idim] = 2*ijk_mirror[idim] - ijk_vec[idim]; + + // Check for components with even symmetry. + // True for E_like and tangential, and B_like and normal + if (E_like ^ is_normal_to_boundary) { sign *= -1._rt; } + + field_value = ( (iside == -1) ? field_lo[idim] : field_hi[idim] ); + set_field = ( (iside == -1) ? set_field_lo[idim]==1 : set_field_hi[idim]==1 ); + +#if (defined WARPX_DIM_RZ) + if (idim == 0 && iside == +1) { + // Upper radial boundary + amrex::Real const rguard = ijk_vec[idim] + 0.5_rt*(1._rt - is_nodal[idim]); + if (icomp == 0) { + // Add radial scale so that the divergence, drFr/dr, is 0. + // This only works for the first guard cell and with + // Fr cell centered in r. + amrex::Real const rmirror = ijk_mirror[idim] + 0.5_rt*(1._rt - is_nodal[idim]); + // Calculate radial scale factor + sign *= rmirror/rguard; + } + if (isInsulatorBoundary) { + // Apply radial scale factor + field_value *= dom_hi[idim]/rguard; + } + } +#endif + } + } // is pec_insulator boundary + } // loop over iside + } // loop over dimensions + + if (isInsulatorBoundary) { + if (is_normal_to_boundary) { + // The value on the boundary is left unmodified + // The values in the guard cells are extrapolated + if (GuardCell) { + field(ijk_vec, n) = 2._rt*field(ijk_mirror, n) - field(ijk_mirrorp1, n); + } + } else if ((OnBoundary || GuardCell) && set_field) { + field(ijk_vec, n) = field_value; + } else if (GuardCell) { + field(ijk_vec, n) = 2._rt*field(ijk_mirror, n) - field(ijk_mirrorp1, n); + } + } else { + if (OnBoundary && (E_like ^ is_normal_to_boundary)) { + // If ijk_vec is on a boundary, set to zero if + // E_like and tangential or B_like and normal + field(ijk_vec,n) = 0._rt; + } else if (GuardCell) { + // Fnormal and Ftangential is set opposite and equal to the value + // in the mirror location, respectively. + field(ijk_vec,n) = sign * field(ijk_mirror,n); + } + } + } +} + + +bool +PEC_Insulator::ReadTangentialFieldParser (amrex::ParmParse const & pp_insulator, + std::unique_ptr & parser, + std::string const & input_name, + std::string const & coord1, + std::string const & coord2) +{ + std::string str = "0"; + bool const specified = utils::parser::Query_parserString(pp_insulator, input_name, str); + parser = std::make_unique(utils::parser::makeParser(str, {coord1, coord2, "t"})); + return specified; +} + +PEC_Insulator::PEC_Insulator () +{ + + amrex::ParmParse const pp_insulator("insulator"); + +#if (AMREX_SPACEDIM > 1) + std::string str_area_x_lo = "0"; + std::string str_area_x_hi = "0"; + utils::parser::Query_parserString( pp_insulator, "area_x_lo(y,z)", str_area_x_lo); + utils::parser::Query_parserString( pp_insulator, "area_x_hi(y,z)", str_area_x_hi); + m_insulator_area_lo.push_back( + std::make_unique(utils::parser::makeParser(str_area_x_lo, {"y", "z"}))); + m_insulator_area_hi.push_back( + std::make_unique(utils::parser::makeParser(str_area_x_hi, {"y", "z"}))); + + m_set_B_x_lo |= ReadTangentialFieldParser(pp_insulator, m_By_x_lo, "By_x_lo(y,z,t)", "y", "z"); + m_set_B_x_lo |= ReadTangentialFieldParser(pp_insulator, m_Bz_x_lo, "Bz_x_lo(y,z,t)", "y", "z"); + m_set_B_x_hi |= ReadTangentialFieldParser(pp_insulator, m_By_x_hi, "By_x_hi(y,z,t)", "y", "z"); + m_set_B_x_hi |= ReadTangentialFieldParser(pp_insulator, m_Bz_x_hi, "Bz_x_hi(y,z,t)", "y", "z"); + + m_set_E_x_lo |= ReadTangentialFieldParser(pp_insulator, m_Ey_x_lo, "Ey_x_lo(y,z,t)", "y", "z"); + m_set_E_x_lo |= ReadTangentialFieldParser(pp_insulator, m_Ez_x_lo, "Ez_x_lo(y,z,t)", "y", "z"); + m_set_E_x_hi |= ReadTangentialFieldParser(pp_insulator, m_Ey_x_hi, "Ey_x_hi(y,z,t)", "y", "z"); + m_set_E_x_hi |= ReadTangentialFieldParser(pp_insulator, m_Ez_x_hi, "Ez_x_hi(y,z,t)", "y", "z"); +#endif +#if defined(WARPX_DIM_3D) + std::string str_area_y_lo = "0"; + std::string str_area_y_hi = "0"; + utils::parser::Query_parserString( pp_insulator, "area_y_lo(x,z)", str_area_y_lo); + utils::parser::Query_parserString( pp_insulator, "area_y_hi(x,z)", str_area_y_hi); + m_insulator_area_lo.push_back( + std::make_unique(utils::parser::makeParser(str_area_y_lo, {"x", "z"}))); + m_insulator_area_hi.push_back( + std::make_unique(utils::parser::makeParser(str_area_y_hi, {"x", "z"}))); + + m_set_B_y_lo |= ReadTangentialFieldParser(pp_insulator, m_Bx_y_lo, "Bx_y_lo(x,z,t)", "x", "z"); + m_set_B_y_lo |= ReadTangentialFieldParser(pp_insulator, m_Bz_y_lo, "Bz_y_lo(x,z,t)", "x", "z"); + m_set_B_y_hi |= ReadTangentialFieldParser(pp_insulator, m_Bx_y_hi, "Bx_y_hi(x,z,t)", "x", "z"); + m_set_B_y_hi |= ReadTangentialFieldParser(pp_insulator, m_Bz_y_hi, "Bz_y_hi(x,z,t)", "x", "z"); + + m_set_E_y_lo |= ReadTangentialFieldParser(pp_insulator, m_Ex_y_lo, "Ex_y_lo(x,z,t)", "x", "z"); + m_set_E_y_lo |= ReadTangentialFieldParser(pp_insulator, m_Ez_y_lo, "Ez_y_lo(x,z,t)", "x", "z"); + m_set_E_y_hi |= ReadTangentialFieldParser(pp_insulator, m_Ex_y_hi, "Ex_y_hi(x,z,t)", "x", "z"); + m_set_E_y_hi |= ReadTangentialFieldParser(pp_insulator, m_Ez_y_hi, "Ez_y_hi(x,z,t)", "x", "z"); +#endif + + std::string str_area_z_lo = "0"; + std::string str_area_z_hi = "0"; + utils::parser::Query_parserString( pp_insulator, "area_z_lo(x,y)", str_area_z_lo); + utils::parser::Query_parserString( pp_insulator, "area_z_hi(x,y)", str_area_z_hi); + m_insulator_area_lo.push_back( + std::make_unique(utils::parser::makeParser(str_area_z_lo, {"x", "y"}))); + m_insulator_area_hi.push_back( + std::make_unique(utils::parser::makeParser(str_area_z_hi, {"x", "y"}))); + + m_set_B_z_lo |= ReadTangentialFieldParser(pp_insulator, m_Bx_z_lo, "Bx_z_lo(x,y,t)", "x", "y"); + m_set_B_z_lo |= ReadTangentialFieldParser(pp_insulator, m_By_z_lo, "By_z_lo(x,y,t)", "x", "y"); + m_set_B_z_hi |= ReadTangentialFieldParser(pp_insulator, m_Bx_z_hi, "Bx_z_hi(x,y,t)", "x", "y"); + m_set_B_z_hi |= ReadTangentialFieldParser(pp_insulator, m_By_z_hi, "By_z_hi(x,y,t)", "x", "y"); + + m_set_E_z_lo |= ReadTangentialFieldParser(pp_insulator, m_Ex_z_lo, "Ex_z_lo(x,y,t)", "x", "y"); + m_set_E_z_lo |= ReadTangentialFieldParser(pp_insulator, m_Ey_z_lo, "Ey_z_lo(x,y,t)", "x", "y"); + m_set_E_z_hi |= ReadTangentialFieldParser(pp_insulator, m_Ex_z_hi, "Ex_z_hi(x,y,t)", "x", "y"); + m_set_E_z_hi |= ReadTangentialFieldParser(pp_insulator, m_Ey_z_hi, "Ey_z_hi(x,y,t)", "x", "y"); + +} + +void +PEC_Insulator::ApplyPEC_InsulatortoEfield ( + std::array Efield, + amrex::Array const & field_boundary_lo, + amrex::Array const & field_boundary_hi, + amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom, + int lev, PatchType patch_type, amrex::Vector const & ref_ratios, + amrex::Real time, + bool split_pml_field) +{ + bool const E_like = true; + ApplyPEC_InsulatortoField(Efield, field_boundary_lo, field_boundary_hi, ng_fieldgather, geom, + lev, patch_type, ref_ratios, time, split_pml_field, + E_like, +#if (AMREX_SPACEDIM > 1) + m_set_E_x_lo, m_set_E_x_hi, + m_Ey_x_lo, m_Ez_x_lo, m_Ey_x_hi, m_Ez_x_hi, +#endif +#if defined(WARPX_DIM_3D) + m_set_E_y_lo, m_set_E_y_hi, + m_Ex_y_lo, m_Ez_y_lo, m_Ex_y_hi, m_Ez_y_hi, +#endif + m_set_E_z_lo, m_set_E_z_hi, + m_Ex_z_lo, m_Ey_z_lo, m_Ex_z_hi, m_Ey_z_hi); +} + + +void +PEC_Insulator::ApplyPEC_InsulatortoBfield ( + std::array Bfield, + amrex::Array const & field_boundary_lo, + amrex::Array const & field_boundary_hi, + amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom, + int lev, PatchType patch_type, amrex::Vector const & ref_ratios, + amrex::Real time) +{ + bool const E_like = false; + bool const split_pml_field = false; + ApplyPEC_InsulatortoField(Bfield, field_boundary_lo, field_boundary_hi, ng_fieldgather, geom, + lev, patch_type, ref_ratios, time, split_pml_field, + E_like, +#if (AMREX_SPACEDIM > 1) + m_set_B_x_lo, m_set_B_x_hi, + m_By_x_lo, m_Bz_x_lo, m_By_x_hi, m_Bz_x_hi, +#endif +#if defined(WARPX_DIM_3D) + m_set_B_y_lo, m_set_B_y_hi, + m_Bx_y_lo, m_Bz_y_lo, m_Bx_y_hi, m_Bz_y_hi, +#endif + m_set_B_z_lo, m_set_B_z_hi, + m_Bx_z_lo, m_By_z_lo, m_Bx_z_hi, m_By_z_hi); +} + + +void +PEC_Insulator::ApplyPEC_InsulatortoField ( + std::array field, + amrex::Array const & field_boundary_lo, + amrex::Array const & field_boundary_hi, + amrex::IntVect const & ng_fieldgather, amrex::Geometry const & geom, + int lev, PatchType patch_type, amrex::Vector const & ref_ratios, + amrex::Real time, + bool split_pml_field, + bool E_like, +#if (AMREX_SPACEDIM > 1) + bool set_F_x_lo, bool set_F_x_hi, + std::unique_ptr const & a_Fy_x_lo, std::unique_ptr const & a_Fz_x_lo, + std::unique_ptr const & a_Fy_x_hi, std::unique_ptr const & a_Fz_x_hi, +#endif +#if defined(WARPX_DIM_3D) + bool set_F_y_lo, bool set_F_y_hi, + std::unique_ptr const & a_Fx_y_lo, std::unique_ptr const & a_Fz_y_lo, + std::unique_ptr const & a_Fx_y_hi, std::unique_ptr const & a_Fz_y_hi, +#endif + bool set_F_z_lo, bool set_F_z_hi, + std::unique_ptr const & a_Fx_z_lo, std::unique_ptr const & a_Fy_z_lo, + std::unique_ptr const & a_Fx_z_hi, std::unique_ptr const & a_Fy_z_hi) +{ + using namespace amrex::literals; + amrex::Box domain_box = geom.Domain(); + if (patch_type == PatchType::coarse && (lev > 0)) { + domain_box.coarsen(ref_ratios[lev-1]); + } + amrex::IntVect const domain_lo = domain_box.smallEnd(); + amrex::IntVect const domain_hi = domain_box.bigEnd(); + amrex::GpuArray fbndry_lo; + amrex::GpuArray fbndry_hi; + for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { + fbndry_lo[idim] = field_boundary_lo[idim]; + fbndry_hi[idim] = field_boundary_hi[idim]; + } + +#if (AMREX_SPACEDIM > 1) + amrex::ParserExecutor<2> const area_parsers_x_lo = m_insulator_area_lo[0]->compile<2>(); + amrex::ParserExecutor<2> const area_parsers_x_hi = m_insulator_area_hi[0]->compile<2>(); +#endif +#if defined(WARPX_DIM_3D) + amrex::ParserExecutor<2> const area_parsers_y_lo = m_insulator_area_lo[1]->compile<2>(); + amrex::ParserExecutor<2> const area_parsers_y_hi = m_insulator_area_hi[1]->compile<2>(); +#endif + amrex::ParserExecutor<2> const area_parsers_z_lo = m_insulator_area_lo[WARPX_ZINDEX]->compile<2>(); + amrex::ParserExecutor<2> const area_parsers_z_hi = m_insulator_area_hi[WARPX_ZINDEX]->compile<2>(); + +#if (AMREX_SPACEDIM > 1) + amrex::ParserExecutor<3> const Fy_x_lo_parser = a_Fy_x_lo->compile<3>(); + amrex::ParserExecutor<3> const Fz_x_lo_parser = a_Fz_x_lo->compile<3>(); + amrex::ParserExecutor<3> const Fy_x_hi_parser = a_Fy_x_hi->compile<3>(); + amrex::ParserExecutor<3> const Fz_x_hi_parser = a_Fz_x_hi->compile<3>(); +#endif +#if defined(WARPX_DIM_3D) + amrex::ParserExecutor<3> const Fx_y_lo_parser = a_Fx_y_lo->compile<3>(); + amrex::ParserExecutor<3> const Fz_y_lo_parser = a_Fz_y_lo->compile<3>(); + amrex::ParserExecutor<3> const Fx_y_hi_parser = a_Fx_y_hi->compile<3>(); + amrex::ParserExecutor<3> const Fz_y_hi_parser = a_Fz_y_hi->compile<3>(); +#endif + amrex::ParserExecutor<3> const Fx_z_lo_parser = a_Fx_z_lo->compile<3>(); + amrex::ParserExecutor<3> const Fy_z_lo_parser = a_Fy_z_lo->compile<3>(); + amrex::ParserExecutor<3> const Fx_z_hi_parser = a_Fx_z_hi->compile<3>(); + amrex::ParserExecutor<3> const Fy_z_hi_parser = a_Fy_z_hi->compile<3>(); + + amrex::IntVect const Fx_nodal = field[0]->ixType().toIntVect(); + amrex::IntVect const Fy_nodal = field[1]->ixType().toIntVect(); + amrex::IntVect const Fz_nodal = field[2]->ixType().toIntVect(); + // For each field multifab, apply boundary condition to ncomponents + // If not split field, the boundary condition is applied to the regular field used in Maxwell's eq. + // If split_pml_field is true, then boundary condition is applied to all the split field components. + int const nComp_x = field[0]->nComp(); + int const nComp_y = field[1]->nComp(); + int const nComp_z = field[2]->nComp(); + + std::array const & dx = WarpX::CellSize(lev); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(*field[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + // Extract field data + amrex::Array4 const & Fx = field[0]->array(mfi); + amrex::Array4 const & Fy = field[1]->array(mfi); + amrex::Array4 const & Fz = field[2]->array(mfi); + + // Extract tileboxes for which to loop + // if split field, the box includes nodal flag + // For E-field used in Maxwell's update, nodal flag plus cells that particles + // gather fields from in the guard-cell region are included. + // Note that for simulations without particles or laser, ng_field_gather is 0 + // and the guard-cell values of the E-field multifab will not be modified. + amrex::Box const & tex = (split_pml_field) ? mfi.tilebox(field[0]->ixType().toIntVect()) + : mfi.tilebox(field[0]->ixType().toIntVect(), ng_fieldgather); + amrex::Box const & tey = (split_pml_field) ? mfi.tilebox(field[1]->ixType().toIntVect()) + : mfi.tilebox(field[1]->ixType().toIntVect(), ng_fieldgather); + amrex::Box const & tez = (split_pml_field) ? mfi.tilebox(field[2]->ixType().toIntVect()) + : mfi.tilebox(field[2]->ixType().toIntVect(), ng_fieldgather); + + const amrex::XDim3 xyzmin_x = WarpX::LowerCorner(tex, lev, 0._rt); + const amrex::XDim3 xyzmin_y = WarpX::LowerCorner(tey, lev, 0._rt); + const amrex::XDim3 xyzmin_z = WarpX::LowerCorner(tez, lev, 0._rt); + amrex::IntVect const lo_x = tex.smallEnd(); + amrex::IntVect const lo_y = tey.smallEnd(); + amrex::IntVect const lo_z = tez.smallEnd(); + + // loop over cells and update fields + amrex::ParallelFor( + tex, nComp_x, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + amrex::ignore_unused(j, k); + + amrex::IntVect const iv(AMREX_D_DECL(i, j, k)); + amrex::Real const x = (AMREX_SPACEDIM > 1 ? xyzmin_x.x + (iv[0] - lo_x[0])*dx[0] : 0._rt); + amrex::Real const y = (AMREX_SPACEDIM == 3 ? xyzmin_x.y + (iv[1] - lo_x[1])*dx[1] : 0._rt); +#if (AMREX_SPACEDIM > 1) + amrex::Real const z = xyzmin_x.z + (iv[WARPX_ZINDEX] - lo_x[WARPX_ZINDEX])*dx[2]; +#endif + + amrex::IntVect is_insulator_lo; + amrex::IntVect is_insulator_hi; + amrex::RealVect F_lo, F_hi; + amrex::IntVect set_field_lo; + amrex::IntVect set_field_hi; +#if (AMREX_SPACEDIM > 1) + is_insulator_lo[0] = (area_parsers_x_lo(y, z) > 0._rt); + is_insulator_hi[0] = (area_parsers_x_hi(y, z) > 0._rt); + F_lo[0] = 0._rt; // Will be unused + F_hi[0] = 0._rt; // Will be unused + set_field_lo[0] = 0; // Will be unused + set_field_hi[0] = 0; // Will be unused +#endif +#if defined(WARPX_DIM_3D) + is_insulator_lo[1] = (area_parsers_y_lo(x, z) > 0._rt); + is_insulator_hi[1] = (area_parsers_y_hi(x, z) > 0._rt); + F_lo[1] = (set_F_y_lo ? Fx_y_lo_parser(x, z, time) : 0._rt); + F_hi[1] = (set_F_y_hi ? Fx_y_hi_parser(x, z, time) : 0._rt); + set_field_lo[1] = set_F_y_lo; + set_field_hi[1] = set_F_y_hi; +#endif + is_insulator_lo[WARPX_ZINDEX] = (area_parsers_z_lo(x, y) > 0._rt); + is_insulator_hi[WARPX_ZINDEX] = (area_parsers_z_hi(x, y) > 0._rt); + F_lo[WARPX_ZINDEX] = (set_F_z_lo ? Fx_z_lo_parser(x, y, time) : 0._rt); + F_hi[WARPX_ZINDEX] = (set_F_z_hi ? Fx_z_hi_parser(x, y, time) : 0._rt); + set_field_lo[WARPX_ZINDEX] = set_F_z_lo; + set_field_hi[WARPX_ZINDEX] = set_F_z_hi; + + int const icomp = 0; + ::SetFieldOnPEC_Insulator(icomp, domain_lo, domain_hi, iv, n, + Fx, E_like, Fx_nodal, is_insulator_lo, is_insulator_hi, + F_lo, F_hi, set_field_lo, set_field_hi, + fbndry_lo, fbndry_hi); + }, + tey, nComp_y, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + amrex::ignore_unused(j, k); + + amrex::IntVect const iv(AMREX_D_DECL(i, j, k)); + amrex::Real const x = (AMREX_SPACEDIM > 1 ? xyzmin_y.x + (iv[0] - lo_y[0])*dx[0] : 0._rt); + amrex::Real const y = (AMREX_SPACEDIM == 3 ? xyzmin_y.y + (iv[1] - lo_y[1])*dx[1] : 0._rt); +#if (AMREX_SPACEDIM > 1) + amrex::Real const z = xyzmin_y.z + (iv[WARPX_ZINDEX] - lo_y[WARPX_ZINDEX])*dx[2]; +#endif + + amrex::IntVect is_insulator_lo; + amrex::IntVect is_insulator_hi; + amrex::RealVect F_lo, F_hi; + amrex::IntVect set_field_lo; + amrex::IntVect set_field_hi; +#if (AMREX_SPACEDIM > 1) + is_insulator_lo[0] = (area_parsers_x_lo(y, z) > 0._rt); + is_insulator_hi[0] = (area_parsers_x_hi(y, z) > 0._rt); + F_lo[0] = (set_F_x_lo ? Fy_x_lo_parser(y, z, time) : 0._rt); + F_hi[0] = (set_F_x_hi ? Fy_x_hi_parser(y, z, time) : 0._rt); + set_field_lo[0] = set_F_x_lo; + set_field_hi[0] = set_F_x_hi; +#endif +#if defined(WARPX_DIM_3D) + is_insulator_lo[1] = (area_parsers_y_lo(x, z) > 0._rt); + is_insulator_hi[1] = (area_parsers_y_hi(x, z) > 0._rt); + F_lo[1] = 0._rt; // Will be unused + F_hi[1] = 0._rt; // Will be unused + set_field_lo[1] = 0; // Will be unused + set_field_hi[1] = 0; // Will be unused +#endif + is_insulator_lo[WARPX_ZINDEX] = (area_parsers_z_lo(x, y) > 0._rt); + is_insulator_hi[WARPX_ZINDEX] = (area_parsers_z_hi(x, y) > 0._rt); + F_lo[WARPX_ZINDEX] = (set_F_z_lo ? Fy_z_lo_parser(x, y, time) : 0._rt); + F_hi[WARPX_ZINDEX] = (set_F_z_hi ? Fy_z_hi_parser(x, y, time) : 0._rt); + set_field_lo[WARPX_ZINDEX] = set_F_z_lo; + set_field_hi[WARPX_ZINDEX] = set_F_z_hi; + + int const icomp = 1; + ::SetFieldOnPEC_Insulator(icomp, domain_lo, domain_hi, iv, n, + Fy, E_like, Fy_nodal, is_insulator_lo, is_insulator_hi, + F_lo, F_hi, set_field_lo, set_field_hi, + fbndry_lo, fbndry_hi); + }, + tez, nComp_z, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + amrex::ignore_unused(j, k); + + amrex::IntVect const iv(AMREX_D_DECL(i, j, k)); + amrex::Real const x = (AMREX_SPACEDIM > 1 ? xyzmin_z.x + (iv[0] - lo_z[0])*dx[0] : 0._rt); + amrex::Real const y = (AMREX_SPACEDIM == 3 ? xyzmin_z.y + (iv[1] - lo_z[1])*dx[1] : 0._rt); +#if (AMREX_SPACEDIM > 1) + amrex::Real const z = xyzmin_z.z + (iv[WARPX_ZINDEX] - lo_z[WARPX_ZINDEX])*dx[2]; +#endif + + amrex::IntVect is_insulator_lo; + amrex::IntVect is_insulator_hi; + amrex::RealVect F_lo, F_hi; + amrex::IntVect set_field_lo; + amrex::IntVect set_field_hi; +#if (AMREX_SPACEDIM > 1) + is_insulator_lo[0] = (area_parsers_x_lo(y, z) > 0._rt); + is_insulator_hi[0] = (area_parsers_x_hi(y, z) > 0._rt); + F_lo[0] = (set_F_x_lo ? Fz_x_lo_parser(y, z, time) : 0._rt); + F_hi[0] = (set_F_x_hi ? Fz_x_hi_parser(y, z, time) : 0._rt); + set_field_lo[0] = set_F_x_lo; + set_field_hi[0] = set_F_x_hi; +#endif +#if defined(WARPX_DIM_3D) + is_insulator_lo[1] = (area_parsers_y_lo(x, z) > 0._rt); + is_insulator_hi[1] = (area_parsers_y_hi(x, z) > 0._rt); + F_lo[1] = (set_F_y_lo ? Fz_y_lo_parser(x, z, time) : 0._rt); + F_hi[1] = (set_F_y_hi ? Fz_y_hi_parser(x, z, time) : 0._rt); + set_field_lo[1] = set_F_y_lo; + set_field_hi[1] = set_F_y_hi; +#endif + is_insulator_lo[WARPX_ZINDEX] = (area_parsers_z_lo(x, y) > 0._rt); + is_insulator_hi[WARPX_ZINDEX] = (area_parsers_z_hi(x, y) > 0._rt); + F_lo[WARPX_ZINDEX] = 0._rt; // Will be unused + F_hi[WARPX_ZINDEX] = 0._rt; // Will be unused + set_field_lo[WARPX_ZINDEX] = 0; // Will be unused + set_field_hi[WARPX_ZINDEX] = 0; // Will be unused + + int const icomp = 2; + ::SetFieldOnPEC_Insulator(icomp, domain_lo, domain_hi, iv, n, + Fz, E_like, Fz_nodal, is_insulator_lo, is_insulator_hi, + F_lo, F_hi, set_field_lo, set_field_hi, + fbndry_lo, fbndry_hi); + } + ); + } +} diff --git a/Source/BoundaryConditions/PEC_Insulator_fwd.H b/Source/BoundaryConditions/PEC_Insulator_fwd.H new file mode 100644 index 00000000000..9b2c1b05307 --- /dev/null +++ b/Source/BoundaryConditions/PEC_Insulator_fwd.H @@ -0,0 +1,13 @@ +/* Copyright 2024 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef PEC_INSULATOR_FWD_H +#define PEC_INSULATOR_FWD_H + +class PEC_Insulator; + +#endif /* PEC_INSULATOR_FWD_H */ diff --git a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp index dc41e95f40f..7566979557e 100644 --- a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp +++ b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp @@ -1,4 +1,5 @@ #include "WarpX.H" +#include "BoundaryConditions/PEC_Insulator.H" #include "BoundaryConditions/PML.H" #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H" @@ -92,6 +93,47 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type) } } + if (::isAnyBoundary(field_boundary_lo, field_boundary_hi)) { + amrex::Real const tnew = gett_new(lev); + if (patch_type == PatchType::fine) { + pec_insulator_boundary->ApplyPEC_InsulatortoEfield( + {m_fields.get(FieldType::Efield_fp,Direction{0},lev), + m_fields.get(FieldType::Efield_fp,Direction{1},lev), + m_fields.get(FieldType::Efield_fp,Direction{2},lev)}, + field_boundary_lo, field_boundary_hi, + get_ng_fieldgather(), Geom(lev), + lev, patch_type, ref_ratio, tnew); + if (::isAnyBoundary(field_boundary_lo, field_boundary_hi)) { + // apply pec on split E-fields in PML region + const bool split_pml_field = true; + pec_insulator_boundary->ApplyPEC_InsulatortoEfield( + m_fields.get_alldirs(FieldType::pml_E_fp, lev), + field_boundary_lo, field_boundary_hi, + get_ng_fieldgather(), Geom(lev), + lev, patch_type, ref_ratio, tnew, + split_pml_field); + } + } else { + pec_insulator_boundary->ApplyPEC_InsulatortoEfield( + {m_fields.get(FieldType::Efield_cp,Direction{0},lev), + m_fields.get(FieldType::Efield_cp,Direction{1},lev), + m_fields.get(FieldType::Efield_cp,Direction{2},lev)}, + field_boundary_lo, field_boundary_hi, + get_ng_fieldgather(), Geom(lev), + lev, patch_type, ref_ratio, tnew); + if (::isAnyBoundary(field_boundary_lo, field_boundary_hi)) { + // apply pec on split E-fields in PML region + const bool split_pml_field = true; + pec_insulator_boundary->ApplyPEC_InsulatortoEfield( + m_fields.get_alldirs(FieldType::pml_E_cp, lev), + field_boundary_lo, field_boundary_hi, + get_ng_fieldgather(), Geom(lev), + lev, patch_type, ref_ratio, tnew, + split_pml_field); + } + } + } + #ifdef WARPX_DIM_RZ if (patch_type == PatchType::fine) { ApplyFieldBoundaryOnAxis(m_fields.get(FieldType::Efield_fp, Direction{0}, lev), @@ -129,6 +171,27 @@ void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_d } } + if (::isAnyBoundary(field_boundary_lo, field_boundary_hi)) { + amrex::Real const tnew = gett_new(lev); + if (patch_type == PatchType::fine) { + pec_insulator_boundary->ApplyPEC_InsulatortoBfield( + {m_fields.get(FieldType::Bfield_fp,Direction{0},lev), + m_fields.get(FieldType::Bfield_fp,Direction{1},lev), + m_fields.get(FieldType::Bfield_fp,Direction{2},lev)}, + field_boundary_lo, field_boundary_hi, + get_ng_fieldgather(), Geom(lev), + lev, patch_type, ref_ratio, tnew); + } else { + pec_insulator_boundary->ApplyPEC_InsulatortoBfield( + {m_fields.get(FieldType::Bfield_cp,Direction{0},lev), + m_fields.get(FieldType::Bfield_cp,Direction{1},lev), + m_fields.get(FieldType::Bfield_cp,Direction{2},lev)}, + field_boundary_lo, field_boundary_hi, + get_ng_fieldgather(), Geom(lev), + lev, patch_type, ref_ratio, tnew); + } + } + // Silver-Mueller boundaries are only applied on the first half-push of B // This is because the formula used for Silver-Mueller assumes that // E and B are staggered in time, which is only true after the first half-push diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp index 1ce7959e7e4..932304d5009 100644 --- a/Source/Python/WarpX.cpp +++ b/Source/Python/WarpX.cpp @@ -7,6 +7,7 @@ #include // see WarpX.cpp - full includes for _fwd.H headers +#include #include #include #include diff --git a/Source/Utils/Parser/ParserUtils.H b/Source/Utils/Parser/ParserUtils.H index 96937abdbbe..19f976c3a6c 100644 --- a/Source/Utils/Parser/ParserUtils.H +++ b/Source/Utils/Parser/ParserUtils.H @@ -88,6 +88,20 @@ namespace utils::parser std::string& stored_string); + /** + * \brief If the input is provided, parse the string (typically a mathematical expression) from the + * input file and store it into a variable, replacing its contents. + * + * \param pp used to read the query_string `pp.=string` + * \param query_string ParmParse.query will look for this string + * \param stored_string variable in which the string to parse is stored + */ + bool Query_parserString( + amrex::ParmParse const& pp, + std::string const& query_string, + std::string& stored_string); + + /** Parse a string and return as a double precision floating point number * * In case the string cannot be interpreted as a double, diff --git a/Source/Utils/Parser/ParserUtils.cpp b/Source/Utils/Parser/ParserUtils.cpp index 0339b766e38..d017a6e019c 100644 --- a/Source/Utils/Parser/ParserUtils.cpp +++ b/Source/Utils/Parser/ParserUtils.cpp @@ -51,6 +51,19 @@ void utils::parser::Store_parserString( } } +bool utils::parser::Query_parserString( + amrex::ParmParse const& pp, + std::string const& query_string, + std::string& stored_string) +{ + bool const input_specified = pp.contains(query_string.c_str()); + if (input_specified) { + stored_string.clear(); + utils::parser::Store_parserString(pp, query_string, stored_string); + } + return input_specified; +} + int utils::parser::query (const amrex::ParmParse& a_pp, std::string const& group, char const * str, std::string& val) { const bool is_specified_without_group = a_pp.contains(str); diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index f67aeddadd0..088ef295364 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -132,6 +132,7 @@ AMREX_ENUM(FieldBoundaryType, Open, // Used in the Integrated Green Function Poisson solver // Note that the solver implicitely assumes open BCs: // no need to enforce them separately + PECInsulator, // Mixed boundary with PEC and insulator Default = PML); /** Particle boundary conditions at the domain boundary diff --git a/Source/WarpX.H b/Source/WarpX.H index a635196d044..4dc3ab8c8be 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -12,6 +12,7 @@ #ifndef WARPX_H_ #define WARPX_H_ +#include "BoundaryConditions/PEC_Insulator_fwd.H" #include "BoundaryConditions/PML_fwd.H" #include "Diagnostics/MultiDiagnostics_fwd.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags_fwd.H" @@ -1435,6 +1436,9 @@ private: #endif amrex::Real v_particle_pml; + // Insulator boundary conditions + std::unique_ptr pec_insulator_boundary; + // External fields parameters std::unique_ptr m_p_ext_field_params; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d1e3108e32a..cb46f3129c8 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -11,6 +11,7 @@ */ #include "WarpX.H" +#include "BoundaryConditions/PEC_Insulator.H" #include "BoundaryConditions/PML.H" #include "Diagnostics/MultiDiagnostics.H" #include "Diagnostics/ReducedDiags/MultiReducedDiags.H" @@ -1685,6 +1686,9 @@ WarpX::ReadParameters () } } + // Setup pec_insulator boundary conditions + pec_insulator_boundary = std::make_unique(); + // for slice generation // { const ParmParse pp_slice("slice");