From b7d37ccc3f1164449ec5681d45e8415a9c039db7 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 16 Mar 2020 19:15:23 -0700 Subject: Use C++ templates for the `EvolveF` and `ComputeDivE` function (#753) * Prepare EvolveE * Cartesian equations without current * Implement Cartesian EvolveE * Progress towards cylindrical solver * Correct typo * Implement cylindrical solver (without on-axis condition) * Fix compilation errors * Add regularization for RZ solver * Added correction term for F * Remove file for nodal stencil * Apply stylistic changes to EvolveE * Fix compilation errors * Correction to avoid out of bound * Remove references to old file * Correct bug in EvolveB * Implement correction on axis for Et * Remove previous field update functions * Implement EvolveF in Cartesian * Add missing file * Update on-axis condition * Correct typo * Add theta derivative for divE calculation * Fix rho component * Fix compilation error * Use C++ templates in computation of divE * Update Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp Co-Authored-By: MaxThevenet * Apply suggestions from code review Co-Authored-By: MaxThevenet * Use comparisons to 0 again Co-authored-by: MaxThevenet --- .../FiniteDifferenceSolver/ComputeDivE.cpp | 176 +++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp (limited to 'Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp') diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp new file mode 100644 index 000000000..57d241659 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp @@ -0,0 +1,176 @@ +/* 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 "WarpXConst.H" +#include "WarpX.H" +#include + +using namespace amrex; + +/** + * \brief Update the F field, over one timestep + */ +void FiniteDifferenceSolver::ComputeDivE ( + const std::array,3>& Efield, + amrex::MultiFab& divEfield ) { + + // 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){ + + ComputeDivECylindrical ( Efield, divEfield ); + +#else + if (m_do_nodal) { + + ComputeDivECartesian ( Efield, divEfield ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + ComputeDivECartesian ( Efield, divEfield ); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { + + ComputeDivECartesian ( Efield, divEfield ); + +#endif + } else { + amrex::Abort("Unknown algorithm"); + } + +} + + +#ifndef WARPX_DIM_RZ + +template +void FiniteDifferenceSolver::ComputeDivECartesian ( + const std::array,3>& Efield, + amrex::MultiFab& divEfield ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(divEfield, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4 const& divE = divEfield.array(mfi); + Array4 const& Ex = Efield[0]->array(mfi); + Array4 const& Ey = Efield[1]->array(mfi); + Array4 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& tdive = mfi.tilebox(divEfield.ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tdive, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + divE(i, j, k) = + T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k) + + T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k) + + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k); + } + + ); + + } + +} + +#else // corresponds to ifndef WARPX_DIM_RZ + +template +void FiniteDifferenceSolver::ComputeDivECylindrical ( + const std::array,3>& Efield, + amrex::MultiFab& divEfield ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(divEfield, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4 divE = divEfield.array(mfi); + Array4 const& Er = Efield[0]->array(mfi); + Array4 const& Et = Efield[1]->array(mfi); + Array4 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& tdive = mfi.tilebox(divEfield.ixType().ixType()); + + // Loop over the cells and update the fields + amrex::ParallelFor(tdive, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const r = rmin + i*dr; // r on a nodal grid (F is nodal in r) + if (r != 0) { // Off-axis, regular equations + divE(i, j, 0, 0) = + T_Algo::DownwardDrr_over_r(Er, r, dr, coefs_r, n_coefs_r, i, j, 0, 0) + + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, 0, 0); + for (int m=1 ; m