diff options
Diffstat (limited to 'Source/FieldSolver')
12 files changed, 959 insertions, 377 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp new file mode 100644 index 000000000..16227febc --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -0,0 +1,221 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "WarpXAlgorithmSelection.H" +#include "FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ +# include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#include <AMReX_Gpu.H> + +using namespace amrex; + +/** + * \brief Update the B field, over one timestep + */ +void FiniteDifferenceSolver::EvolveB ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > 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_do_nodal) { + + EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, dt ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, dt ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + + EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, dt ); + +#endif + } else { + amrex::Abort("Unknown algorithm"); + } + +} + + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveBCartesian ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > 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 + Array4<Real> const& Bx = Bfield[0]->array(mfi); + Array4<Real> const& By = Bfield[1]->array(mfi); + Array4<Real> const& Bz = Bfield[2]->array(mfi); + Array4<Real> const& Ex = Efield[0]->array(mfi); + Array4<Real> const& Ey = Efield[1]->array(mfi); + Array4<Real> const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = m_stencil_coefs_x.size(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // Extract tileboxes for which to loop + Box const& tbx = mfi.tilebox(Bfield[0]->ixType().ixType()); + Box const& tby = mfi.tilebox(Bfield[1]->ixType().ixType()); + Box const& 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 ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > 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 + Array4<Real> const& Br = Bfield[0]->array(mfi); + Array4<Real> const& Bt = Bfield[1]->array(mfi); + Array4<Real> const& Bz = Bfield[2]->array(mfi); + Array4<Real> const& Er = Efield[0]->array(mfi); + Array4<Real> const& Et = Efield[1]->array(mfi); + Array4<Real> const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_r = m_stencil_coefs_r.dataPtr(); + int const n_coefs_r = m_stencil_coefs_r.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_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 + Box const& tbr = mfi.tilebox(Bfield[0]->ixType().ixType()); + Box const& tbt = mfi.tilebox(Bfield[1]->ixType().ixType()); + Box const& 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) + if (r != 0) { // Off-axis, regular Maxwell equations + 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 * Ez(i, j, 0, 2*m )/r ); // Real part + Br(i, j, 0, 2*m ) += dt*( + T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m ) + + m * Ez(i, j, 0, 2*m-1)/r ); // Imaginary part + } + } else { // r==0: On-axis corrections + // Ensure that Br remains 0 on axis (except for m=1) + Br(i, j, 0, 0) = 0.; // Mode m=0 + for (int m=1; m<nmodes; m++) { // Higher-order modes + if (m == 1){ + // For m==1, Ez is linear in r, for small r + // Therefore, the formula below regularizes the singularity + Br(i, j, 0, 2*m-1) += dt*( + T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m-1) + - m * Ez(i+1, j, 0, 2*m )/dr ); // Real part + Br(i, j, 0, 2*m ) += dt*( + T_Algo::UpwardDz(Et, coefs_z, n_coefs_z, i, j, 0, 2*m ) + + m * Ez(i+1, j, 0, 2*m-1)/dr ); // Imaginary part + } else { + 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/CartesianCKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H new file mode 100644 index 000000000..fa5dd073d --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H @@ -0,0 +1,221 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_ + +#include <AMReX_REAL.H> +#include <AMReX_Array4.H> +#include <AMReX_Gpu.H> + +/** + * This struct contains only static functions to initialize the stencil coefficients + * and to compute finite-difference derivatives for the Cartesian CKC algorithm. + */ +struct CartesianCKCAlgorithm { + + 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 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 + 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; +#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; +#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; + } + + /** + /* 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, + amrex::Real const * const coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real const alphax = coefs_x[1]; + amrex::Real const betaxy = coefs_x[2]; + 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)); +#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 + }; + + /** + /* 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, + amrex::Real const * const coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real const inv_dx = coefs_x[0]; + return inv_dx*( F(i,j,k) - F(i-1,j,k) ); + }; + + /** + /* 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 ) { + +#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)); +#elif (defined WARPX_DIM_XZ) + return 0; // 2D Cartesian: derivative along y is 0 +#endif + }; + + /** + /* 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 ) { + +#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) ); +#elif (defined WARPX_DIM_XZ) + return 0; // 2D Cartesian: derivative along y is 0 +#endif + }; + + /** + /* 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]; +#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 + }; + + /** + /* 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, + amrex::Real const * const coefs_z, int const n_coefs_z, + int const i, int const j, int const k ) { + + 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) ); +#elif (defined WARPX_DIM_XZ) + return inv_dz*( F(i,j,k) - F(i,j-1,k) ); +#endif + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H new file mode 100644 index 000000000..69622c5fe --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H @@ -0,0 +1,130 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_NODAL_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_NODAL_H_ + +#include <AMReX_REAL.H> +#include <AMReX_Array4.H> +#include <AMReX_Gpu.H> + +/** + * This struct contains only static functions to initialize the stencil coefficients + * and to compute finite-difference derivatives for the Cartesian nodal algorithm. + */ +struct CartesianNodalAlgorithm { + + 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]; + } + + /** + /* Perform derivative along x + /* (For a solver on a staggered grid, `UpwardDx` and `DownwardDx` take into + /* account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDx ( + amrex::Array4<amrex::Real> const& F, + amrex::Real const * const coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real const inv_dx = coefs_x[0]; + return 0.5*inv_dx*( F(i+1,j,k) - F(i-1,j,k) ); + }; + + /** + /* Perform derivative along x + /* (For a solver on a staggered grid, `UpwardDx` and `DownwardDx` take into + /* account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDx ( + amrex::Array4<amrex::Real> const& F, + amrex::Real const * const coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + return UpwardDx( F, coefs_x, n_coefs_x, i, j, k ); + // For CartesianNodalAlgorithm, UpwardDx and DownwardDx are equivalent + }; + + /** + /* Perform derivative along y + /* (For a solver on a staggered grid, `UpwardDy` and `DownwardDy` take into + /* account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ + 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 ) { + +#if defined WARPX_DIM_3D + amrex::Real const inv_dy = coefs_y[0]; + return 0.5*inv_dy*( F(i,j+1,k) - F(i,j-1,k) ); +#elif (defined WARPX_DIM_XZ) + return 0; // 2D Cartesian: derivative along y is 0 +#endif + }; + + /** + /* Perform derivative along y + /* (For a solver on a staggered grid, `UpwardDy` and `DownwardDy` take into + /* account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ + 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 ) { + + return UpwardDy( F, coefs_y, n_coefs_y, i, j, k ); + // For CartesianNodalAlgorithm, UpwardDy and DownwardDy are equivalent + }; + + /** + /* Perform derivative along z + /* (For a solver on a staggered grid, `UpwardDz` and `DownwardDz` take into + /* account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ + 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 inv_dz = coefs_z[0]; +#if defined WARPX_DIM_3D + return 0.5*inv_dz*( F(i,j,k+1) - F(i,j,k-1) ); +#elif (defined WARPX_DIM_XZ) + return 0.5*inv_dz*( F(i,j+1,k) - F(i,j-1,k) ); +#endif + }; + + /** + /* Perform derivative along z + /* (For a solver on a staggered grid, `UpwardDz` and `DownwardDz` take into + /* account the staggering; but for `CartesianNodalAlgorithm`, they are equivalent) */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDz ( + 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 ) { + + return UpwardDz( F, coefs_z, n_coefs_z, i, j, k ); + // For CartesianNodalAlgorithm, UpwardDz and DownwardDz are equivalent + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_NODAL_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H new file mode 100644 index 000000000..268c5aa89 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -0,0 +1,126 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ +#define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ + +#include <AMReX_REAL.H> +#include <AMReX_Array4.H> +#include <AMReX_Gpu.H> + +/** + * This struct contains only static functions to initialize the stencil coefficients + * and to compute finite-difference derivatives for the Cartesian Yee algorithm. + */ +struct CartesianYeeAlgorithm { + + 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]; + } + + /** + /* 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, + amrex::Real const * const coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real const inv_dx = coefs_x[0]; + return inv_dx*( F(i+1,j,k) - F(i,j,k) ); + }; + + /** + /* 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, + amrex::Real const * const coefs_x, int const n_coefs_x, + int const i, int const j, int const k ) { + + amrex::Real const inv_dx = coefs_x[0]; + return inv_dx*( F(i,j,k) - F(i-1,j,k) ); + }; + + /** + /* 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 ) { + +#if defined WARPX_DIM_3D + amrex::Real const 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 + }; + + /** + /* 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 ) { + +#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) ); +#elif (defined WARPX_DIM_XZ) + return 0; // 2D Cartesian: derivative along y is 0 +#endif + }; + + /** + /* 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 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 + }; + + /** + /* 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, + amrex::Real const * const coefs_z, int const n_coefs_z, + int const i, int const j, int const k ) { + + 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) ); +#elif (defined WARPX_DIM_XZ) + return inv_dz*( F(i,j,k) - F(i,j-1,k) ); +#endif + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H new file mode 100644 index 000000000..ab32c8bcb --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H @@ -0,0 +1,113 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#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> + +/** + * This struct contains only static functions to initialize the stencil coefficients + * and to compute finite-difference derivatives for the Cartesian Yee algorithm. + */ +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 a *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, amrex::Real const dr, + amrex::Real const * const coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real const inv_dr = coefs_r[0]; + return 1./r * inv_dr*( (r+0.5*dr)*F(i+1,j,k,comp) - (r-0.5*dr)*F(i,j,k,comp) ); + }; + + /** Applies the differential operator `1/r * d(rF)/dr`, + * where `F` is on a *cell-centered* grid in `r` + * and the differential operator is evaluated on a *nodal* grid. + * The input parameter `r` is given at the cell-centered position */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDrr_over_r ( + amrex::Array4<amrex::Real> const& F, + amrex::Real const r, amrex::Real const dr, + amrex::Real const * const coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real const inv_dr = coefs_r[0]; + return 1./r * inv_dr*( (r+0.5*dr)*F(i,j,k,comp) - (r-0.5*dr)*F(i-1,j,k,comp) ); + }; + + /** + /* Perform derivative along r on a cell-centered grid, from a nodal field `F`*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real UpwardDr ( + amrex::Array4<amrex::Real> const& F, + amrex::Real const * const coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real const inv_dr = coefs_r[0]; + return inv_dr*( F(i+1,j,k,comp) - F(i,j,k,comp) ); + }; + + /** + /* Perform derivative along r on a nodal grid, from a cell-centered field `F`*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real DownwardDr ( + amrex::Array4<amrex::Real> const& F, + amrex::Real const * const coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { + + amrex::Real const inv_dr = coefs_r[0]; + return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) ); + }; + + /** + /* 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, int const comp ) { + + amrex::Real const inv_dz = coefs_z[0]; + return inv_dz*( F(i,j+1,k,comp) - F(i,j,k,comp) ); + }; + + /** + /* 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, + amrex::Real const * const coefs_z, int const n_coefs_z, + int const i, int const j, int const k, int const comp ) { + + amrex::Real const inv_dz = coefs_z[0]; + return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) ); + }; + +}; + +#endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H new file mode 100644 index 000000000..4bf88077f --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -0,0 +1,76 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_ +#define WARPX_FINITE_DIFFERENCE_SOLVER_H_ + +#include <AMReX_MultiFab.H> + +/** + * \brief Top-level class for the electromagnetic finite-difference solver + * + * Stores the coefficients of the finite-difference stencils, + * and has member functions to update fields over one time step. + */ +class FiniteDifferenceSolver +{ + public: + + // Constructor + /** \brief Initialize the finite-difference Maxwell solver (for a given refinement level) + * + * This function initializes the stencil coefficients for the chosen finite-difference algorithm + * + * \param fdtd_algo Identifies the chosen algorithm, as defined in WarpXAlgorithmSelection.H + * \param cell_size Cell size along each dimension, for the chosen refinement level + * \param do_nodal Whether the solver is applied to a nodal or staggered grid + */ + FiniteDifferenceSolver ( + int const fdtd_algo, + std::array<amrex::Real,3> cell_size, + bool const do_nodal ); + + void EvolveB ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, + amrex::Real const dt ); + private: + + int m_fdtd_algo; + bool m_do_nodal; + +#ifdef WARPX_DIM_RZ + amrex::Real m_dr, m_rmin; + amrex::Real m_nmodes; + amrex::Gpu::ManagedVector<amrex::Real> m_stencil_coefs_r; + amrex::Gpu::ManagedVector<amrex::Real> m_stencil_coefs_z; +#else + amrex::Gpu::ManagedVector<amrex::Real> m_stencil_coefs_x; + amrex::Gpu::ManagedVector<amrex::Real> m_stencil_coefs_y; + amrex::Gpu::ManagedVector<amrex::Real> m_stencil_coefs_z; +#endif + + public: + // The member functions below contain extended __device__ lambda. + // In order to compile with nvcc, they need to be public. + +#ifdef WARPX_DIM_RZ + template< typename T_Algo > + void EvolveBCylindrical ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, + amrex::Real const dt ); +#else + template< typename T_Algo > + void EvolveBCartesian ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > 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..ea7af6677 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp @@ -0,0 +1,58 @@ +/* Copyright 2020 Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "WarpXAlgorithmSelection.H" +#ifdef WARPX_DIM_RZ +# include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" +#else +# include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +# include "FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H" +#endif +#include "FiniteDifferenceSolver.H" +#include "WarpX.H" + +/* This function initializes the stencil coefficients for the chosen finite-difference algorithm */ +FiniteDifferenceSolver::FiniteDifferenceSolver ( + int const fdtd_algo, + std::array<amrex::Real,3> cell_size, + bool do_nodal ) { + + // Register the type of finite-difference algorithm + m_fdtd_algo = fdtd_algo; + m_do_nodal = do_nodal; + + // 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, + m_stencil_coefs_r, m_stencil_coefs_z ); +#else + if (do_nodal) { + + CartesianNodalAlgorithm::InitializeStencilCoefficients( cell_size, + m_stencil_coefs_x, m_stencil_coefs_y, m_stencil_coefs_z ); + + } else if (fdtd_algo == MaxwellSolverAlgo::Yee) { + + CartesianYeeAlgorithm::InitializeStencilCoefficients( cell_size, + m_stencil_coefs_x, m_stencil_coefs_y, m_stencil_coefs_z ); + + } else if (fdtd_algo == MaxwellSolverAlgo::CKC) { + + CartesianCKCAlgorithm::InitializeStencilCoefficients( cell_size, + m_stencil_coefs_x, m_stencil_coefs_y, m_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 diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package index 018cfbfba..522c4c07a 100644 --- a/Source/FieldSolver/Make.package +++ b/Source/FieldSolver/Make.package @@ -7,6 +7,7 @@ ifeq ($(USE_PSATD),TRUE) include $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package endif endif +include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 085f34441..3c4b4439c 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -118,129 +118,16 @@ WarpX::EvolveB (int lev, amrex::Real a_dt) void WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { - 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]; - const Real dxinv = 1./dx[0]; - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; - if (patch_type == PatchType::fine) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], a_dt ); + } else { + m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], a_dt ); } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); - } - - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - - // xmin is only used by the kernel for cylindrical geometry, - // in which case it is actually rmin. - const Real xmin = Geom(0).ProbLo(0); - - // 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 ) - { - Real wt = amrex::second(); - - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - auto const& Bxfab = Bx->array(mfi); - auto const& Byfab = By->array(mfi); - auto const& Bzfab = Bz->array(mfi); - auto const& Exfab = Ex->array(mfi); - auto const& Eyfab = Ey->array(mfi); - auto const& Ezfab = Ez->array(mfi); - if (do_nodal) { - amrex::ParallelFor(tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_bx_nodal(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_by_nodal(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); - }); - } else if (WarpX::maxwell_fdtd_solver_id == 0) { - const long nmodes = n_rz_azimuthal_modes; - amrex::ParallelFor(tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdx,dtsdy,dtsdz,dxinv,xmin,nmodes); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz,nmodes); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin,nmodes); - }); - } else if (WarpX::maxwell_fdtd_solver_id == 1) { - Real betaxy, betaxz, betayx, betayz, betazx, betazy; - Real gammax, gammay, gammaz; - Real alphax, alphay, alphaz; - warpx_calculate_ckc_coefficients(dtsdx, dtsdy, dtsdz, - betaxy, betaxz, betayx, betayz, betazx, betazy, - gammax, gammay, gammaz, - alphax, alphay, alphaz); - amrex::ParallelFor(tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_bx_ckc(j,k,l,Bxfab,Eyfab,Ezfab, - betaxy, betaxz, betayx, betayz, betazx, betazy, - gammax, gammay, gammaz, - alphax, alphay, alphaz); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_by_ckc(j,k,l,Byfab,Exfab,Ezfab, - betaxy, betaxz, betayx, betayz, betazx, betazy, - gammax, gammay, gammaz, - alphax, alphay, alphaz); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) - { - warpx_push_bz_ckc(j,k,l,Bzfab,Exfab,Eyfab, - betaxy, betaxz, betayx, betayz, betazx, betazy, - gammax, gammay, gammaz, - alphax, alphay, alphaz); - }); - } - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (patch_type == PatchType::coarse) cbx.refine(rr); - wt = (amrex::second() - wt) / cbx.d_numPts(); - auto costfab = cost->array(mfi); - amrex::ParallelFor(cbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - costfab(i,j,k) += wt; - }); - } - } + 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]; if (do_pml && pml[lev]->ok()) { diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 4ad251264..8a2202059 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -10,105 +10,6 @@ #include <AMReX_FArrayBox.H> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_bx_yee(int i, int j, int k, - amrex::Array4<amrex::Real> const& Bx, - amrex::Array4<amrex::Real const> const& Ey, - amrex::Array4<amrex::Real const> const& Ez, - amrex::Real dtsdx, amrex::Real dtsdy, amrex::Real dtsdz, - amrex::Real dxinv, amrex::Real rmin, const long nmodes) -{ -#if defined WARPX_DIM_3D - Bx(i,j,k) += - dtsdy * (Ez(i,j+1,k ) - Ez(i,j,k)) - + dtsdz * (Ey(i,j ,k+1) - Ey(i,j,k)); -#elif (defined WARPX_DIM_XZ) - Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0)); -#elif (defined WARPX_DIM_RZ) - if (i != 0 || rmin != 0.) { - Bx(i,j,0,0) += + dtsdz * (Ey(i,j+1,0,0) - Ey(i,j,0,0)); - } else { - Bx(i,j,0,0) = 0.; - } - for (int imode=1 ; imode < nmodes ; imode++) { - if (i == 0 && rmin == 0) { - if (imode == 1) { - // For the mode m = 1, the bulk equation diverges on axis - // (due to the 1/r terms). The following expressions regularize - // these divergences by assuming, on axis : - // Ez/r = 0/r + dEz/dr - Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i+1,j,0,2*imode) - + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); - Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i+1,j,0,2*imode-1) - + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); - } else { - Bx(i,j,0,2*imode-1) = 0.; - Bx(i,j,0,2*imode) = 0.; - } - } else { - // Br(i,j,m) = Br(i,j,m) + I*m*dt*Ez(i,j,m)/r + dtsdz*(Et(i,j+1,m) - Et(i,j,m)) - const amrex::Real r = rmin*dxinv + i; - Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i,j,0,2*imode)/r - + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); - Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i,j,0,2*imode-1)/r - + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); - } - } -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_by_yee(int i, int j, int k, - amrex::Array4<amrex::Real> const& By, - amrex::Array4<amrex::Real const> const& Ex, - amrex::Array4<amrex::Real const> const& Ez, - amrex::Real dtsdx, amrex::Real dtsdz, - const long nmodes) -{ -#if defined WARPX_DIM_3D - By(i,j,k) += + dtsdx * (Ez(i+1,j,k ) - Ez(i,j,k)) - - dtsdz * (Ex(i ,j,k+1) - Ex(i,j,k)); -#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - // Note that the 2D Cartesian and RZ mode 0 are the same - By(i,j,0,0) += + dtsdx * (Ez(i+1,j ,0,0) - Ez(i,j,0,0)) - - dtsdz * (Ex(i ,j+1,0,0) - Ex(i,j,0,0)); -#if (defined WARPX_DIM_RZ) - for (int imode=1 ; imode < nmodes ; imode++) { - // Bt(i,j,m) = Bt(i,j,m) + dtsdr*(Ez(i+1,j,m) - Ez(i,j,m)) - dtsdz*(Er(i,j+1,m) - Er(i,j,m)) - By(i,j,0,2*imode-1) += + dtsdx*(Ez(i+1,j ,0,2*imode-1) - Ez(i,j,0,2*imode-1)) - - dtsdz*(Ex(i ,j+1,0,2*imode-1) - Ex(i,j,0,2*imode-1)); - By(i,j,0,2*imode) += + dtsdx*(Ez(i+1,j ,0,2*imode) - Ez(i,j,0,2*imode)) - - dtsdz*(Ex(i ,j+1,0,2*imode) - Ex(i,j,0,2*imode)); - } -#endif -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_bz_yee(int i, int j, int k, - amrex::Array4<amrex::Real> const& Bz, - amrex::Array4<amrex::Real const> const& Ex, - amrex::Array4<amrex::Real const> const& Ey, - amrex::Real dtsdx, amrex::Real dtsdy, - amrex::Real dxinv, amrex::Real rmin, const long nmodes) -{ -#if defined WARPX_DIM_3D - Bz(i,j,k) += - dtsdx * (Ey(i+1,j ,k) - Ey(i,j,k)) - + dtsdy * (Ex(i ,j+1,k) - Ex(i,j,k)); -#elif defined WARPX_DIM_XZ - Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0)); -#elif defined WARPX_DIM_RZ - const amrex::Real r = rmin*dxinv + i + 0.5; - const amrex::Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); - const amrex::Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); - Bz(i,j,0,0) += - dtsdx*(ru*Ey(i+1,j,0,0) - rd*Ey(i,j,0,0)); - for (int imode=1 ; imode < nmodes ; imode++) { - // Bz(i,j,m) = Bz(i,j,m) - dtsdr*(ru*Et(i+1,j,m) - rd*Et(i,j,m)) - I*m*dt*Er(i,j,m)/r - Bz(i,j,0,2*imode-1) += - dtsdx*(ru*Ey(i+1,j,0,2*imode-1) - rd*Ey(i,j,0,2*imode-1)) + imode*dtsdx*Ex(i,j,0,2*imode)/r; - Bz(i,j,0,2*imode) += - dtsdx*(ru*Ey(i+1,j,0,2*imode) - rd*Ey(i,j,0,2*imode)) - imode*dtsdx*Ex(i,j,0,2*imode-1)/r; - } -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ex_yee(int i, int j, int k, amrex::Array4<amrex::Real> const& Ex, amrex::Array4<amrex::Real const> const& By, @@ -342,117 +243,6 @@ static void warpx_calculate_ckc_coefficients(amrex::Real dtsdx, amrex::Real dtsd } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_bx_ckc(int j, int k, int l, - amrex::Array4<amrex::Real> const& Bx, - amrex::Array4<amrex::Real const> const& Ey, - amrex::Array4<amrex::Real const> const& Ez, - amrex::Real betaxy, amrex::Real betaxz, amrex::Real betayx, - amrex::Real betayz, amrex::Real betazx, amrex::Real betazy, - amrex::Real gammax, amrex::Real gammay, amrex::Real gammaz, - amrex::Real alphax, amrex::Real alphay, amrex::Real alphaz) -{ -#if defined WARPX_DIM_3D - Bx(j,k,l) += - alphay * (Ez(j ,k+1,l ) - Ez(j, k ,l )) - - betayx * (Ez(j+1,k+1,l ) - Ez(j+1,k ,l ) - + Ez(j-1,k+1,l ) - Ez(j-1,k ,l )) - - betayz * (Ez(j ,k+1,l+1) - Ez(j ,k ,l+1) - + Ez(j ,k+1,l-1) - Ez(j ,k ,l-1)) - - gammay * (Ez(j+1,k+1,l+1) - Ez(j+1,k ,l+1) - + Ez(j-1,k+1,l+1) - Ez(j-1,k ,l+1) - + Ez(j+1,k+1,l-1) - Ez(j+1,k ,l-1) - + Ez(j-1,k+1,l-1) - Ez(j-1,k ,l-1)) - + alphaz * (Ey(j ,k ,l+1) - Ey(j, k, l )) - + betazx * (Ey(j+1,k ,l+1) - Ey(j+1,k ,l ) - + Ey(j-1,k ,l+1) - Ey(j-1,k ,l )) - + betazy * (Ey(j ,k+1,l+1) - Ey(j ,k+1,l ) - + Ey(j ,k-1,l+1) - Ey(j ,k-1,l )) - + gammaz * (Ey(j+1,k+1,l+1) - Ey(j+1,k+1,l ) - + Ey(j-1,k+1,l+1) - Ey(j-1,k+1,l ) - + Ey(j+1,k-1,l+1) - Ey(j+1,k-1,l ) - + Ey(j-1,k-1,l+1) - Ey(j-1,k-1,l )); -#elif defined WARPX_DIM_XZ - Bx(j,k,0) += + alphaz * (Ey(j ,k+1,0) - Ey(j, k,0)) - + betazx * (Ey(j+1,k+1,0) - Ey(j+1,k,0) - + Ey(j-1,k+1,0) - Ey(j-1,k,0)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_by_ckc(int j, int k, int l, - amrex::Array4<amrex::Real> const& By, - amrex::Array4<amrex::Real const> const& Ex, - amrex::Array4<amrex::Real const> const& Ez, - amrex::Real betaxy, amrex::Real betaxz, amrex::Real betayx, - amrex::Real betayz, amrex::Real betazx, amrex::Real betazy, - amrex::Real gammax, amrex::Real gammay, amrex::Real gammaz, - amrex::Real alphax, amrex::Real alphay, amrex::Real alphaz) -{ -#if defined WARPX_DIM_3D - By(j,k,l) += + alphax * (Ez(j+1,k ,l ) - Ez(j, k, l )) - + betaxy * (Ez(j+1,k+1,l ) - Ez(j ,k+1,l ) - + Ez(j+1,k-1,l ) - Ez(j ,k-1,l )) - + betaxz * (Ez(j+1,k ,l+1) - Ez(j ,k ,l+1) - + Ez(j+1,k ,l-1) - Ez(j ,k ,l-1)) - + gammax * (Ez(j+1,k+1,l+1) - Ez(j ,k+1,l+1) - + Ez(j+1,k-1,l+1) - Ez(j ,k-1,l+1) - + Ez(j+1,k+1,l-1) - Ez(j ,k+1,l-1) - + Ez(j+1,k-1,l-1) - Ez(j ,k-1,l-1)) - - alphaz * (Ex(j ,k ,l+1) - Ex(j ,k ,l )) - - betazx * (Ex(j+1,k ,l+1) - Ex(j+1,k ,l ) - + Ex(j-1,k ,l+1) - Ex(j-1,k ,l )) - - betazy * (Ex(j ,k+1,l+1) - Ex(j ,k+1,l ) - + Ex(j ,k-1,l+1) - Ex(j ,k-1,l )) - - gammaz * (Ex(j+1,k+1,l+1) - Ex(j+1,k+1,l ) - + Ex(j-1,k+1,l+1) - Ex(j-1,k+1,l ) - + Ex(j+1,k-1,l+1) - Ex(j+1,k-1,l ) - + Ex(j-1,k-1,l+1) - Ex(j-1,k-1,l )); -#elif defined WARPX_DIM_XZ - By(j,k,0) += + alphax * (Ez(j+1,k ,0) - Ez(j,k ,0)) - + betaxz * (Ez(j+1,k+1,0) - Ez(j,k+1,0) - + Ez(j+1,k-1,0) - Ez(j,k-1,0)) - - alphaz * (Ex(j ,k+1,0) - Ex(j,k ,0)) - - betazx * (Ex(j+1,k+1,0) - Ex(j+1,k,0) - + Ex(j-1,k+1,0) - Ex(j-1,k,0)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void warpx_push_bz_ckc(int j, int k, int l, - amrex::Array4<amrex::Real> const& Bz, - amrex::Array4<amrex::Real const> const& Ex, - amrex::Array4<amrex::Real const> const& Ey, - amrex::Real betaxy, amrex::Real betaxz, amrex::Real betayx, - amrex::Real betayz, amrex::Real betazx, amrex::Real betazy, - amrex::Real gammax, amrex::Real gammay, amrex::Real gammaz, - amrex::Real alphax, amrex::Real alphay, amrex::Real alphaz) -{ -#if defined WARPX_DIM_3D - Bz(j,k,l) += - alphax * (Ey(j+1,k ,l ) - Ey(j ,k ,l )) - - betaxy * (Ey(j+1,k+1,l ) - Ey(j ,k+1,l ) - + Ey(j+1,k-1,l ) - Ey(j ,k-1,l )) - - betaxz * (Ey(j+1,k ,l+1) - Ey(j ,k ,l+1) - + Ey(j+1,k ,l-1) - Ey(j ,k ,l-1)) - - gammax * (Ey(j+1,k+1,l+1) - Ey(j ,k+1,l+1) - + Ey(j+1,k-1,l+1) - Ey(j ,k-1,l+1) - + Ey(j+1,k+1,l-1) - Ey(j ,k+1,l-1) - + Ey(j+1,k-1,l-1) - Ey(j ,k-1,l-1)) - + alphay * (Ex(j ,k+1,l ) - Ex(j ,k ,l )) - + betayx * (Ex(j+1,k+1,l ) - Ex(j+1,k ,l ) - + Ex(j-1,k+1,l ) - Ex(j-1,k ,l )) - + betayz * (Ex(j ,k+1,l+1) - Ex(j ,k ,l+1) - + Ex(j ,k+1,l-1) - Ex(j ,k ,l-1)) - + gammay * (Ex(j+1,k+1,l+1) - Ex(j+1,k ,l+1) - + Ex(j-1,k+1,l+1) - Ex(j-1,k ,l+1) - + Ex(j+1,k+1,l-1) - Ex(j+1,k ,l-1) - + Ex(j-1,k+1,l-1) - Ex(j-1,k ,l-1)); -#elif defined WARPX_DIM_XZ - Bz(j,k,0) += - alphax * (Ey(j+1,k ,0) - Ey(j,k ,0)) - - betaxz * (Ey(j+1,k+1,0) - Ey(j,k+1,0) - + Ey(j+1,k-1,0) - Ey(j,k-1,0)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ex_f_ckc(int j, int k, int l, amrex::Array4<amrex::Real> const& Ex, amrex::Array4<amrex::Real const> const& F, diff --git a/Source/FieldSolver/WarpX_K.H b/Source/FieldSolver/WarpX_K.H index a4f657909..226c0abc6 100644 --- a/Source/FieldSolver/WarpX_K.H +++ b/Source/FieldSolver/WarpX_K.H @@ -10,53 +10,6 @@ #include <AMReX_FArrayBox.H> AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_push_bx_nodal (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bx, - amrex::Array4<amrex::Real const> const& Ey, - amrex::Array4<amrex::Real const> const& Ez, - amrex::Real dtsdy, amrex::Real dtsdz) -{ -#if (AMREX_SPACEDIM == 3) - Bx(j,k,l) = Bx(j,k,l) - 0.5*dtsdy * (Ez(j,k+1,l ) - Ez(j,k-1,l )) - + 0.5*dtsdz * (Ey(j,k ,l+1) - Ey(j,k ,l-1)); -#else - Bx(j,k,0) = Bx(j,k,0) + 0.5*dtsdz * (Ey(j,k+1,0) - Ey(j,k-1,0)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_push_by_nodal (int j, int k, int l, - amrex::Array4<amrex::Real> const& By, - amrex::Array4<amrex::Real const> const& Ex, - amrex::Array4<amrex::Real const> const& Ez, - amrex::Real dtsdx, - amrex::Real dtsdz) -{ -#if (AMREX_SPACEDIM == 3) - By(j,k,l) = By(j,k,l) + 0.5*dtsdx * (Ez(j+1,k,l ) - Ez(j-1,k,l )) - - 0.5*dtsdz * (Ex(j ,k,l+1) - Ex(j ,k,l-1)); -#else - By(j,k,0) = By(j,k,0) + 0.5*dtsdx * (Ez(j+1,k ,0) - Ez(j-1,k ,0)) - - 0.5*dtsdz * (Ex(j ,k+1,0) - Ex(j ,k-1,0)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_push_bz_nodal (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bz, - amrex::Array4<amrex::Real const> const& Ex, - amrex::Array4<amrex::Real const> const& Ey, - amrex::Real dtsdx, amrex::Real dtsdy) -{ -#if (AMREX_SPACEDIM == 3) - Bz(j,k,l) = Bz(j,k,l) - 0.5*dtsdx * (Ey(j+1,k ,l) - Ey(j-1,k ,l)) - + 0.5*dtsdy * (Ex(j ,k+1,l) - Ex(j ,k-1,l)); -#else - Bz(j,k,0) = Bz(j,k,0) - 0.5*dtsdx * (Ey(j+1,k ,0) - Ey(j-1,k ,0)); -#endif -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ex_nodal (int j, int k, int l, amrex::Array4<amrex::Real> const& Ex, amrex::Array4<amrex::Real const> const& By, |