diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
3 files changed, 73 insertions, 69 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index a1f35f903..c08dee04f 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -110,10 +110,9 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield, 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 dr = m_dr; + 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(); @@ -126,12 +125,13 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield, 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; + 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, inv_dz, i, j, 0, 2*m-1) - - m * T_Algo::DivideByR(Ez, r, dr, i, j, 0, 2*m )); // Real part + - 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, inv_dz, i, j, 0, 2*m ) - + m * T_Algo::DivideByR(Ez, r, dr, i, j, 0, 2*m-1)); // Imaginary part + + 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 @@ -144,13 +144,24 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield, }, [=] 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); + Bt(i, j, 0, 0) += dt*( T_Algo::UpwardDr(Ez, inv_dr, i, j, 0, 0) + - T_Algo::UpwardDz(Er, inv_dz, 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, inv_dr, i, j, 0, 2*m-1) + - T_Algo::UpwardDz(Er, inv_dz, i, j, 0, 2*m-1)); // Real part + Bt(i, j, 0, 2*m ) += dt*( T_Algo::UpwardDr(Ez, inv_dr, i, j, 0, 2*m ) + - T_Algo::UpwardDz(Er, inv_dz, i, j, 0, 2*m )); // Imaginary part + } }, [=] 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); + 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, 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, 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, i, j, 0, 2*m )); // Imaginary part } ); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H index 3f408479c..9258fdff4 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H @@ -5,97 +5,77 @@ #include <AMReX_Array4.H> #include <AMReX_Gpu.H> -struct YeeAlgorithm { +struct CylindricalYeeAlgorithm { 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_r, 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_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]; + stencil_coefs_z[0] = 1./cell_size[2]; // 1./dz } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static amrex::Real UpwardDx( + static amrex::Real UpwardDr( 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 const* coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { - amrex::Real inv_dx = coefs_x[0]; - return inv_dx*( F(i+1,j,k) - F(i,j,k) ); + 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 DownwardDx( + static amrex::Real DownwardDr( 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 const* coefs_r, int const n_coefs_r, + int const i, int const j, int const k, int const comp ) { - 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::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 i, int const j, int const k, int const comp ) { 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 + 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 i, int const j, int const k, int const comp ) { 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 + 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; + } + } }; }; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index c08810f85..deb89b1d6 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -27,12 +27,19 @@ class FiniteDifferenceSolver m_fdtd_algo = fdtd_algo; // Calculate coefficients of finite-difference stencil +#ifdef WARPX_DIM_RZ + m_dr = cell_size[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"); } @@ -45,9 +52,15 @@ class FiniteDifferenceSolver int m_fdtd_algo; +#ifdef WARPX_DIM_RZ + amrex::Real m_dr; + 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 > |