diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
3 files changed, 205 insertions, 8 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index a78fc8602..a1f35f903 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -1,7 +1,11 @@ #include "WarpXAlgorithmSelection.H" -#include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" -#include "FiniteDifferenceAlgorithms/CKCAlgorithm.H" #include "FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ + #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" +#else + #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" + #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H" +#endif #include <AMReX_Gpu.H> using namespace amrex; @@ -9,19 +13,27 @@ using namespace amrex; void FiniteDifferenceSolver::EvolveB ( VectorField& Bfield, VectorField const& Efield, amrex::Real const dt ) { + +#ifdef WARPX_DIM_RZ + EvolveBCylindrical( Bfield, Efield, dt ); +#else // Select algorithm (The choice of algorithm is a runtime option, // but we compile code for each algorithm, using templates) if (m_fdtd_algo == MaxwellSolverAlgo::Yee){ - EvolveBwithAlgo <YeeAlgorithm> ( Bfield, Efield, dt ); + EvolveBCartesian <YeeAlgorithm> ( Bfield, Efield, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { - EvolveBwithAlgo <CKCAlgorithm> ( Bfield, Efield, dt ); + EvolveBCartesian <CKCAlgorithm> ( Bfield, Efield, dt ); } else { amrex::Abort("Unknown algorithm"); } +#endif + } +#ifndef WARPX_DIM_RZ + template<typename T_Algo> -void FiniteDifferenceSolver::EvolveBwithAlgo ( VectorField& Bfield, +void FiniteDifferenceSolver::EvolveBCartesian ( VectorField& Bfield, VectorField const& Efield, amrex::Real const dt ) { @@ -75,3 +87,76 @@ void FiniteDifferenceSolver::EvolveBwithAlgo ( VectorField& Bfield, } } + +#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_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& 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){ + Br(i, j, 0, 0) += dt * T_Algo::UpwardDz(Et, inv_dz, 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, inv_dz, i, j, 0, 2*m-1) + - m * T_Algo::DivideByR(Ez, r, dr, i, j, 0, 2*m )); // Real part + Br(i, j, 0, 2*m ) += dt*( T_Algo::UpwardDz(Et, inv_dz, i, j, 0, 2*m ) + + m * T_Algo::DivideByR(Ez, r, dr, 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::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); + } + + ); + + } + +} + +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H new file mode 100644 index 000000000..3f408479c --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H @@ -0,0 +1,103 @@ +#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 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_CYLINDRICAL_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 4c0e1955d..c08810f85 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -2,8 +2,10 @@ #define WARPX_FINITE_DIFFERENCE_SOLVER_H_ #include "WarpXAlgorithmSelection.H" -#include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" -#include "FiniteDifferenceAlgorithms/CKCAlgorithm.H" +#ifndef WARPX_DIM_RZ + #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H" + #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H" +#endif #include <AMReX_MultiFab.H> /** @@ -47,10 +49,17 @@ class FiniteDifferenceSolver amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_y; amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z; +#ifdef WARPX_DIM_RZ template< typename T_Algo > - void EvolveBwithAlgo ( VectorField& Bfield, + 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 }; |