aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H
diff options
context:
space:
mode:
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H')
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H180
1 files changed, 93 insertions, 87 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H
index fa5dd073d..6a6bc6d4c 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H
@@ -12,6 +12,9 @@
#include <AMReX_Array4.H>
#include <AMReX_Gpu.H>
+#include <algorithm>
+#include <array>
+
/**
* This struct contains only static functions to initialize the stencil coefficients
* and to compute finite-difference derivatives for the Cartesian CKC algorithm.
@@ -28,42 +31,42 @@ struct CartesianCKCAlgorithm {
// Compute Cole-Karkkainen-Cowan coefficients according
// to Cowan - PRST-AB 16, 041303 (2013)
- Real const inv_dx = 1./cell_size[0];
- Real const inv_dy = 1./cell_size[1];
- Real const inv_dz = 1./cell_size[2];
+ Real const inv_dx = 1._rt/cell_size[0];
+ Real const inv_dy = 1._rt/cell_size[1];
+ Real const inv_dz = 1._rt/cell_size[2];
#if defined WARPX_DIM_3D
- Real const delta = std::max( { inv_dx,inv_dy,inv_dz } );
- Real const rx = (inv_dx/delta)*(inv_dx/delta);
- Real const ry = (inv_dy/delta)*(inv_dy/delta);
- Real const rz = (inv_dz/delta)*(inv_dz/delta);
- Real const beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry));
- Real const betaxy = ry*beta*inv_dx;
- Real const betaxz = rz*beta*inv_dx;
- Real const betayx = rx*beta*inv_dy;
- Real const betayz = rz*beta*inv_dy;
- Real const betazx = rx*beta*inv_dz;
- Real const betazy = ry*beta*inv_dz;
- Real const gamma = (0.0625 - 0.125*ry*rz/(ry*rz + rz*rx + rx*ry));
- Real const gammax = ry*rz*gamma;
- Real const gammay = rx*rz*gamma;
- Real const gammaz = rx*ry*gamma;
- Real const alphax = (1. - 2.*ry*beta - 2.*rz*beta - 4.*ry*rz*gamma)*inv_dx;
- Real const alphay = (1. - 2.*rx*beta - 2.*rz*beta - 4.*rx*rz*gamma)*inv_dy;
- Real const alphaz = (1. - 2.*rx*beta - 2.*ry*beta - 4.*rx*ry*gamma)*inv_dz;
+ Real const delta = std::max( { inv_dx,inv_dy,inv_dz } );
+ Real const rx = (inv_dx/delta)*(inv_dx/delta);
+ Real const ry = (inv_dy/delta)*(inv_dy/delta);
+ Real const rz = (inv_dz/delta)*(inv_dz/delta);
+ Real const beta = 0.125_rt*(1._rt - rx*ry*rz/(ry*rz + rz*rx + rx*ry));
+ Real const betaxy = ry*beta*inv_dx;
+ Real const betaxz = rz*beta*inv_dx;
+ Real const betayx = rx*beta*inv_dy;
+ Real const betayz = rz*beta*inv_dy;
+ Real const betazx = rx*beta*inv_dz;
+ Real const betazy = ry*beta*inv_dz;
+ Real const gamma = (0.0625_rt - 0.125_rt*ry*rz/(ry*rz + rz*rx + rx*ry));
+ Real const gammax = ry*rz*gamma;
+ Real const gammay = rx*rz*gamma;
+ Real const gammaz = rx*ry*gamma;
+ Real const alphax = (1._rt - 2._rt*ry*beta - 2._rt*rz*beta - 4._rt*ry*rz*gamma)*inv_dx;
+ Real const alphay = (1._rt - 2._rt*rx*beta - 2._rt*rz*beta - 4._rt*rx*rz*gamma)*inv_dy;
+ Real const alphaz = (1._rt - 2._rt*rx*beta - 2._rt*ry*beta - 4._rt*rx*ry*gamma)*inv_dz;
#elif defined WARPX_DIM_XZ
- Real const delta = std::max(inv_dx,inv_dz);
- Real const rx = (inv_dx/delta)*(inv_dx/delta);
- Real const rz = (inv_dz/delta)*(inv_dz/delta);
- Real const beta = 0.125;
- Real const betaxz = beta*rz*inv_dx;
- Real const betazx = beta*rx*inv_dz;
- Real const alphax = (1. - 2.*rz*beta)*inv_dx;
- Real const alphaz = (1. - 2.*rx*beta)*inv_dz;
- // Other coefficients are 0 in 2D Cartesian
- // (and will actually not be used in the stencil)
- Real const gammax=0, gammay=0, gammaz=0;
- Real const betaxy=0, betazy=0, betayx=0, betayz=0;
- Real const alphay=0;
+ Real const delta = std::max(inv_dx,inv_dz);
+ Real const rx = (inv_dx/delta)*(inv_dx/delta);
+ Real const rz = (inv_dz/delta)*(inv_dz/delta);
+ constexpr Real beta = 0.125_rt;
+ Real const betaxz = beta*rz*inv_dx;
+ Real const betazx = beta*rx*inv_dz;
+ Real const alphax = (1._rt - 2._rt*rz*beta)*inv_dx;
+ Real const alphaz = (1._rt - 2._rt*rx*beta)*inv_dz;
+ // Other coefficients are 0 in 2D Cartesian
+ // (and will actually not be used in the stencil)
+ constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt;
+ constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt;
+ constexpr Real alphay=0._rt;
#endif
// Store the coefficients in array `stencil_coefs`, in prescribed order
@@ -88,7 +91,7 @@ struct CartesianCKCAlgorithm {
}
/**
- /* Perform derivative along x on a cell-centered grid, from a nodal field `F`*/
+ * Perform derivative along x on a cell-centered grid, from a nodal field `F` */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real UpwardDx (
amrex::Array4<amrex::Real> const& F,
@@ -100,24 +103,24 @@ struct CartesianCKCAlgorithm {
amrex::Real const betaxz = coefs_x[3];
amrex::Real const gammax = coefs_x[4];
#if defined WARPX_DIM_3D
- return alphax * (F(i+1,j ,k ) - F(i, j, k ))
- + betaxy * (F(i+1,j+1,k ) - F(i ,j+1,k )
- + F(i+1,j-1,k ) - F(i ,j-1,k ))
- + betaxz * (F(i+1,j ,k+1) - F(i ,j ,k+1)
- + F(i+1,j ,k-1) - F(i ,j ,k-1))
- + gammax * (F(i+1,j+1,k+1) - F(i ,j+1,k+1)
- + F(i+1,j-1,k+1) - F(i ,j-1,k+1)
- + F(i+1,j+1,k-1) - F(i ,j+1,k-1)
- + F(i+1,j-1,k-1) - F(i ,j-1,k-1));
+ return alphax * (F(i+1,j ,k ) - F(i, j, k ))
+ + betaxy * (F(i+1,j+1,k ) - F(i ,j+1,k )
+ + F(i+1,j-1,k ) - F(i ,j-1,k ))
+ + betaxz * (F(i+1,j ,k+1) - F(i ,j ,k+1)
+ + F(i+1,j ,k-1) - F(i ,j ,k-1))
+ + gammax * (F(i+1,j+1,k+1) - F(i ,j+1,k+1)
+ + F(i+1,j-1,k+1) - F(i ,j-1,k+1)
+ + F(i+1,j+1,k-1) - F(i ,j+1,k-1)
+ + F(i+1,j-1,k-1) - F(i ,j-1,k-1));
#elif (defined WARPX_DIM_XZ)
- return alphax * (F(i+1,j ,k ) - F(i, j, k ))
- + betaxz * (F(i+1,j+1,k ) - F(i ,j+1,k )
- + F(i+1,j-1,k ) - F(i ,j-1,k ));
+ return alphax * (F(i+1,j ,k ) - F(i, j, k ))
+ + betaxz * (F(i+1,j+1,k ) - F(i ,j+1,k )
+ + F(i+1,j-1,k ) - F(i ,j-1,k ));
#endif
};
/**
- /* Perform derivative along x on a nodal grid, from a cell-centered field `F`*/
+ * Perform derivative along x on a nodal grid, from a cell-centered field `F` */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DownwardDx (
amrex::Array4<amrex::Real> const& F,
@@ -129,79 +132,82 @@ struct CartesianCKCAlgorithm {
};
/**
- /* Perform derivative along y on a cell-centered grid, from a nodal field `F`*/
+ * Perform derivative along y on a cell-centered grid, from a nodal field `F` */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real UpwardDy (
amrex::Array4<amrex::Real> const& F,
amrex::Real const * const coefs_y, int const n_coefs_y,
int const i, int const j, int const k ) {
+ using namespace amrex;
#if defined WARPX_DIM_3D
- amrex::Real const alphay = coefs_y[1];
- amrex::Real const betayz = coefs_y[2];
- amrex::Real const betayx = coefs_y[3];
- amrex::Real const gammay = coefs_y[4];
- return alphay * (F(i ,j+1,k ) - F(i ,j ,k ))
- + betayx * (F(i+1,j+1,k ) - F(i+1,j ,k )
- + F(i-1,j+1,k ) - F(i-1,j ,k ))
- + betayz * (F(i ,j+1,k+1) - F(i ,j ,k+1)
- + F(i ,j+1,k-1) - F(i ,j ,k-1))
- + gammay * (F(i+1,j+1,k+1) - F(i+1,j ,k+1)
- + F(i-1,j+1,k+1) - F(i-1,j ,k+1)
- + F(i+1,j+1,k-1) - F(i+1,j ,k-1)
- + F(i-1,j+1,k-1) - F(i-1,j ,k-1));
+ Real const alphay = coefs_y[1];
+ Real const betayz = coefs_y[2];
+ Real const betayx = coefs_y[3];
+ Real const gammay = coefs_y[4];
+ return alphay * (F(i ,j+1,k ) - F(i ,j ,k ))
+ + betayx * (F(i+1,j+1,k ) - F(i+1,j ,k )
+ + F(i-1,j+1,k ) - F(i-1,j ,k ))
+ + betayz * (F(i ,j+1,k+1) - F(i ,j ,k+1)
+ + F(i ,j+1,k-1) - F(i ,j ,k-1))
+ + gammay * (F(i+1,j+1,k+1) - F(i+1,j ,k+1)
+ + F(i-1,j+1,k+1) - F(i-1,j ,k+1)
+ + F(i+1,j+1,k-1) - F(i+1,j ,k-1)
+ + F(i-1,j+1,k-1) - F(i-1,j ,k-1));
#elif (defined WARPX_DIM_XZ)
- return 0; // 2D Cartesian: derivative along y is 0
+ return 0._rt; // 2D Cartesian: derivative along y is 0
#endif
};
/**
- /* Perform derivative along y on a nodal grid, from a cell-centered field `F`*/
+ * Perform derivative along y on a nodal grid, from a cell-centered field `F` */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DownwardDy (
amrex::Array4<amrex::Real> const& F,
amrex::Real const * const coefs_y, int const n_coefs_y,
int const i, int const j, int const k ) {
+ using namespace amrex;
#if defined WARPX_DIM_3D
- amrex::Real const inv_dy = coefs_y[0];
- return inv_dy*( F(i,j,k) - F(i,j-1,k) );
+ Real const inv_dy = coefs_y[0];
+ return inv_dy*( F(i,j,k) - F(i,j-1,k) );
#elif (defined WARPX_DIM_XZ)
- return 0; // 2D Cartesian: derivative along y is 0
+ return 0._rt; // 2D Cartesian: derivative along y is 0
#endif
};
/**
- /* Perform derivative along z on a cell-centered grid, from a nodal field `F`*/
+ * Perform derivative along z on a cell-centered grid, from a nodal field `F` */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real UpwardDz (
amrex::Array4<amrex::Real> const& F,
amrex::Real const * const coefs_z, int const n_coefs_z,
int const i, int const j, int const k ) {
- amrex::Real const alphaz = coefs_z[1];
- amrex::Real const betazx = coefs_z[2];
- amrex::Real const betazy = coefs_z[3];
- amrex::Real const gammaz = coefs_z[4];
+ using namespace amrex;
+ Real const alphaz = coefs_z[1];
+ Real const betazx = coefs_z[2];
+ Real const betazy = coefs_z[3];
+ Real const gammaz = coefs_z[4];
#if defined WARPX_DIM_3D
- return alphaz * (F(i ,j ,k+1) - F(i ,j ,k ))
- + betazx * (F(i+1,j ,k+1) - F(i+1,j ,k )
- + F(i-1,j ,k+1) - F(i-1,j ,k ))
- + betazy * (F(i ,j+1,k+1) - F(i ,j+1,k )
- + F(i ,j-1,k+1) - F(i ,j-1,k ))
- + gammaz * (F(i+1,j+1,k+1) - F(i+1,j+1,k )
- + F(i-1,j+1,k+1) - F(i-1,j+1,k )
- + F(i+1,j-1,k+1) - F(i+1,j-1,k )
- + F(i-1,j-1,k+1) - F(i-1,j-1,k ));
+ return alphaz * (F(i ,j ,k+1) - F(i ,j ,k ))
+ + betazx * (F(i+1,j ,k+1) - F(i+1,j ,k )
+ + F(i-1,j ,k+1) - F(i-1,j ,k ))
+ + betazy * (F(i ,j+1,k+1) - F(i ,j+1,k )
+ + F(i ,j-1,k+1) - F(i ,j-1,k ))
+ + gammaz * (F(i+1,j+1,k+1) - F(i+1,j+1,k )
+ + F(i-1,j+1,k+1) - F(i-1,j+1,k )
+ + F(i+1,j-1,k+1) - F(i+1,j-1,k )
+ + F(i-1,j-1,k+1) - F(i-1,j-1,k ));
#elif (defined WARPX_DIM_XZ)
- return alphaz * (F(i ,j+1,k ) - F(i ,j ,k ))
- + betazx * (F(i+1,j+1,k ) - F(i+1,j ,k )
- + F(i-1,j+1,k ) - F(i-1,j ,k ));
+ return alphaz * (F(i ,j+1,k ) - F(i ,j ,k ))
+ + betazx * (F(i+1,j+1,k ) - F(i+1,j ,k )
+ + F(i-1,j+1,k ) - F(i-1,j ,k ));
#endif
};
/**
- /* Perform derivative along z on a nodal grid, from a cell-centered field `F`*/
+ * Perform derivative along z on a nodal grid, from a cell-centered field `F` */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DownwardDz (
amrex::Array4<amrex::Real> const& F,
@@ -210,9 +216,9 @@ struct CartesianCKCAlgorithm {
amrex::Real const inv_dz = coefs_z[0];
#if defined WARPX_DIM_3D
- return inv_dz*( F(i,j,k) - F(i,j,k-1) );
+ return inv_dz*( F(i,j,k) - F(i,j,k-1) );
#elif (defined WARPX_DIM_XZ)
- return inv_dz*( F(i,j,k) - F(i,j-1,k) );
+ return inv_dz*( F(i,j,k) - F(i,j-1,k) );
#endif
};
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-10 20:05:44 +0000