diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
7 files changed, 691 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp new file mode 100644 index 000000000..72d3a1135 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -0,0 +1,184 @@ +#include "WarpXAlgorithmSelection.H" +#include "FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ + #include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else + #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" + #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H" +#endif +#include <AMReX_Gpu.H> + +using namespace amrex; + +void FiniteDifferenceSolver::EvolveB ( VectorField& Bfield, + VectorField const& Efield, + amrex::Real const dt ) { + + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) +#ifdef WARPX_DIM_RZ + if (m_fdtd_algo == MaxwellSolverAlgo::Yee){ + EvolveBCylindrical <CylindricalYeeAlgorithm> ( Bfield, Efield, dt ); +#else + if (m_fdtd_algo == MaxwellSolverAlgo::Yee){ + EvolveBCartesian <YeeAlgorithm> ( Bfield, Efield, dt ); + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + EvolveBCartesian <CKCAlgorithm> ( Bfield, Efield, dt ); +#endif + } else { + amrex::Abort("Unknown algorithm"); + } + +} + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveBCartesian ( VectorField& Bfield, + VectorField const& Efield, + amrex::Real const dt ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + auto const& Bx = Bfield[0]->array(mfi); + auto const& By = Bfield[1]->array(mfi); + auto const& Bz = Bfield[2]->array(mfi); + auto const& Ex = Efield[0]->array(mfi); + auto const& Ey = Efield[1]->array(mfi); + auto const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const* AMREX_RESTRICT coefs_x = stencil_coefs_x.dataPtr(); + int const n_coefs_x = stencil_coefs_x.size(); + Real const* AMREX_RESTRICT coefs_y = stencil_coefs_y.dataPtr(); + int const n_coefs_y = stencil_coefs_y.size(); + Real const* AMREX_RESTRICT coefs_z = stencil_coefs_z.dataPtr(); + int const n_coefs_z = stencil_coefs_z.size(); + + // Extract tileboxes for which to loop + const Box& tbx = mfi.tilebox(Bfield[0]->ixType().ixType()); + const Box& tby = mfi.tilebox(Bfield[1]->ixType().ixType()); + const Box& tbz = mfi.tilebox(Bfield[2]->ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tbx, tby, tbz, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bx(i, j, k) += dt * T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k) + - dt * T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + By(i, j, k) += dt * T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k) + - dt * T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bz(i, j, k) += dt * T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k) + - dt * T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k); + } + + ); + + } + +} + +#else // corresponds to ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield, + VectorField const& Efield, + amrex::Real const dt ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + auto const& Br = Bfield[0]->array(mfi); + auto const& Bt = Bfield[1]->array(mfi); + auto const& Bz = Bfield[2]->array(mfi); + auto const& Er = Efield[0]->array(mfi); + auto const& Et = Efield[1]->array(mfi); + auto const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const* AMREX_RESTRICT coefs_r = stencil_coefs_r.dataPtr(); + int const n_coefs_r = stencil_coefs_r.size(); + Real const* AMREX_RESTRICT coefs_z = stencil_coefs_z.dataPtr(); + int const n_coefs_z = stencil_coefs_z.size(); + + // Extract cylindrical specific parameters + Real const dr = m_dr; + int const nmodes = m_nmodes; + Real const rmin = m_rmin; + + // Extract tileboxes for which to loop + const Box& tbr = mfi.tilebox(Bfield[0]->ixType().ixType()); + const Box& tbt = mfi.tilebox(Bfield[1]->ixType().ixType()); + const Box& tbz = mfi.tilebox(Bfield[2]->ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tbr, tbt, tbz, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const r = rmin + i*dr; // r on nodal point (Br is nodal in r) + Br(i, j, 0, 0) += dt * T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 0); // Mode m=0 + for (int m=1; m<nmodes; m++) { // Higher-order modes + Br(i, j, 0, 2*m-1) += dt*( + T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m-1) + - m * T_Algo::DivideByR(Ez, r, dr, m, i, j, 0, 2*m )); // Real part + Br(i, j, 0, 2*m ) += dt*( + T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m ) + + m * T_Algo::DivideByR(Ez, r, dr, m, i, j, 0, 2*m-1)); // Imaginary part + } + // Ensure that Br remains 0 on axis (except for m=1) + if (r==0) { // On axis + Br(i, j, 0, 0) = 0.; // Mode m=0 + for (int m=2; m<nmodes; m++) { // Higher-order modes (but not m=1) + Br(i, j, 0, 2*m-1) = 0.; + Br(i, j, 0, 2*m ) = 0.; + } + } + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bt(i, j, 0, 0) += dt*( + T_Algo::UpwardDr(Ez, coefs_r, n_coefs_r, i, j, 0, 0) + - T_Algo::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 0)); // Mode m=0 + for (int m=1 ; m<nmodes ; m++) { // Higher-order modes + Bt(i, j, 0, 2*m-1) += dt*( + T_Algo::UpwardDr(Ez, coefs_r, n_coefs_r, i, j, 0, 2*m-1) + - T_Algo::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 2*m-1)); // Real part + Bt(i, j, 0, 2*m ) += dt*( + T_Algo::UpwardDr(Ez, coefs_r, n_coefs_r, i, j, 0, 2*m ) + - T_Algo::UpwardDz(Er, coefs_z, n_coefs_z, i, j, 0, 2*m )); // Imaginary part + } + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const r = rmin + (i + 0.5)*dr; // r on a cell-centered grid (Bz is cell-centered in r) + Bz(i, j, 0, 0) += dt*( - T_Algo::UpwardDrr_over_r(Et, r, dr, coefs_r, n_coefs_r, i, j, 0, 0)); + for (int m=1 ; m<nmodes ; m++) { // Higher-order modes + Bz(i, j, 0, 2*m-1) += dt*( m * Er(i, j, 0, 2*m )/r + - T_Algo::UpwardDrr_over_r(Et, r, dr, coefs_r, n_coefs_r, i, j, 0, 2*m-1)); // Real part + Bz(i, j, 0, 2*m ) += dt*(-m * Er(i, j, 0, 2*m-1)/r + - T_Algo::UpwardDrr_over_r(Et, r, dr, coefs_r, n_coefs_r, i, j, 0, 2*m )); // Imaginary part + } + } + + ); + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H new file mode 100644 index 000000000..5b0d5e718 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H @@ -0,0 +1,210 @@ +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CKC_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_CKC_H_ + +#include <AMReX_REAL.H> +#include <AMReX_Array4.H> +#include <AMReX_Gpu.H> + +struct CKCAlgorithm { + + static void InitializeStencilCoefficients( + std::array<amrex::Real,3>& cell_size, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_x, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_y, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_z ) { + + using namespace amrex; + + // Compute Cole-Karkkainen-Cowan coefficients according + // to Cowan - PRST-AB 16, 041303 (2013) + Real alphax=0, alphay=0, alphaz=0; + Real betaxy=0, betaxz=0, betayx=0, betayz=0, betazx=0, betazy=0; + Real gammax=0, gammay=0, gammaz=0; + Real const inv_dx = 1./cell_size[0]; + Real const inv_dy = 1./cell_size[1]; + Real const inv_dz = 1./cell_size[2]; + #if defined WARPX_DIM_3D + const Real delta = std::max( { inv_dx,inv_dy,inv_dz } ); + const Real rx = (inv_dx/delta)*(inv_dx/delta); + const Real ry = (inv_dy/delta)*(inv_dy/delta); + const Real rz = (inv_dz/delta)*(inv_dz/delta); + const Real beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry)); + betaxy = ry*beta; + betaxz = rz*beta; + betayx = rx*beta; + betayz = rz*beta; + betazx = rx*beta; + betazy = ry*beta; + gammax = ry*rz*(0.0625 - 0.125*ry*rz/(ry*rz + rz*rx + rx*ry)); + gammay = rx*rz*(0.0625 - 0.125*rx*rz/(ry*rz + rz*rx + rx*ry)); + gammaz = rx*ry*(0.0625 - 0.125*rx*ry/(ry*rz + rz*rx + rx*ry)); + alphax = 1. - 2.*betaxy - 2.*betaxz - 4.*gammax; + alphay = 1. - 2.*betayx - 2.*betayz - 4.*gammay; + alphaz = 1. - 2.*betazx - 2.*betazy - 4.*gammaz; + betaxy *= inv_dx; + betaxz *= inv_dx; + betayx *= inv_dy; + betayz *= inv_dy; + betazx *= inv_dz; + betazy *= inv_dz; + alphax *= inv_dx; + alphay *= inv_dy; + alphaz *= inv_dz; + gammax *= inv_dx; + gammay *= inv_dy; + gammaz *= inv_dz; + #elif defined WARPX_DIM_XZ + const Real delta = std::max(inv_dx,inv_dz); + const Real rx = (inv_dx/delta)*(inv_dx/delta); + const Real rz = (inv_dz/delta)*(inv_dz/delta); + betaxz = 0.125*rz; + betazx = 0.125*rx; + alphax = 1. - 2.*betaxz; + alphaz = 1. - 2.*betazx; + betaxz *= inv_dx; + betazx *= inv_dz; + alphax *= inv_dx; + alphaz *= inv_dz; + #endif + + // Store the coefficients in array `stencil_coefs`, in prescribed order + stencil_coefs_x.resize(6); + stencil_coefs_x[0] = inv_dx; + stencil_coefs_x[1] = alphax; + stencil_coefs_x[2] = betaxy; + stencil_coefs_x[3] = betaxz; + stencil_coefs_x[4] = gammax; + stencil_coefs_y.resize(6); + stencil_coefs_y[0] = inv_dy; + stencil_coefs_y[1] = alphay; + stencil_coefs_y[2] = betayz; + stencil_coefs_y[3] = betayx; + stencil_coefs_y[4] = gammay; + stencil_coefs_z.resize(6); + stencil_coefs_z[0] = inv_dz; + stencil_coefs_z[1] = alphaz; + stencil_coefs_z[2] = betazx; + stencil_coefs_z[3] = betazy; + stencil_coefs_z[4] = gammaz; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDx( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real alphax = coefs_x[1]; + amrex::Real betaxy = coefs_x[2]; + amrex::Real betaxz = coefs_x[3]; + amrex::Real 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)); + #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 )); + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDx( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real inv_dx = coefs_x[0]; + return inv_dx*( F(i,j,k) - F(i-1,j,k) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDy( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_y, int const n_coefs_y, + int const i, int const j, int const k ) { + + #if defined WARPX_DIM_3D + amrex::Real alphay = coefs_y[1]; + amrex::Real betayz = coefs_y[2]; + amrex::Real betayx = coefs_y[3]; + amrex::Real 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 + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDy( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_y, int const n_coefs_y, + int const i, int const j, int const k ) { + + #if defined WARPX_DIM_3D + amrex::Real 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 + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDz( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_z, int const n_coefs_z, + int const i, int const j, int const k ) { + + amrex::Real alphaz = coefs_z[1]; + amrex::Real betazx = coefs_z[2]; + amrex::Real betazy = coefs_z[3]; + amrex::Real 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 )); + #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 )); + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDz( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_z, int const n_coefs_z, + int const i, int const j, int const k ) { + + amrex::Real inv_dz = coefs_z[0]; + #if defined WARPX_DIM_3D + 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) ); + #endif + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CKC_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H new file mode 100644 index 000000000..5d160d865 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H @@ -0,0 +1,98 @@ +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_ + +#include <AMReX_REAL.H> +#include <AMReX_Array4.H> +#include <AMReX_Gpu.H> + +struct CylindricalYeeAlgorithm { + + static void InitializeStencilCoefficients( + std::array<amrex::Real,3>& cell_size, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_r, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_z ) { + + // Store the inverse cell size along each direction in the coefficients + stencil_coefs_r.resize(1); + stencil_coefs_r[0] = 1./cell_size[0]; // 1./dr + stencil_coefs_z.resize(1); + stencil_coefs_z[0] = 1./cell_size[2]; // 1./dz + } + + /** Applies the differential operator `1/r * d(rF)/dr`, + * where `F` is on a *nodal* grid in `r` + * and the differential operator is evaluated on *cell-centered* grid. + * The input parameter `r` is given at the cell-centered position*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDrr_over_r( + amrex::Array4<amrex::Real> const& F, + amrex::Real const r, int const dr, + amrex::Real const* coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real inv_dr = coefs_r[0]; + return inv_dr*( (r+0.5*dr)*F(i+1,j,k,comp) - (r-0.5*dr)*F(i,j,k,comp) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDr( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real inv_dr = coefs_r[0]; + return inv_dr*( F(i+1,j,k,comp) - F(i,j,k,comp) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDr( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real inv_dr = coefs_r[0]; + return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDz( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_z, int const n_coefs_z, + int const i, int const j, int const k, int const comp ) { + + amrex::Real inv_dz = coefs_z[0]; + return inv_dz*( F(i,j+1,k,comp) - F(i,j,k,comp) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDz( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_z, int const n_coefs_z, + int const i, int const j, int const k, int const comp ) { + + amrex::Real inv_dz = coefs_z[0]; + return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DivideByR( + amrex::Array4<amrex::Real> const& F, + amrex::Real const r, int const dr, int const m, + int const i, int const j, int const k, int const comp) { + + if (r != 0) { + return F(i,j,k,comp)/r; + } else { // r==0 ; singularity when dividing by r + if (m==1) { + // For m==1, F is linear in r, for small r + // Therefore, the formula below regularizes the singularity + return F(i+1,j,k,comp)/dr; + } else { + return 0; + } + } + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H new file mode 100644 index 000000000..54057091d --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H @@ -0,0 +1,103 @@ +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_ + +#include <AMReX_REAL.H> +#include <AMReX_Array4.H> +#include <AMReX_Gpu.H> + +struct YeeAlgorithm { + + static void InitializeStencilCoefficients( + std::array<amrex::Real,3>& cell_size, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_x, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_y, + amrex::Gpu::ManagedVector<amrex::Real>& stencil_coefs_z ) { + + // Store the inverse cell size along each direction in the coefficients + stencil_coefs_x.resize(1); + stencil_coefs_x[0] = 1./cell_size[0]; + stencil_coefs_y.resize(1); + stencil_coefs_y[0] = 1./cell_size[1]; + stencil_coefs_z.resize(1); + stencil_coefs_z[0] = 1./cell_size[2]; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDx( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real inv_dx = coefs_x[0]; + return inv_dx*( F(i+1,j,k) - F(i,j,k) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDx( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real inv_dx = coefs_x[0]; + return inv_dx*( F(i,j,k) - F(i-1,j,k) ); + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDy( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_y, int const n_coefs_y, + int const i, int const j, int const k ) { + + #if defined WARPX_DIM_3D + amrex::Real inv_dy = coefs_y[0]; + return inv_dy*( F(i,j+1,k) - F(i,j,k) ); + #elif (defined WARPX_DIM_XZ) + return 0; // 2D Cartesian: derivative along y is 0 + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDy( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_y, int const n_coefs_y, + int const i, int const j, int const k ) { + + #if defined WARPX_DIM_3D + amrex::Real 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 + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDz( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_z, int const n_coefs_z, + int const i, int const j, int const k ) { + + amrex::Real inv_dz = coefs_z[0]; + #if defined WARPX_DIM_3D + return inv_dz*( F(i,j,k+1) - F(i,j,k) ); + #elif (defined WARPX_DIM_XZ) + return inv_dz*( F(i,j+1,k) - F(i,j,k) ); + #endif + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDz( + amrex::Array4<amrex::Real> const& F, + amrex::Real const* coefs_z, int const n_coefs_z, + int const i, int const j, int const k ) { + + amrex::Real inv_dz = coefs_z[0]; + #if defined WARPX_DIM_3D + 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) ); + #endif + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H new file mode 100644 index 000000000..9ce910e3d --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -0,0 +1,60 @@ +#ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_ +#define WARPX_FINITE_DIFFERENCE_SOLVER_H_ + +#include "WarpXAlgorithmSelection.H" +#ifdef WARPX_DIM_RZ + #include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else + #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" + #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H" +#endif +#include <AMReX_MultiFab.H> + +/** + * \brief Top-level class for the electromagnetic finite-difference solver + * + * TODO + */ +class FiniteDifferenceSolver +{ + public: + + using VectorField = std::array< std::unique_ptr<amrex::MultiFab>, 3 >; + + // Constructor + FiniteDifferenceSolver ( int const fdtd_algo, + std::array<amrex::Real,3> cell_size ); + + void EvolveB ( VectorField& Bfield, + VectorField const& Efield, + amrex::Real const dt ); + private: + + int m_fdtd_algo; + +#ifdef WARPX_DIM_RZ + amrex::Real m_dr, m_rmin; + amrex::Real m_nmodes; + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_r; + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z; +#else + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_x; + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_y; + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z; +#endif + +#ifdef WARPX_DIM_RZ + template< typename T_Algo > + void EvolveBCylindrical ( VectorField& Bfield, + VectorField const& Efield, + amrex::Real const dt ); +#else + template< typename T_Algo > + void EvolveBCartesian ( VectorField& Bfield, + VectorField const& Efield, + amrex::Real const dt ); +#endif + +}; + +#endif // WARPX_FINITE_DIFFERENCE_SOLVER_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp new file mode 100644 index 000000000..d89e6e9d3 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp @@ -0,0 +1,30 @@ +#include "FiniteDifferenceSolver.H" +#include "WarpX.H" + +// Constructor +FiniteDifferenceSolver::FiniteDifferenceSolver ( int const fdtd_algo, + std::array<amrex::Real,3> cell_size ) { + + // Register the type of finite-difference algorithm + m_fdtd_algo = fdtd_algo; + + // Calculate coefficients of finite-difference stencil +#ifdef WARPX_DIM_RZ + m_dr = cell_size[0]; + m_nmodes = WarpX::GetInstance().n_rz_azimuthal_modes; + m_rmin = WarpX::GetInstance().Geom(0).ProbLo(0); + if (fdtd_algo == MaxwellSolverAlgo::Yee){ + CylindricalYeeAlgorithm::InitializeStencilCoefficients( cell_size, + stencil_coefs_r, stencil_coefs_z ); +#else + if (fdtd_algo == MaxwellSolverAlgo::Yee){ + YeeAlgorithm::InitializeStencilCoefficients( cell_size, + stencil_coefs_x, stencil_coefs_y, stencil_coefs_z ); + } else if (fdtd_algo == MaxwellSolverAlgo::CKC) { + CKCAlgorithm::InitializeStencilCoefficients( cell_size, + stencil_coefs_x, stencil_coefs_y, stencil_coefs_z ); +#endif + } else { + amrex::Abort("Unknown algorithm"); + } +}; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package new file mode 100644 index 000000000..289ed98f2 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += FiniteDifferenceSolver.H +CEXE_sources += FiniteDifferenceSolver.cpp +CEXE_sources += EvolveB.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver |