WarpX_PEC.H

Go to the documentation of this file.

#ifndef WARPX_PEC_KERNELS_H_
#define WARPX_PEC_KERNELS_H_
 
#include "WarpX.H"
#include "Utils/WarpXAlgorithmSelection.H"
 
#include <AMReX_Array.H>
#include <AMReX_Array4.H>
#include <AMReX_Config.H>
#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_IntVect.H>
#include <AMReX_REAL.H>
 
#include <AMReX_BaseFwd.H>
 
#include <array>
#include <memory>
 
namespace PEC {
using namespace amrex;
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    bool is_boundary_PEC (amrex::GpuArray<int, 3> const& fboundary, int dir) {
        return ( fboundary[dir] == FieldBoundaryType::PEC );
    }
 
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    bool is_boundary_reflecting (amrex::GpuArray<ParticleBoundaryType, 3> const& pboundary, int dir) {
        return ( pboundary[dir] == ParticleBoundaryType::Reflecting );
    }
 
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    int get_cell_count_to_boundary (const amrex::IntVect& dom_lo,
        const amrex::IntVect& dom_hi, const amrex::IntVect& ijk_vec,
        const amrex::IntVect& is_nodal, const int idim, const int iside)
    {
        return ((iside == 0) ? (dom_lo[idim] - ijk_vec[idim])
                             : (ijk_vec[idim] - (dom_hi[idim] + is_nodal[idim])));
    }
 
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void SetEfieldOnPEC (const int icomp, const amrex::IntVect & dom_lo,
                                const amrex::IntVect &dom_hi,
                                const amrex::IntVect &ijk_vec, const int n,
                                amrex::Array4<amrex::Real> const& Efield,
                                const amrex::IntVect& is_nodal,
                                amrex::GpuArray<int, 3> const& fbndry_lo,
                                amrex::GpuArray<int, 3> const& fbndry_hi )
    {
        // Tangential Efield components in guard cells set equal and opposite to cells
        // in the mirror locations across the PEC boundary, whereas normal E-field
        // components are set equal to values in the mirror locations across the PEC
        // boundary. Here we just initialize it.
        amrex::IntVect ijk_mirror = ijk_vec;
        bool OnPECBoundary = false;
        bool GuardCell = false;
        amrex::Real sign = 1._rt;
        // Loop over all the dimensions
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
            // Loop over sides, iside = 0 (lo), iside = 1 (hi)
            for (int iside = 0; iside < 2; ++iside) {
                const bool isPECBoundary = ( (iside == 0)
                                        ? is_boundary_PEC(fbndry_lo, idim)
                                        : is_boundary_PEC(fbndry_hi, idim) );
#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
                // For 2D : for icomp==1, (Ey in XZ, Etheta in RZ),
                //          icomp=1 is tangential to both x and z boundaries
                //          The logic below ensures that the flags are set right for 2D
                const bool is_tangent_to_PEC = (icomp != AMREX_SPACEDIM*idim);
#elif (defined WARPX_DIM_1D_Z)
                // For 1D : icomp=0 and icomp=1 (Ex and Ey are tangential to the z boundary)
                //          The logic below ensures that the flags are set right for 1D
                const bool is_tangent_to_PEC = (icomp != idim+2);
#else
                const bool is_tangent_to_PEC = (icomp != idim);
#endif
                if (isPECBoundary) {
                    // grid point ijk_vec is ig number of points pass the
                    // domain boundary in direction, idim
                    const int ig = get_cell_count_to_boundary(
                        dom_lo, dom_hi, ijk_vec, is_nodal, idim, iside);
 
                    if (ig == 0) {
                        if (is_tangent_to_PEC && is_nodal[idim] == 1) {
                            OnPECBoundary = true;
                        }
                    } else if (ig > 0) {
                        // Find mirror location across PEC boundary
                        ijk_mirror[idim] = ( ( iside == 0)
                                        ? (dom_lo[idim] + ig - (1 - is_nodal[idim]))
                                        : (dom_hi[idim] + 1 - ig));
                        GuardCell = true;
                        // tangential components are inverted across PEC boundary
                        if (is_tangent_to_PEC) { sign *= -1._rt; }
#if (defined WARPX_DIM_RZ)
                        if (icomp == 0 && idim == 0 && iside == 1) {
                            // Add radial scale so that drEr/dr = 0.
                            // This only works for the first guard cell and with
                            // Er cell centered in r.
                            const amrex::Real rguard = ijk_vec[idim] + 0.5_rt*(1._rt - is_nodal[idim]);
                            const amrex::Real rmirror = ijk_mirror[idim] + 0.5_rt*(1._rt - is_nodal[idim]);
                            sign *= rmirror/rguard;
                        }
#endif
                    }
                } // is PEC boundary
            } // loop over iside
        } // loop over dimensions
        if (OnPECBoundary) {
            // if ijk_vec is on a PEC boundary in any direction, set Etangential to 0.
            Efield(ijk_vec,n) = 0._rt;
        } else if (GuardCell) {
            Efield(ijk_vec,n) = sign * Efield(ijk_mirror,n);
        }
    }
 
 
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void SetBfieldOnPEC (const int icomp, const amrex::IntVect & dom_lo,
                           const amrex::IntVect & dom_hi,
                           const amrex::IntVect & ijk_vec, const int n,
                           amrex::Array4<amrex::Real> const& Bfield,
                           const amrex::IntVect & is_nodal,
                           amrex::GpuArray<int, 3> const& fbndry_lo,
                           amrex::GpuArray<int, 3> const& fbndry_hi )
    {
        amrex::IntVect ijk_mirror = ijk_vec;
        bool OnPECBoundary = false;
        bool GuardCell = false;
        amrex::Real sign = 1._rt;
        // Loop over all dimensions
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
            // Loop over sides, iside = 0 (lo), iside = 1 (hi)
            for (int iside = 0; iside < 2; ++iside) {
                const bool isPECBoundary = ( (iside == 0 )
                                        ? is_boundary_PEC(fbndry_lo, idim)
                                        : is_boundary_PEC(fbndry_hi, idim) );
                if (isPECBoundary) {
#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ)
                    // For 2D : for icomp==1, (By in XZ, Btheta in RZ),
                    //          icomp=1 is not normal to x or z boundary
                    //          The logic below ensures that the flags are set right for 2D
                    const bool is_normal_to_PEC = (icomp == (AMREX_SPACEDIM*idim));
#elif (defined WARPX_DIM_1D_Z)
                    // For 1D : icomp=0 and icomp=1 (Bx and By are not normal to the z boundary)
                    //          The logic below ensures that the flags are set right for 1D
                    const bool is_normal_to_PEC = (icomp == (idim+2));
#else
                    const bool is_normal_to_PEC = (icomp == idim);
#endif
 
                    // grid point ijk_vec is ig number of points pass the
                    // domain boundary in direction, idim
                    const int ig = get_cell_count_to_boundary(
                        dom_lo, dom_hi, ijk_vec, is_nodal, idim, iside);
 
                    if (ig == 0) {
                        // Only normal component is set to 0
                        if (is_normal_to_PEC && is_nodal[idim]==1) {
                            OnPECBoundary = true;
                        }
                    } else if ( ig > 0) {
                        // Mirror location inside the domain by "ig" number of cells
                        // across PEC boundary in direction, idim, and side, iside
                        ijk_mirror[idim] = ( (iside == 0)
                                        ? (dom_lo[idim] + ig - (1 - is_nodal[idim]))
                                        : (dom_hi[idim] + 1 - ig));
                        GuardCell = true;
                        // Sign of the normal component in guard cell is inverted
                        if (is_normal_to_PEC) { sign *= -1._rt; }
#if (defined WARPX_DIM_RZ)
                        if (icomp == 0 && idim == 0 && iside == 1) {
                            // Add radial scale so that drBr/dr = 0.
                            const amrex::Real rguard = ijk_vec[idim] + 0.5_rt*(1._rt - is_nodal[idim]);
                            const amrex::Real rmirror = ijk_mirror[idim] + 0.5_rt*(1._rt - is_nodal[idim]);
                            sign *= rmirror/rguard;
                        }
#endif
                    }
                } // if PEC Boundary
            } // loop over sides
        } // loop of dimensions
 
        if (OnPECBoundary) {
            // if ijk_vec is on a PEC boundary in any direction, set Bnormal to 0.
            Bfield(ijk_vec,n) = 0._rt;
        } else if (GuardCell) {
            // Bnormal and Btangential is set opposite and equal to the value
            // in the mirror location, respectively.
            Bfield(ijk_vec,n) = sign * Bfield(ijk_mirror,n);
        }
    }
 
 
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void SetRhoOrJfieldFromPEC (const int n,
                                const amrex::IntVect & ijk_vec,
                                amrex::Array4<amrex::Real> const& field,
                                amrex::GpuArray<GpuArray<int, 2>, AMREX_SPACEDIM> const& mirrorfac,
                                amrex::GpuArray<GpuArray<amrex::Real, 2>, AMREX_SPACEDIM> const& psign,
                                amrex::GpuArray<GpuArray<bool, 2>, AMREX_SPACEDIM> const& is_pec,
                                amrex::GpuArray<bool, AMREX_SPACEDIM> const& tangent_to_bndy,
                                amrex::Box const& fabbox)
    {
        // The boundary is handled in 2 steps:
        // 1) The cells internal to the domain are updated using the
        //    current deposited in the guard cells
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
        {
            for (int iside = 0; iside < 2; ++iside)
            {
                if (!is_pec[idim][iside]) { continue; }
 
                // Get the mirror guard cell index
                amrex::IntVect iv_mirror = ijk_vec;
                iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim];
 
                // On the PEC boundary the charge/current density is set to 0
                if (ijk_vec == iv_mirror) {
                    field(ijk_vec, n) = 0._rt;
                // otherwise update the internal cell if the mirror guard cell exists
                } else if (fabbox.contains(iv_mirror)) {
                    field(ijk_vec,n) += psign[idim][iside] * field(iv_mirror,n);
                }
            }
        }
        // 2) The guard cells are updated with the appropriate image
        //    charge based on the charge/current in the valid cells
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
        {
            for (int iside = 0; iside < 2; ++iside)
            {
                if (!is_pec[idim][iside]) { continue; }
 
                amrex::IntVect iv_mirror = ijk_vec;
                iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim];
                if (ijk_vec != iv_mirror && fabbox.contains(iv_mirror))
                {
                    if (tangent_to_bndy[idim]) {
                        field(iv_mirror, n) = -field(ijk_vec, n);
                    } else {
                        field(iv_mirror, n) = field(ijk_vec, n);
                    }
                }
            }
        }
    }
 
 
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void SetNeumannOnPEC (const int n,
                          const amrex::IntVect & ijk_vec,
                          amrex::Array4<amrex::Real> const& field,
                          amrex::GpuArray<GpuArray<int, 2>, AMREX_SPACEDIM> const& mirrorfac,
                          amrex::GpuArray<GpuArray<bool, 2>, AMREX_SPACEDIM> const& is_pec,
                          amrex::Box const& fabbox )
    {
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
        {
            for (int iside = 0; iside < 2; ++iside)
            {
                if (!is_pec[idim][iside]) { continue; }
 
                // Get the mirror guard cell index
                amrex::IntVect iv_mirror = ijk_vec;
                iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim];
 
                // On the PEC boundary the field value is set equal to the
                // first value in the domain (nodal fields)
                if (ijk_vec == iv_mirror) {
                    iv_mirror[idim] += (iside == 0) ? 1 : -1;
                    if (fabbox.contains(iv_mirror)) { field(ijk_vec, n) = field(iv_mirror, n); }
                }
                // otherwise set the mirror guard cell equal to the internal cell value
                else if (fabbox.contains(iv_mirror))
                {
                    field(iv_mirror, n) = field(ijk_vec, n);
                }
            }
        }
    }
 
 
    bool isAnyBoundaryPEC();
    void ApplyPECtoEfield ( std::array<amrex::MultiFab*, 3> Efield,
                            int lev, PatchType patch_type,
                            bool split_pml_field = false);
    void ApplyPECtoBfield ( std::array<amrex::MultiFab*, 3> Bfield,
                            int lev, PatchType patch_type);
 
    void ApplyPECtoRhofield(amrex::MultiFab* rho, int lev,
                            PatchType patch_type);
 
    void ApplyPECtoJfield(amrex::MultiFab* Jx, amrex::MultiFab* Jy,
                          amrex::MultiFab* Jz, int lev,
                          PatchType patch_type);
 
    void ApplyPECtoElectronPressure (amrex::MultiFab* Pefield,
                                     int lev, PatchType patch_type);
}
#endif // WarpX_PEC_KERNELS_H_