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
    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 componentes 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)
            // 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)
                                             ? false : true );
#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) ? false : true );
#else
                const bool is_tangent_to_PEC = ( ( icomp == idim) ? false : true );
#endif
                if (isPECBoundary == true) {
                    // guard cell outside the domain by "ig" number cells, in direction, idim
                    const int ig = ( (iside == 0)
                                 ? (dom_lo[idim] - ijk_vec[idim])
                                 : (ijk_vec[idim] - (dom_hi[idim] + is_nodal[idim]) ) );
                    if (ig == 0) {
                        if (is_tangent_to_PEC == true and 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)
                                        : (dom_hi[idim] + is_nodal[idim] - ig));
                        GuardCell = true;
                        // tangential components are inverted across PEC boundary
                        if (is_tangent_to_PEC == true) sign *= -1._rt;
                    }
                } // is PEC boundary
            } // loop over iside
        } // loop over dimensions
        if (OnPECBoundary == true) {
            // if ijk_vec is on a PEC boundary in any direction, set Etangential to 0.
            Efield(ijk_vec,n) = 0._rt;
        } else if (GuardCell == true) {
            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 == true) {
#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)
                                                ? true : false );
#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) ? true : false );
#else
                    const bool is_normal_to_PEC = ( ( icomp == idim) ? true : false );
#endif
                    // guard cell outside the domain by "ig" number cells, in direction, idim
                    const int ig = ( (iside == 0)
                                   ? (dom_lo[idim] - ijk_vec[idim])
                                   : (ijk_vec[idim] - (dom_hi[idim] + is_nodal[idim]) ) );
                    if (ig == 0) {
                        // Only normal component is set to 0
                        if (is_normal_to_PEC == true and 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)
                                       : (dom_hi[idim] + is_nodal[idim] - ig));
                        GuardCell = true;
                        // Sign of the normal component in guard cell is inverted
                        if (is_normal_to_PEC == true) sign *= -1._rt;
                    }
                } // if PEC Boundary
            } // loop over sides
        } // loop of dimensions

        if (OnPECBoundary == true) {
            // if ijk_vec is on a PEC boundary in any direction, set Bnormal to 0.
            Bfield(ijk_vec,n) = 0._rt;
        } else if (GuardCell == true) {
            // 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);
        }
    }

    bool isAnyBoundaryPEC();
    void ApplyPECtoEfield ( std::array<amrex::MultiFab*, 3> Efield,
                            const int lev, PatchType patch_type,
                            const bool split_pml_field = false);
    void ApplyPECtoBfield ( std::array<amrex::MultiFab*, 3> Bfield,
                            const int lev, PatchType patch_type);
}

#endif // WarpX_PEC_KERNELS_H_