diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
5 files changed, 140 insertions, 7 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index 4b1388f42..2455d8907 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -7,6 +7,7 @@ target_sources(WarpX EvolveEPML.cpp EvolveF.cpp EvolveFPML.cpp + EvolveG.cpp FiniteDifferenceSolver.cpp MacroscopicEvolveE.cpp ApplySilverMuellerBoundary.cpp diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index 3ddfc66b0..3d2dc42d8 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -25,6 +25,7 @@ using namespace amrex; void FiniteDifferenceSolver::EvolveB ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, + std::unique_ptr<amrex::MultiFab> const& Gfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, int lev, amrex::Real const dt ) { @@ -38,21 +39,20 @@ void FiniteDifferenceSolver::EvolveB ( #else if (m_do_nodal) { - EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, face_areas, lev, dt ); + EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { - EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, face_areas, lev, dt ); + EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { - EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, face_areas, lev, dt ); + EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt ); #endif } else { amrex::Abort("EvolveB: Unknown algorithm"); } - } @@ -62,6 +62,7 @@ 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, + std::unique_ptr<amrex::MultiFab> const& Gfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, int lev, amrex::Real const dt ) { @@ -71,6 +72,8 @@ void FiniteDifferenceSolver::EvolveBCartesian ( amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev); + Real constexpr c2 = PhysConst::c * PhysConst::c; + // Loop through the grids, and over the tiles within each grid #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -138,9 +141,32 @@ void FiniteDifferenceSolver::EvolveBCartesian ( 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); } - ); + // div(B) cleaning correction for errors in magnetic Gauss law (div(B) = 0) + if (Gfield) + { + // Extract field data for this grid/tile + Array4<Real> G = Gfield->array(mfi); + + // Loop over cells and update G + amrex::ParallelFor(tbx, tby, tbz, + + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Bx(i,j,k) += c2 * dt * T_Algo::DownwardDx(G, coefs_x, n_coefs_x, i, j, k); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + By(i,j,k) += c2 * dt * T_Algo::DownwardDy(G, coefs_y, n_coefs_y, i, j, k); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Bz(i,j,k) += c2 * dt * T_Algo::DownwardDz(G, coefs_z, n_coefs_z, i, j, k); + } + ); + } + if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) { amrex::Gpu::synchronize(); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveG.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveG.cpp new file mode 100644 index 000000000..f69fb57c4 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveG.cpp @@ -0,0 +1,95 @@ +/* Copyright 2020 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "Utils/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 "Utils/WarpXConst.H" +#include "WarpX.H" +#include <AMReX_Gpu.H> + +using namespace amrex; + +void FiniteDifferenceSolver::EvolveG ( + std::unique_ptr<amrex::MultiFab>& Gfield, + std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield, + amrex::Real const dt) +{ +#ifdef WARPX_DIM_RZ + // TODO Implement G update equation in RZ geometry + amrex::ignore_unused(Gfield, Bfield, dt); +#else + // Select algorithm + if (m_do_nodal) + { + EvolveGCartesian<CartesianNodalAlgorithm>(Gfield, Bfield, dt); + } + else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) + { + EvolveGCartesian<CartesianYeeAlgorithm>(Gfield, Bfield, dt); + } + else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) + { + EvolveGCartesian<CartesianCKCAlgorithm>(Gfield, Bfield, dt); + } + else + { + amrex::Abort("EvolveG: unknown FDTD algorithm"); + } +#endif +} + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +void FiniteDifferenceSolver::EvolveGCartesian ( + std::unique_ptr<amrex::MultiFab>& Gfield, + std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield, + amrex::Real const dt) +{ +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + + // Loop over grids and over tiles within each grid + for (amrex::MFIter mfi(*Gfield, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Extract field data for this grid/tile + amrex::Array4<amrex::Real> const& G = Gfield->array(mfi); + amrex::Array4<amrex::Real> const& Bx = Bfield[0]->array(mfi); + amrex::Array4<amrex::Real> const& By = Bfield[1]->array(mfi); + amrex::Array4<amrex::Real> const& Bz = Bfield[2]->array(mfi); + + // Extract stencil coefficients + amrex::Real const* const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + amrex::Real const* const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + amrex::Real const* const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + + const int n_coefs_x = m_stencil_coefs_x.size(); + const int n_coefs_y = m_stencil_coefs_y.size(); + const int n_coefs_z = m_stencil_coefs_z.size(); + + // Extract tilebox to loop over + amrex::Box const& tf = mfi.tilebox(Gfield->ixType().toIntVect()); + + // Loop over cells and update G + amrex::ParallelFor(tf, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + G(i,j,k) += dt * (T_Algo::UpwardDx(Bx, coefs_x, n_coefs_x, i, j, k) + + T_Algo::UpwardDy(By, coefs_y, n_coefs_y, i, j, k) + + T_Algo::UpwardDz(Bz, coefs_z, n_coefs_z, i, j, k)); + }); + } +} + +#endif diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index f434ad3d0..92812d4e0 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -38,6 +38,7 @@ class FiniteDifferenceSolver void EvolveB ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, + std::unique_ptr<amrex::MultiFab> const& Gfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, int lev, amrex::Real const dt ); @@ -54,6 +55,10 @@ class FiniteDifferenceSolver int const rhocomp, amrex::Real const dt ); + void EvolveG (std::unique_ptr<amrex::MultiFab>& Gfield, + std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield, + amrex::Real const dt); + void ApplySilverMuellerBoundary( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield, std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, @@ -146,12 +151,12 @@ class FiniteDifferenceSolver void ComputeDivECylindrical ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); - #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, + std::unique_ptr<amrex::MultiFab> const& Gfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, int lev, amrex::Real const dt ); @@ -173,6 +178,12 @@ class FiniteDifferenceSolver amrex::Real const dt ); template< typename T_Algo > + void EvolveGCartesian ( + std::unique_ptr<amrex::MultiFab>& Gfield, + std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield, + amrex::Real const dt); + + template< typename T_Algo > void ComputeDivECartesian ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); @@ -205,7 +216,6 @@ class FiniteDifferenceSolver void EvolveFPMLCartesian ( amrex::MultiFab* Ffield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); - #endif }; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index ea94f2c66..e1ddee9ac 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -2,6 +2,7 @@ CEXE_sources += FiniteDifferenceSolver.cpp CEXE_sources += EvolveB.cpp CEXE_sources += EvolveE.cpp CEXE_sources += EvolveF.cpp +CEXE_sources += EvolveG.cpp CEXE_sources += ComputeDivE.cpp CEXE_sources += MacroscopicEvolveE.cpp CEXE_sources += EvolveBPML.cpp |