From 6d0438121539950f8c21a556443307965661bf51 Mon Sep 17 00:00:00 2001 From: roelof-groenewald Date: Tue, 27 Aug 2024 11:22:12 -0700 Subject: [PATCH] major refactor of magnetostatic solver --- Python/pywarpx/fields.py | 6 +- .../FiniteDifferenceSolver/EvolveB.cpp | 188 ++++++++++++++ .../FiniteDifferenceSolver.H | 32 +++ .../FiniteDifferenceSolver.cpp | 8 +- .../MagnetostaticSolver/MagnetostaticSolver.H | 30 +-- .../MagnetostaticSolver.cpp | 230 +++++++----------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 21 ++ Source/Initialization/WarpXInitData.cpp | 6 +- Source/WarpX.H | 10 +- Source/WarpX.cpp | 40 +-- 10 files changed, 360 insertions(+), 211 deletions(-) diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py index b765424cf57..b56c3407ad2 100644 --- a/Python/pywarpx/fields.py +++ b/Python/pywarpx/fields.py @@ -701,7 +701,7 @@ def GFPWrapper(level=0, include_ghosts=False): def AxFPWrapper(level=0, include_ghosts=False): return _MultiFABWrapper( - mf_name="vector_potential_fp_nodal[x]", + mf_name="Afield_fp_nodal[x]", level=level, include_ghosts=include_ghosts, ) @@ -709,7 +709,7 @@ def AxFPWrapper(level=0, include_ghosts=False): def AyFPWrapper(level=0, include_ghosts=False): return _MultiFABWrapper( - mf_name="vector_potential_fp_nodal[y]", + mf_name="Afield_fp_nodal[y]", level=level, include_ghosts=include_ghosts, ) @@ -717,7 +717,7 @@ def AyFPWrapper(level=0, include_ghosts=False): def AzFPWrapper(level=0, include_ghosts=False): return _MultiFABWrapper( - mf_name="vector_potential_fp_nodal[z]", + mf_name="Afield_fp_nodal[z]", level=level, include_ghosts=include_ghosts, ) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index fbc1397b413..53c865410c3 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -484,3 +484,191 @@ void FiniteDifferenceSolver::EvolveBCylindrical ( } #endif // corresponds to ifndef WARPX_DIM_RZ + + +void FiniteDifferenceSolver::CalcBfromVectorPotential ( + amrex::Array, 3>& Bfield, + amrex::Array, 3> const& Afield, + std::array, 3> const& edge_lengths, + const int lev ) +{ + if (m_fdtd_algo != ElectromagneticSolverAlgo::CKC && m_fdtd_algo != ElectromagneticSolverAlgo::PSATD) { +#ifdef WARPX_DIM_RZ + CalcBfromVectorPotentialCylindrical ( + Bfield, Afield, edge_lengths, lev + ); +#else + CalcBfromVectorPotentialCartesian ( + Bfield, Afield, edge_lengths, lev + ); +#endif + } else { + WARPX_ABORT_WITH_MESSAGE("CalcBfromVectorPotential: Unknown algorithm"); + } +} + +#ifdef WARPX_DIM_RZ +template +void FiniteDifferenceSolver::CalcBfromVectorPotentialCylindrical ( + amrex::Array, 3>& Bfield, + amrex::Array, 3> const& Afield, + [[maybe_unused]] std::array< std::unique_ptr, 3 > const& edge_lengths, + [[maybe_unused]] const int lev +) +{ + // Loop through the grids, and over the tiles within each grid +#ifdef AMREX_USE_OMP +#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 const& Br = Bfield[0]->array(mfi); + Array4 const& Bt = Bfield[1]->array(mfi); + Array4 const& Bz = Bfield[2]->array(mfi); + Array4 const& Ar = Afield[0]->array(mfi); + Array4 const& At = Afield[1]->array(mfi); + Array4 const& Az = Afield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_r = m_stencil_coefs_r.dataPtr(); + auto const n_coefs_r = static_cast(m_stencil_coefs_r.size()); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + auto const n_coefs_z = static_cast(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().toIntVect()); + Box const& tbt = mfi.tilebox(Bfield[1]->ixType().toIntVect()); + Box const& tbz = mfi.tilebox(Bfield[2]->ixType().toIntVect()); + + // 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) += T_Algo::UpwardDz(At, coefs_z, n_coefs_z, i, j, 0, 0); // Mode m=0 + for (int m=1; m +void FiniteDifferenceSolver::CalcBfromVectorPotentialCartesian ( + amrex::Array, 3>& Bfield, + amrex::Array, 3> const& Afield, + [[maybe_unused]] std::array< std::unique_ptr, 3 > const& edge_lengths, + [[maybe_unused]] const int lev +) +{ +// Loop through the grids, and over the tiles within each grid +#ifdef AMREX_USE_OMP +#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 const& Bx = Bfield[0]->array(mfi); + Array4 const& By = Bfield[1]->array(mfi); + Array4 const& Bz = Bfield[2]->array(mfi); + Array4 const& Ax = Afield[0]->array(mfi); + Array4 const& Ay = Afield[1]->array(mfi); + Array4 const& Az = Afield[2]->array(mfi); + + // Extract stencil coefficients + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + auto const n_coefs_x = static_cast(m_stencil_coefs_x.size()); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + auto const n_coefs_y = static_cast(m_stencil_coefs_y.size()); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + auto const n_coefs_z = static_cast(m_stencil_coefs_z.size()); + + // Extract tileboxes for which to loop + Box const& tbx = mfi.tilebox(Bfield[0]->ixType().toIntVect()); + Box const& tby = mfi.tilebox(Bfield[1]->ixType().toIntVect()); + Box const& tbz = mfi.tilebox(Bfield[2]->ixType().toIntVect()); + + // 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) += T_Algo::UpwardDz(Ay, coefs_z, n_coefs_z, i, j, k) + - T_Algo::UpwardDy(Az, coefs_y, n_coefs_y, i, j, k); + + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + + By(i, j, k) += T_Algo::UpwardDx(Az, coefs_x, n_coefs_x, i, j, k) + - T_Algo::UpwardDz(Ax, coefs_z, n_coefs_z, i, j, k); + + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + + Bz(i, j, k) += T_Algo::UpwardDy(Ax, coefs_y, n_coefs_y, i, j, k) + - T_Algo::UpwardDx(Ay, coefs_x, n_coefs_x, i, j, k); + + } + ); + } +} +#endif // corresponds to ifndef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 00e87525a75..86d8276a648 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -177,6 +177,22 @@ class FiniteDifferenceSolver std::array< std::unique_ptr, 3 > const& edge_lengths, int lev ); + /** + * \brief Calculation of the magnetic field given a vector potential + * field from B = ∇ x A. + * + * \param[out] Bfield vector of B-field MultiFabs at a given level + * \param[in] Afield vector of vector potential MultiFabs at a given level + * \param[in] edge_lengths length of edges along embedded boundaries + * \param[in] lev level number for the calculation + */ + void CalcBfromVectorPotential ( + amrex::Array, 3>& Bfield, + amrex::Array, 3> const& Afield, + std::array< std::unique_ptr, 3 > const& edge_lengths, + int lev + ); + private: int m_fdtd_algo; @@ -255,6 +271,14 @@ class FiniteDifferenceSolver int lev ); + template + void CalcBfromVectorPotentialCylindrical ( + amrex::Array, 3>& Bfield, + amrex::Array, 3> const& Afield, + std::array< std::unique_ptr, 3 > const& edge_lengths, + int lev + ); + #else template< typename T_Algo > void EvolveBCartesian ( @@ -359,6 +383,14 @@ class FiniteDifferenceSolver std::array< std::unique_ptr, 3 > const& edge_lengths, int lev ); + + template + void CalcBfromVectorPotentialCartesian ( + amrex::Array, 3>& Bfield, + amrex::Array, 3> const& Afield, + std::array< std::unique_ptr, 3 > const& edge_lengths, + int lev + ); #endif }; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp index 9af610031c0..f7034dcb114 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp @@ -36,7 +36,7 @@ FiniteDifferenceSolver::FiniteDifferenceSolver ( m_grid_type{grid_type} { // return if not FDTD - if (fdtd_algo == ElectromagneticSolverAlgo::None || fdtd_algo == ElectromagneticSolverAlgo::PSATD) { + if (fdtd_algo == ElectromagneticSolverAlgo::PSATD) { return; } @@ -45,7 +45,8 @@ FiniteDifferenceSolver::FiniteDifferenceSolver ( m_dr = cell_size[0]; m_nmodes = WarpX::n_rz_azimuthal_modes; m_rmin = WarpX::GetInstance().Geom(0).ProbLo(0); - if (fdtd_algo == ElectromagneticSolverAlgo::Yee || + if (fdtd_algo == ElectromagneticSolverAlgo::None || + fdtd_algo == ElectromagneticSolverAlgo::Yee || fdtd_algo == ElectromagneticSolverAlgo::HybridPIC ) { CylindricalYeeAlgorithm::InitializeStencilCoefficients( cell_size, m_h_stencil_coefs_r, m_h_stencil_coefs_z ); @@ -68,7 +69,8 @@ FiniteDifferenceSolver::FiniteDifferenceSolver ( CartesianNodalAlgorithm::InitializeStencilCoefficients( cell_size, m_h_stencil_coefs_x, m_h_stencil_coefs_y, m_h_stencil_coefs_z ); - } else if (fdtd_algo == ElectromagneticSolverAlgo::Yee || + } else if (fdtd_algo == ElectromagneticSolverAlgo::None || + fdtd_algo == ElectromagneticSolverAlgo::Yee || fdtd_algo == ElectromagneticSolverAlgo::ECT || fdtd_algo == ElectromagneticSolverAlgo::HybridPIC) { diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H index a8bbc954e29..cbde5cdd56c 100644 --- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H +++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H @@ -1,7 +1,9 @@ -/* Copyright 2022 S. Eric Clark, LLNL +/* Copyright 2022 The WarpX Community * * This file is part of WarpX. * + * Authors: S. Eric Clark (LLNL) + * * License: BSD-3-Clause-LBNL */ #ifndef WARPX_MAGNETOSTATICSOLVER_H_ @@ -12,6 +14,7 @@ #include #include + namespace MagnetostaticSolver { /** Boundary Handler for the Vector Potential Poisson Solver @@ -27,31 +30,6 @@ namespace MagnetostaticSolver { void defineVectorPotentialBCs (); }; - - /** use amrex to directly calculate the magnetic field since with EB's the - * - * simple finite difference scheme in WarpX::computeE sometimes fails - */ - class EBCalcBfromVectorPotentialPerLevel { - private: - const amrex::Vector, 3>>& m_b_field; - const amrex::Vector, 3>>& m_grad_buf_e_stag; - const amrex::Vector, 3>>& m_grad_buf_b_stag; - - public: - EBCalcBfromVectorPotentialPerLevel(const amrex::Vector, 3>>& b_field, - const amrex::Vector, 3>>& grad_buf_e_stag, - const amrex::Vector, 3>>& grad_buf_b_stag) - : m_b_field(b_field), - m_grad_buf_e_stag(grad_buf_e_stag), - m_grad_buf_b_stag(grad_buf_b_stag) - {} - - void operator()(amrex::Array,3> & mlmg, int lev); - - // Function to perform interpolation from cell edges to cell faces - void doInterp(const std::unique_ptr &src, const std::unique_ptr &dst); - }; } // namespace MagnetostaticSolver #endif //WARPX_MAGNETOSTATICSOLVER_H_ diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp index c0d6d0ff6e3..76193f339d0 100644 --- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp +++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp @@ -1,7 +1,9 @@ -/* Copyright 2022 S. Eric Clark, LLNL +/* Copyright 2022-2024 The WarpX Community * * This file is part of WarpX. * + * Authors: S. Eric Clark (LLNL), Roelof Groenewald (TAE Technologies) + * * License: BSD-3-Clause-LBNL */ #include "WarpX.H" @@ -9,9 +11,7 @@ #include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H" #include "Parallelization/GuardCellManager.H" #include "Particles/MultiParticleContainer.H" -#include "Particles/WarpXParticleContainer.H" #include "Python/callbacks.H" -#include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" #include "Utils/TextMsg.H" #include "Utils/WarpXUtil.H" @@ -62,54 +62,63 @@ WarpX::ComputeMagnetostaticField() // Fields have been reset in Electrostatic solver for this time step, these fields // are added into the B fields after electrostatic solve - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(this->max_level == 0, "Magnetostatic solver not implemented with mesh refinement."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + this->max_level == 0, + "Magnetostatic solver not implemented with mesh refinement." + ); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + do_magnetostatic_solve, + "Magnetostatic solver called but warpx.do_magnetostatic is false." + ); - AddMagnetostaticFieldLabFrame(); + ComputeMagnetostaticFieldLabFrame(); } void -WarpX::AddMagnetostaticFieldLabFrame() +WarpX::ComputeMagnetostaticFieldLabFrame() { - WARPX_PROFILE("WarpX::AddMagnetostaticFieldLabFrame"); - - // Store the boundary conditions for the field solver if they haven't been - // stored yet - if (!m_vector_poisson_boundary_handler.bcs_set) { - m_vector_poisson_boundary_handler.defineVectorPotentialBCs(); - } + WARPX_PROFILE("WarpX::ComputeMagnetostaticFieldLabFrame"); #ifdef WARPX_DIM_RZ WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, "Error: RZ magnetostatic only implemented for a single mode"); #endif - // reset current_fp before depositing current density for this step - for (int lev = 0; lev <= max_level; lev++) { - for (int dim=0; dim < 3; dim++) { - current_fp[lev][dim]->setVal(0.); + // Perform current deposition at t_{n+1/2} (source of Poisson solver). + mypc->DepositCurrent(current_fp, dt[0], -0.5_rt * dt[0]); + + // Synchronize J and rho: + // filter (if used), exchange guard cells, interpolate across MR levels + // and apply boundary conditions + SyncCurrent(current_fp, current_cp, current_buf); + for (int lev = 0; lev <= finest_level; ++lev) + { + ApplyJfieldBoundary( + lev, current_fp[lev][0].get(), current_fp[lev][1].get(), + current_fp[lev][2].get(), PatchType::fine + ); + if (lev > 0) { + ApplyJfieldBoundary( + lev, current_cp[lev][0].get(), current_cp[lev][1].get(), + current_cp[lev][2].get(), PatchType::coarse + ); } } - // Deposit current density (source of Poisson solver) - for (int ispecies=0; ispeciesnSpecies(); ispecies++){ - WarpXParticleContainer& species = mypc->GetParticleContainer(ispecies); - if (!species.do_not_deposit) { - species.DepositCurrent(current_fp, dt[0], -0.5_rt*dt[0]); + // SyncCurrent does not include a call to FillBoundary + for (int lev = 0; lev <= finest_level; ++lev) { + for (int idim = 0; idim < 3; ++idim) { + current_fp[lev][idim]->FillBoundary(Geom(lev).periodicity()); } } -#ifdef WARPX_DIM_RZ - for (int lev = 0; lev <= max_level; lev++) { - ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), - current_fp[lev][1].get(), - current_fp[lev][2].get(), lev); + // Store the boundary conditions for the field solver if they haven't been + // stored yet + if (!m_vector_poisson_boundary_handler.bcs_set) { + m_vector_poisson_boundary_handler.defineVectorPotentialBCs(); } -#endif - - SyncCurrent(current_fp, current_cp, current_buf); // Apply filter, perform MPI exchange, interpolate across levels - // set the boundary and current density potentials - setVectorPotentialBC(vector_potential_fp_nodal); + setVectorPotentialBC(Afield_fp_nodal); // Compute the vector potential A, by solving the Poisson equation WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !IsPythonCallbackInstalled("poissonsolver"), @@ -117,9 +126,50 @@ WarpX::AddMagnetostaticFieldLabFrame() const amrex::Real magnetostatic_absolute_tolerance = self_fields_absolute_tolerance*PhysConst::c; - computeVectorPotential( current_fp, vector_potential_fp_nodal, self_fields_required_precision, - magnetostatic_absolute_tolerance, self_fields_max_iters, - self_fields_verbosity ); + computeVectorPotential( + current_fp, Afield_fp_nodal, + self_fields_required_precision, magnetostatic_absolute_tolerance, + self_fields_max_iters, self_fields_verbosity + ); + + // Grab Interpolation Coefficients + // Order of finite-order centering of fields + const int fg_nox = WarpX::field_centering_nox; + const int fg_noy = WarpX::field_centering_noy; + const int fg_noz = WarpX::field_centering_noz; + + // Device vectors of stencil coefficients used for finite-order centering of fields + amrex::Real const * stencil_coeffs_x = WarpX::device_field_centering_stencil_coeffs_x.data(); + amrex::Real const * stencil_coeffs_y = WarpX::device_field_centering_stencil_coeffs_y.data(); + amrex::Real const * stencil_coeffs_z = WarpX::device_field_centering_stencil_coeffs_z.data(); + + for (int lev = 0; lev <= finest_level; ++lev) + { + // interpolate Afield_fp_nodal[lev][i] to Afield_fp[lev][i] + for (int i = 0; i < 3; ++i) { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Afield_fp[lev][i], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + IntVect const src_stag = Afield_fp_nodal[lev][i]->ixType().toIntVect(); + IntVect const dst_stag = Afield_fp[lev][i]->ixType().toIntVect(); + + Array4 const& src_arr = Afield_fp_nodal[lev][i]->const_array(mfi); + Array4 const& dst_arr = Afield_fp[lev][i]->array(mfi); + + const Box bx = mfi.tilebox(); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, dst_arr, src_arr, dst_stag, src_stag, fg_nox, fg_noy, fg_noz, + stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); + }); + } + } + } + + ComputeBfromVectorPotential(); } /* Compute the vector potential `A` by solving the Poisson equation with `J` as @@ -167,9 +217,6 @@ WarpX::computeVectorPotential (const amrex::Vector > eb_farray_box_factory; #endif - const std::optional post_A_calculation( - {Bfield_fp, vector_potential_grad_buf_e_stag, vector_potential_grad_buf_b_stag} - ); ablastr::fields::computeVectorPotential( sorted_curr, @@ -184,11 +231,10 @@ WarpX::computeVectorPotential (const amrex::Vectorm_vector_poisson_boundary_handler, WarpX::do_single_precision_comms, this->ref_ratio, - post_A_calculation, + std::nullopt, // post_A_calculation, gett_new(0), eb_farray_box_factory ); - } @@ -346,107 +392,3 @@ void MagnetostaticSolver::VectorPoissonBoundaryHandler::defineVectorPotentialBCs } bcs_set = true; } - -void MagnetostaticSolver::EBCalcBfromVectorPotentialPerLevel::doInterp(const std::unique_ptr &src, - const std::unique_ptr &dst) -{ - WarpX &warpx = WarpX::GetInstance(); - - // Grab Interpolation Coefficients - // Order of finite-order centering of fields - const int fg_nox = WarpX::field_centering_nox; - const int fg_noy = WarpX::field_centering_noy; - const int fg_noz = WarpX::field_centering_noz; - - // Device vectors of stencil coefficients used for finite-order centering of fields - amrex::Real const * stencil_coeffs_x = warpx.device_field_centering_stencil_coeffs_x.data(); - amrex::Real const * stencil_coeffs_y = warpx.device_field_centering_stencil_coeffs_y.data(); - amrex::Real const * stencil_coeffs_z = warpx.device_field_centering_stencil_coeffs_z.data(); - - // Synchronize the ghost cells, do halo exchange - ablastr::utils::communication::FillBoundary(*src, - src->nGrowVect(), - WarpX::do_single_precision_comms); - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(*dst, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - IntVect const src_stag = src->ixType().toIntVect(); - IntVect const dst_stag = dst->ixType().toIntVect(); - - Array4 const& src_arr = src->const_array(mfi); - Array4 const& dst_arr = dst->array(mfi); - - const Box bx = mfi.tilebox(); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept - { - warpx_interp(j, k, l, dst_arr, src_arr, dst_stag, src_stag, fg_nox, fg_noy, fg_noz, - stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); - }); - } -} - -void MagnetostaticSolver::EBCalcBfromVectorPotentialPerLevel::operator()(amrex::Array,3> & mlmg, int const lev) -{ - using namespace amrex::literals; - - // This operator gets the gradient solution on the cell edges, aligned with E field staggered grid - // This routine interpolates to the B-field staggered grid, - - const amrex::Array buf_ptr = - { -#if defined(WARPX_DIM_3D) - m_grad_buf_e_stag[lev][0].get(), - m_grad_buf_e_stag[lev][1].get(), - m_grad_buf_e_stag[lev][2].get() -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - m_grad_buf_e_stag[lev][0].get(), - m_grad_buf_e_stag[lev][2].get() -#endif - }; - - // This will grab the gradient values for Ax - mlmg[0]->getGradSolution({buf_ptr}); - - // Interpolate dAx/dz to By grid buffer, then add to By - this->doInterp(m_grad_buf_e_stag[lev][2], - m_grad_buf_b_stag[lev][1]); - MultiFab::Add(*(m_b_field[lev][1]), *(m_grad_buf_b_stag[lev][1]), 0, 0, 1, 0 ); - - // Interpolate dAx/dy to Bz grid buffer, then subtract from Bz - this->doInterp(m_grad_buf_e_stag[lev][1], - m_grad_buf_b_stag[lev][2]); - m_grad_buf_b_stag[lev][2]->mult(-1._rt); - MultiFab::Add(*(m_b_field[lev][2]), *(m_grad_buf_b_stag[lev][2]), 0, 0, 1, 0 ); - - // This will grab the gradient values for Ay - mlmg[1]->getGradSolution({buf_ptr}); - - // Interpolate dAy/dx to Bz grid buffer, then add to Bz - this->doInterp(m_grad_buf_e_stag[lev][0], - m_grad_buf_b_stag[lev][2]); - MultiFab::Add(*(m_b_field[lev][2]), *(m_grad_buf_b_stag[lev][2]), 0, 0, 1, 0 ); - - // Interpolate dAy/dz to Bx grid buffer, then subtract from Bx - this->doInterp(m_grad_buf_e_stag[lev][2], - m_grad_buf_b_stag[lev][0]); - m_grad_buf_b_stag[lev][0]->mult(-1._rt); - MultiFab::Add(*(m_b_field[lev][0]), *(m_grad_buf_b_stag[lev][0]), 0, 0, 1, 0 ); - - // This will grab the gradient values for Az - mlmg[2]->getGradSolution({buf_ptr}); - - // Interpolate dAz/dy to Bx grid buffer, then add to Bx - this->doInterp(m_grad_buf_e_stag[lev][1], - m_grad_buf_b_stag[lev][0]); - MultiFab::Add(*(m_b_field[lev][0]), *(m_grad_buf_b_stag[lev][0]), 0, 0, 1, 0 ); - - // Interpolate dAz/dx to By grid buffer, then subtract from By - this->doInterp(m_grad_buf_e_stag[lev][0], - m_grad_buf_b_stag[lev][1]); - m_grad_buf_b_stag[lev][1]->mult(-1._rt); - MultiFab::Add(*(m_b_field[lev][1]), *(m_grad_buf_b_stag[lev][1]), 0, 0, 1, 0 ); -} diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 602a2666b27..295ccb5a0ea 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -845,6 +845,27 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_typ } +void +WarpX::ComputeBfromVectorPotential () +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + ComputeBfromVectorPotential(lev); + } +} + + +void +WarpX::ComputeBfromVectorPotential (const int lev) +{ + // Calculate B = ∇ x A. + m_fdtd_solver_fp[lev]->CalcBfromVectorPotential( + Bfield_fp[lev], Afield_fp[lev], m_edge_lengths[lev], lev + ); + ApplyBfieldBoundary(lev, PatchType::fine, DtType::Full); +} + + void WarpX::EvolveE (amrex::Real a_dt) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 836fdc25bcb..a722b5e0845 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -269,13 +269,13 @@ WarpX::PrintMainPICparameters () amrex::Print() << " | - laboratory frame, semi-implicit" << "\n"; } if (poisson_solver_id == PoissonSolverAlgo::IntegratedGreenFunction){ - amrex::Print() << "Poisson solver: | FFT-based" << "\n"; + amrex::Print() << "Poisson solver: | - FFT-based" << "\n"; } else if(poisson_solver_id == PoissonSolverAlgo::Multigrid){ - amrex::Print() << "Poisson solver: | multigrid" << "\n"; + amrex::Print() << "Poisson solver: | - multigrid" << "\n"; } if (do_magnetostatic_solve){ - amrex::Print() << " | magnetostatic solve included" << "\n"; + amrex::Print() << " | - magnetostatic solve included" << "\n"; } } else{ diff --git a/Source/WarpX.H b/Source/WarpX.H index 2720355c20d..07414b6f7a5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1008,7 +1008,7 @@ public: // Magnetostatic Solver Interface MagnetostaticSolver::VectorPoissonBoundaryHandler m_vector_poisson_boundary_handler; void ComputeMagnetostaticField (); - void AddMagnetostaticFieldLabFrame (); + void ComputeMagnetostaticFieldLabFrame (); void computeVectorPotential (const amrex::Vector, 3> >& curr, amrex::Vector, 3> >& A, amrex::Real required_precision=amrex::Real(1.e-11), @@ -1018,6 +1018,9 @@ public: void setVectorPotentialBC (amrex::Vector, 3> >& A) const; + void ComputeBfromVectorPotential (); + void ComputeBfromVectorPotential (const int lev); + /** * \brief * This function initializes the E and B fields on each level @@ -1507,9 +1510,8 @@ private: // Memory buffers for computing magnetostatic fields // Vector Potential A and previous step. Time buffer needed for computing dA/dt to first order - amrex::Vector, 3 > > vector_potential_fp_nodal; - amrex::Vector, 3 > > vector_potential_grad_buf_e_stag; - amrex::Vector, 3 > > vector_potential_grad_buf_b_stag; + amrex::Vector, 3 > > Afield_fp; + amrex::Vector, 3 > > Afield_fp_nodal; // Same as Bfield_fp/Efield_fp for reading external field data amrex::Vector, 3 > > Efield_fp_external; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c7403518366..c5b118f92b3 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -331,9 +331,8 @@ WarpX::WarpX () // Only allocate vector potential arrays when using the Magnetostatic Solver if (do_magnetostatic_solve) { - vector_potential_fp_nodal.resize(nlevs_max); - vector_potential_grad_buf_e_stag.resize(nlevs_max); - vector_potential_grad_buf_b_stag.resize(nlevs_max); + Afield_fp_nodal.resize(nlevs_max); + Afield_fp.resize(nlevs_max); } // Initialize the semi-implicit electrostatic solver if required @@ -2103,9 +2102,8 @@ WarpX::ClearLevel (int lev) if (do_magnetostatic_solve) { - vector_potential_fp_nodal[lev][i].reset(); - vector_potential_grad_buf_e_stag[lev][i].reset(); - vector_potential_grad_buf_b_stag[lev][i].reset(); + Afield_fp_nodal[lev][i].reset(); + Afield_fp[lev][i].reset(); } current_cp[lev][i].reset(); @@ -2222,7 +2220,6 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d m_accelerator_lattice[lev] = std::make_unique(); m_accelerator_lattice[lev]->InitElementFinder(lev, ba, dm); - } void @@ -2272,6 +2269,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm jy_nodal_flag = IntVect(1,0,1); jz_nodal_flag = IntVect(1,1,0); #endif + if (do_magnetostatic_solve) { jx_nodal_flag = IntVect::TheNodeVector(); @@ -2360,26 +2358,12 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (do_magnetostatic_solve) { - AllocInitMultiFab(vector_potential_fp_nodal[lev][0], amrex::convert(ba, rho_nodal_flag), - dm, ncomps, ngRho, lev, "vector_potential_fp_nodal[x]", 0.0_rt); - AllocInitMultiFab(vector_potential_fp_nodal[lev][1], amrex::convert(ba, rho_nodal_flag), - dm, ncomps, ngRho, lev, "vector_potential_fp_nodal[y]", 0.0_rt); - AllocInitMultiFab(vector_potential_fp_nodal[lev][2], amrex::convert(ba, rho_nodal_flag), - dm, ncomps, ngRho, lev, "vector_potential_fp_nodal[z]", 0.0_rt); - - AllocInitMultiFab(vector_potential_grad_buf_e_stag[lev][0], amrex::convert(ba, Ex_nodal_flag), - dm, ncomps, ngEB, lev, "vector_potential_grad_buf_e_stag[x]", 0.0_rt); - AllocInitMultiFab(vector_potential_grad_buf_e_stag[lev][1], amrex::convert(ba, Ey_nodal_flag), - dm, ncomps, ngEB, lev, "vector_potential_grad_buf_e_stag[y]", 0.0_rt); - AllocInitMultiFab(vector_potential_grad_buf_e_stag[lev][2], amrex::convert(ba, Ez_nodal_flag), - dm, ncomps, ngEB, lev, "vector_potential_grad_buf_e_stag[z]", 0.0_rt); - - AllocInitMultiFab(vector_potential_grad_buf_b_stag[lev][0], amrex::convert(ba, Bx_nodal_flag), - dm, ncomps, ngEB, lev, "vector_potential_grad_buf_b_stag[x]", 0.0_rt); - AllocInitMultiFab(vector_potential_grad_buf_b_stag[lev][1], amrex::convert(ba, By_nodal_flag), - dm, ncomps, ngEB, lev, "vector_potential_grad_buf_b_stag[y]", 0.0_rt); - AllocInitMultiFab(vector_potential_grad_buf_b_stag[lev][2], amrex::convert(ba, Bz_nodal_flag), - dm, ncomps, ngEB, lev, "vector_potential_grad_buf_b_stag[z]", 0.0_rt); + AllocInitMultiFab(Afield_fp_nodal[lev][0], amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, lev, "Afield_fp_nodal[x]", 0.0_rt); + AllocInitMultiFab(Afield_fp_nodal[lev][1], amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, lev, "Afield_fp_nodal[y]", 0.0_rt); + AllocInitMultiFab(Afield_fp_nodal[lev][2], amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, lev, "Afield_fp_nodal[z]", 0.0_rt); + AllocInitMultiFab(Afield_fp[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, lev, "Afield_fp[x]", 0.0_rt); + AllocInitMultiFab(Afield_fp[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, lev, "Afield_fp[y]", 0.0_rt); + AllocInitMultiFab(Afield_fp[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, lev, "Afield_fp[z]", 0.0_rt); } // Allocate extra multifabs needed by the semi-implicit electrostatic algorithm. @@ -3528,7 +3512,7 @@ WarpX::getFieldPointerUnchecked (const FieldType field_type, const int lev, cons field_pointer = phi_fp[lev].get(); break; case FieldType::vector_potential_fp : - field_pointer = vector_potential_fp_nodal[lev][direction].get(); + field_pointer = Afield_fp_nodal[lev][direction].get(); break; case FieldType::Efield_cp : field_pointer = Efield_cp[lev][direction].get();