diff options
Diffstat (limited to 'Source/FieldSolver')
5 files changed, 208 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp new file mode 100644 index 000000000..e5514c697 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -0,0 +1,68 @@ +// TODO include statements + +FiniteDifferenceSolver::EvolveB ( VectorField Bfield, + ConstVectorField Efield, + amrex::Real dt ) { + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) + if (fdtd_algo == MaxwellSolverAlgo::Yee){ + EvolveB <YeeAlgorithm> ( Bfield, Efield, dt ); + } else if (fdtd_algo == MaxwellSolverAlgo::CKC) { + EvolveB <CKCAlgorithm> ( Bfield, Efield, dt ); + } else { + amrex::Abort("Unknown algorithm"); + } +) + +template<typename algo> +FiniteDifferenceSolver::EvolveB ( VectorField Bfield, + ConstVectorField Efield, + amrex::Real dt ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bx, 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(); + Real const* AMREX_RESTRICT coefs_y = stencil_coefs_y.dataPtr(); + Real const* AMREX_RESTRICT coefs_z = stencil_coefs_z.dataPtr(); + + // Extract tileboxes for which to loop + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + // 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 * algo::UpwardDz( Ey, i, j, k, coefs_z) + - dt * algo::UpwardDy( Ez, i, j, k, coefs_y); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + By(i, j, k) += dt * algo::UpwardDx( Ez, i, j, k, coefs_x) + - dt * algo::UpwardDz( Ex, i, j, k, coefs_z); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bz(i, j, k) += dt * algo::UpwardDy( Ex, i, j, k, coefs_y) + - dt * algo::UpwardDx( Ey, i, j, k, coefs_x); + } + + ); + + } + +}; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H new file mode 100644 index 000000000..7c9572298 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H @@ -0,0 +1,80 @@ +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_YEE_H_ + +struct YeeAlgorithm { + + void InitializeStencilCoefficients( std::array<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 + amrex::Real UpwardDx( + Array4 const F, int const i, int const j, int const k, + Real const* coefs_x, int const n_coefs_x ) { + + 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 + amrex::Real DownwardDx( + Array4 const F, int const i, int const j, int const k, + Real const* coefs_x, int const n_coefs_x ) { + + 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 + amrex::Real UpwardDy( + Array4 const F, int const i, int const j, int const k, + Real const* coefs_y, int const n_coefs_y ) { + + #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 + amrex::Real DownwardDy( + Array4 const F, int const i, int const j, int const k, + Real const* coefs_y, int const n_coefs_y ) { + + #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 + amrex::Real UpwardDz( + Array4 const F, int const i, int const j, int const k, + Real const* coefs_z, int const n_coefs_z ) { + + 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 + }; + + +} + +#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..c1fdef5dd --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -0,0 +1,38 @@ +#ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_ +#define WARPX_FINITE_DIFFERENCE_SOLVER_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 >; + using ConstVectorField = std::array< std::unique_ptr<amrex::MultiFab const>, 3 >; + + // Constructor + void FiniteDifferenceSolver( std::array<Real,3> dx ); + + void EvolveB( VectorField Bfield, + ConstVectorField Efield, + amrex::Real dt ) const; + + }; + + private: + + int fdtd_algo; + + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_x; + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_y; + amrex::Gpu::ManagedVector<amrex::Real> stencil_coefs_z; + + template< typename fdtd_algo > + void EvolveB + +}; + +#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..a146b66b1 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp @@ -0,0 +1,15 @@ + + +void FiniteDifferenceSolver::FiniteDifferenceSolver ( std::array<Real,3> cell_size ) { + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) + 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 ); + } else { + amrex::Abort("Unknown algorithm"); + } +}; diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 4848b051e..1137da741 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -111,6 +111,13 @@ WarpX::EvolveB (int lev, amrex::Real a_dt) void WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { + if (patch_type == PatchType::fine) { + fdtd_solver_fp->EvolveB( Bfield_fp[lev], Efield_fp[lev], a_dt ); + } else { + fdtd_solver_cp->EvolveB( Bfield_cp[lev], Efield_cp[lev], a_dt ); + } + + // Goes into initializer const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array<Real,3>& dx = WarpX::CellSize(patch_level); const Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; |