diff options
21 files changed, 88 insertions, 1189 deletions
diff --git a/Docs/Doxyfile b/Docs/Doxyfile index 88935e3df..6875f7311 100644 --- a/Docs/Doxyfile +++ b/Docs/Doxyfile @@ -2061,8 +2061,7 @@ PREDEFINED = AMREX_Linux=1 \ WARPX_DIM_XZ=1 \ WARPX_USE_GPU=1 \ WARPX_USE_OPENPMD=1 \ - WARPX_USE_PSATD=1 \ - WARPX_USE_PSATD_HYBRID=1 + WARPX_USE_PSATD=1 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this # tag can be used to specify a list of macro names that should be expanded. The diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 70ea888c2..904741ebf 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -1031,18 +1031,6 @@ Numerics and algorithms Therefore, all the approximations that are usually made when using local FFTs with guard cells (for problems with multiple boxes) become exact in the case of the periodic, single-box FFT without guard cells. -* ``psatd.hybrid_mpi_decomposition`` (`0` or `1`; default: 0) - Whether to use a different MPI decomposition for the particle-grid operations - (deposition and gather) and for the PSATD solver. If `1`, the FFT will - be performed over MPI groups. - -* ``psatd.ngroups_fft`` (`integer`) - The number of MPI groups that are created for the FFT, when using the code compiled with a PSATD solver - (and only if `hybrid_mpi_decomposition` is `1`). - The FFTs are global within one MPI group and use guard cell exchanges in between MPI groups. - (If ``ngroups_fft`` is larger than the number of MPI ranks used, - than the actual number of MPI ranks is used instead.) - * ``psatd.fftw_plan_measure`` (`0` or `1`) Defines whether the parameters of FFTW plans will be initialized by measuring and optimizing performance (``FFTW_MEASURE`` mode; activated by default here). diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index f9c4f0cae..966de58f8 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -431,25 +431,6 @@ analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi.py analysisOutputImage = langmuir_multi_analysis.png tolerance = 5.e-11 -[Langmuir_multi_psatd_hybrid] -buildDir = . -inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt -runtime_params = psatd.fftw_plan_measure=0 psatd.hybrid_mpi_decomposition=1 psatd.ngroups_fft=2 -dim = 3 -addToCompileString = USE_PSATD=TRUE USE_PSATD_PICSAR=TRUE -restartTest = 0 -useMPI = 1 -numprocs = 2 -useOMP = 1 -numthreads = 1 -compileTest = 0 -doVis = 0 -compareParticles = 1 -particleTypes = electrons positrons -analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi.py -analysisOutputImage = langmuir_multi_analysis.png -tolerance = 5.e-11 - [Langmuir_multi_psatd_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index b73447cab..770c7cadd 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -578,8 +578,6 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, const BoxArray& ba = amrex::convert(boxArray(lev),IntVect::TheUnitVector()); MultiFab divE( ba, DistributionMap(lev), 1, 0 ); #ifdef WARPX_USE_PSATD - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( fft_hybrid_mpi_decomposition == false, - "Spectral computation of div(E) not implemented with PICSAR spectral solver" ); spectral_solver_fp[lev]->ComputeSpectralDivE( Efield_aux[lev], divE ); #else m_fdtd_solver_fp[lev]->ComputeDivE( Efield_aux[lev], divE ); diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index dc5342a78..caaac36ff 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -586,7 +586,7 @@ WarpXOpenPMDPlot::WriteOpenPMDFields( //const std::string& filename, auto meshes = series_iteration.meshes; meshes.setAttribute( "fieldSolver", [](){ #ifdef WARPX_USE_PSATD - return "PSATD"; // TODO double-check if WARPX_USE_PSATD_HYBRID is covered + return "PSATD"; #else switch( WarpX::particle_pusher_algo ) { case MaxwellSolverAlgo::Yee : return "Yee"; diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package index 03c2f830a..ba5eb1ac3 100644 --- a/Source/FieldSolver/Make.package +++ b/Source/FieldSolver/Make.package @@ -3,9 +3,6 @@ CEXE_sources += ElectrostaticSolver.cpp CEXE_sources += WarpX_QED_Field_Pushers.cpp ifeq ($(USE_PSATD),TRUE) include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/Make.package - ifeq ($(USE_PSATD_PICSAR),TRUE) - include $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package - endif endif include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/Make.package diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package b/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package deleted file mode 100644 index 22d38025c..000000000 --- a/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package +++ /dev/null @@ -1,4 +0,0 @@ -CEXE_sources += PicsarHybridSpectralSolver.cpp -F90EXE_sources += picsar_hybrid_spectral.F90 - -VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H deleted file mode 100644 index 7e037cf8d..000000000 --- a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H +++ /dev/null @@ -1,43 +0,0 @@ -/* Copyright 2019 Remi Lehe - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef WARPX_PICSAR_HYBRID_FFTDATA_H_ -#define WARPX_PICSAR_HYBRID_FFTDATA_H_ - -#ifdef WARPX_USE_PSATD -# include <fftw3.h> - - -// FFTData is a stuct containing a 22 pointers to auxiliary arrays -// 1-11: padded arrays in real space (required by FFTW); 12-22: arrays in spectral space -struct FFTData -{ - static constexpr int N = 22; - void* data[N] = { nullptr }; - - ~FFTData () { - for (int i = 0; i < N; ++i) { // The memory is allocated with fftw_alloc. - fftw_free(data[i]); - data[i] = nullptr; - } - } - - FFTData () = default; - - FFTData (FFTData && rhs) noexcept { - for (int i = 0; i < N; ++i) { - data[i] = rhs.data[i]; - rhs.data[i] = nullptr; - } - } - - FFTData (FFTData const&) = delete; - void operator= (FFTData const&) = delete; - void operator= (FFTData&&) = delete; -}; - -#endif // WARPX_USE_PSATD -#endif // WARPX_PICSAR_HYBRID_FFTDATA_H_ diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp deleted file mode 100644 index db6942229..000000000 --- a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp +++ /dev/null @@ -1,459 +0,0 @@ -/* Copyright 2019 Axel Huebl, Maxence Thevenet, Remi Lehe - * Revathi Jambunathan, Weiqun Zhang - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#include "PicsarHybridFFTData.H" -#include "WarpX.H" -#include "FortranInterface/WarpX_f.H" -#include <AMReX_iMultiFab.H> - - -#ifdef WARPX_USE_PSATD - -using namespace amrex; - -constexpr int FFTData::N; - -namespace { -static std::unique_ptr<FFTData> nullfftdata; // This for process with nz_fft=0 - -/** \brief Returns an "owner mask" which 1 for all cells, except - * for the duplicated (physical) cells of a nodal grid. - * - * More precisely, for these cells (which are represented on several grids) - * the owner mask is 1 only if these cells are at the lower left end of - * the local grid - or if these cells are at the end of the physical domain - * Therefore, there for these cells, there will be only one grid for - * which the owner mask is non-zero. - */ -static iMultiFab -BuildFFTOwnerMask (const MultiFab& mf, const Geometry& geom) -{ - const BoxArray& ba = mf.boxArray(); - const DistributionMapping& dm = mf.DistributionMap(); - iMultiFab mask(ba, dm, 1, 0); - const int owner = 1; - const int nonowner = 0; - mask.setVal(owner); - - const Box& domain_box = amrex::convert(geom.Domain(), ba.ixType()); - - AMREX_ASSERT(ba.complementIn(domain_box).isEmpty()); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(mask); mfi.isValid(); ++mfi) - { - IArrayBox& fab = mask[mfi]; - const Box& bx = fab.box(); - Box bx2 = bx; - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - // Detect nodal dimensions - if (bx2.type(idim) == IndexType::NODE) { - // Make sure that this grid does not touch the end of - // the physical domain. - if (bx2.bigEnd(idim) < domain_box.bigEnd(idim)) { - bx2.growHi(idim, -1); - } - } - } - const BoxList& bl = amrex::boxDiff(bx, bx2); - // Set owner mask in these cells - for (const auto& b : bl) { - fab.setVal(nonowner, b, 0, 1); - } - - } - - return mask; -} - -/** \brief Copy the data from the FFT grid to the regular grid - * - * Because, for nodal grid, some cells are duplicated on several boxes, - * special care has to be taken in order to have consistent values on - * each boxes when copying this data. Here this is done by setting a - * mask, where, for these duplicated cells, the mask is non-zero on only - * one box. - */ -static void -CopyDataFromFFTToValid (MultiFab& mf, const MultiFab& mf_fft, const BoxArray& ba_valid_fft, const Geometry& geom) -{ - auto idx_type = mf_fft.ixType(); - MultiFab mftmp(amrex::convert(ba_valid_fft,idx_type), mf_fft.DistributionMap(), 1, 0); - - const iMultiFab& mask = BuildFFTOwnerMask(mftmp, geom); - - // Local copy: whenever an MPI rank owns both the data from the FFT - // grid and from the regular grid, for overlapping region, copy it locally -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(mftmp,true); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - FArrayBox& dstfab = mftmp[mfi]; - - const FArrayBox& srcfab = mf_fft[mfi]; - const Box& srcbox = srcfab.box(); - - if (srcbox.contains(bx)) - { - // Copy the interior region (without guard cells) - dstfab.copy(srcfab, bx, 0, bx, 0, 1); - // Set the value to 0 whenever the mask is 0 - // (i.e. for nodal duplicated cells, there is a single box - // for which the mask is different than 0) - // if mask == 0, set value to zero - dstfab.setValIfNot(0.0, bx, mask[mfi], 0, 1); - } - } - - // Global copy: Get the remaining the data from other procs - // Use ParallelAdd instead of ParallelCopy, so that the value from - // the cell that has non-zero mask is the one which is retained. - mf.setVal(0.0, 0); - mf.ParallelAdd(mftmp); - - -} - -} - -void -WarpX::AllocLevelDataFFT (int lev) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lev == 0, "PSATD doesn't work with mesh refinement yet"); - - static_assert(std::is_standard_layout<FFTData>::value, "FFTData must have standard layout"); - static_assert(sizeof(FFTData) == sizeof(void*)*FFTData::N, "sizeof FFTData is wrong"); - - - - InitFFTComm(lev); - - BoxArray ba_fp_fft; - DistributionMapping dm_fp_fft; - FFTDomainDecomposition(lev, ba_fp_fft, dm_fp_fft, ba_valid_fp_fft[lev], domain_fp_fft[lev], - geom[lev].Domain()); - - // rho2 has one extra ghost cell, so that it's safe to deposit charge density after - // pushing particle. - - Efield_fp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_fp_fft,Ex_nodal_flag), - dm_fp_fft, 1, 0)); - Efield_fp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_fp_fft,Ey_nodal_flag), - dm_fp_fft, 1, 0)); - Efield_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,Ez_nodal_flag), - dm_fp_fft, 1, 0)); - Bfield_fp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_fp_fft,Bx_nodal_flag), - dm_fp_fft, 1, 0)); - Bfield_fp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_fp_fft,By_nodal_flag), - dm_fp_fft, 1, 0)); - Bfield_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,Bz_nodal_flag), - dm_fp_fft, 1, 0)); - current_fp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_fp_fft,jx_nodal_flag), - dm_fp_fft, 1, 0)); - current_fp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_fp_fft,jy_nodal_flag), - dm_fp_fft, 1, 0)); - current_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,jz_nodal_flag), - dm_fp_fft, 1, 0)); - rho_fp_fft[lev].reset(new MultiFab(amrex::convert(ba_fp_fft,IntVect::TheNodeVector()), - dm_fp_fft, 2, 0)); - - dataptr_fp_fft[lev].reset(new LayoutData<FFTData>(ba_fp_fft, dm_fp_fft)); - - if (lev > 0) - { - BoxArray ba_cp_fft; - DistributionMapping dm_cp_fft; - FFTDomainDecomposition(lev, ba_cp_fft, dm_cp_fft, ba_valid_cp_fft[lev], domain_cp_fft[lev], - amrex::coarsen(geom[lev].Domain(),2)); - - Efield_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,Ex_nodal_flag), - dm_cp_fft, 1, 0)); - Efield_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,Ey_nodal_flag), - dm_cp_fft, 1, 0)); - Efield_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,Ez_nodal_flag), - dm_cp_fft, 1, 0)); - Bfield_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,Bx_nodal_flag), - dm_cp_fft, 1, 0)); - Bfield_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,By_nodal_flag), - dm_cp_fft, 1, 0)); - Bfield_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,Bz_nodal_flag), - dm_cp_fft, 1, 0)); - current_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,jx_nodal_flag), - dm_cp_fft, 1, 0)); - current_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,jy_nodal_flag), - dm_cp_fft, 1, 0)); - current_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,jz_nodal_flag), - dm_cp_fft, 1, 0)); - rho_cp_fft[lev].reset(new MultiFab(amrex::convert(ba_cp_fft,IntVect::TheNodeVector()), - dm_cp_fft, 2, 0)); - - dataptr_cp_fft[lev].reset(new LayoutData<FFTData>(ba_cp_fft, dm_cp_fft)); - } - - InitFFTDataPlan(lev); -} - -/** \brief Create MPI sub-communicators for each FFT group, - * and put them in PICSAR module - * - * These communicators are passed to the parallel FFTW library, in order - * to perform a global FFT within each FFT group. - */ -void -WarpX::InitFFTComm (int lev) -{ - int nprocs = ParallelDescriptor::NProcs(); - ngroups_fft = std::min(ngroups_fft, nprocs); - - // # of processes in the subcommunicator - int np_fft = nprocs / ngroups_fft; - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(np_fft*ngroups_fft == nprocs, - "Number of processes must be divisible by number of FFT groups"); - - int myproc = ParallelDescriptor::MyProc(); - // my color in ngroups_fft subcommunicators. 0 <= color_fft < ngroups_fft - color_fft[lev] = myproc / np_fft; - MPI_Comm_split(ParallelDescriptor::Communicator(), color_fft[lev], myproc, &comm_fft[lev]); - - int fcomm = MPI_Comm_c2f(comm_fft[lev]); - // Set the communicator of the PICSAR module to the one we just created - warpx_fft_mpi_init(fcomm); -} - -/** \brief Perform domain decomposition for the FFTW - * - * Attribute one (unique) box to each proc, in such a way that: - * - The global domain is divided among FFT groups, - * with additional guard cells around each FFT group - * - The domain associated to an FFT group (with its guard cells) - * is further divided in sub-subdomains along z, so as to distribute - * it among the procs within an FFT group - * - * The attribution is done by setting (within this function): - * - ba_fft: the BoxArray representing the final set of sub-domains for the FFT - * (includes/covers the guard cells of the FFT groups) - * - dm_fft: the mapping between these sub-domains and the corresponding proc - * (imposes one unique box for each proc) - * - ba_valid: the BoxArray that contains valid part of the sub-domains of ba_fft - * (i.e. does not include/cover the guard cells of the FFT groups) - * - domain_fft: a Box that represent the domain of the FFT group for the current proc - */ -void -WarpX::FFTDomainDecomposition (int lev, BoxArray& ba_fft, DistributionMapping& dm_fft, - BoxArray& ba_valid, Box& domain_fft, const Box& domain) -{ - - IntVect nguards_fft(AMREX_D_DECL(nox_fft/2,noy_fft/2,noz_fft/2)); - - int nprocs = ParallelDescriptor::NProcs(); - - BoxList bl(domain, ngroups_fft); // This does a multi-D domain decomposition for groups - AMREX_ALWAYS_ASSERT(bl.size() == ngroups_fft); - const Vector<Box>& bldata = bl.data(); - - // This is the domain for the FFT sub-group (including guard cells) - domain_fft = amrex::grow(bldata[color_fft[lev]], nguards_fft); - // Ask FFTW to chop the current FFT sub-group domain in the z-direction - // and give a chunk to each MPI rank in the current sub-group. - int nz_fft, z0_fft; - - warpx_fft_domain_decomp(&nz_fft, &z0_fft, WARPX_TO_FORTRAN_BOX(domain_fft)); - // Each MPI rank adds a box with its chunk of the FFT grid - // (given by the above decomposition) to the list `bx_fft`, - // then list is shared among all MPI ranks via AllGather - Vector<Box> bx_fft; - if (nz_fft > 0) { - Box b = domain_fft; - b.setRange(AMREX_SPACEDIM-1, z0_fft+domain_fft.smallEnd(AMREX_SPACEDIM-1), nz_fft); - bx_fft.push_back(b); - } else { - // Add empty box for the AllGather call - bx_fft.push_back(Box()); - } - amrex::AllGatherBoxes(bx_fft); - AMREX_ASSERT(bx_fft.size() == ParallelDescriptor::NProcs()); - // Build pmap and bx_fft without the empty boxes - Vector<int> pmap; - for (int i = 0; i < bx_fft.size(); ++i) { - if (bx_fft[i].ok()) { - pmap.push_back(i); - } - } - bx_fft.erase(std::remove_if(bx_fft.begin(),bx_fft.end(), - [](Box const& b) { return b.isEmpty(); }), - bx_fft.end()); - AMREX_ASSERT(bx_fft.size() == pmap.size()); - - // Define the AMReX objects for the FFT grid: BoxArray and DistributionMapping - ba_fft.define(BoxList(std::move(bx_fft))); - dm_fft.define(std::move(pmap)); - - // For communication between WarpX normal domain and FFT domain, we need to create a - // special BoxArray ba_valid - const Box foobox(-nguards_fft-2, -nguards_fft-2); - - BoxList bl_valid; // List of boxes: will be filled by the valid part of the subdomains of ba_fft - bl_valid.reserve(ba_fft.size()); - int np_fft = nprocs / ngroups_fft; - for (int i = 0; i < ba_fft.size(); ++i) - { - int igroup = dm_fft[i] / np_fft; // This should be consistent with InitFFTComm - const Box& bx = ba_fft[i] & bldata[igroup]; // Intersection with the domain of - // the FFT group *without* guard cells - if (bx.ok()) - { - bl_valid.push_back(bx); - } - else - { - bl_valid.push_back(foobox); - } - } - - ba_valid.define(std::move(bl_valid)); -} - -/** /brief Set all the flags and metadata of the PICSAR FFT module. - * Allocate the auxiliary arrays of `fft_data` - * - * Note: dataptr_data is a stuct containing 22 pointers to arrays - * 1-11: padded arrays in real space ; 12-22 arrays for the fields in Fourier space - */ -void -WarpX::InitFFTDataPlan (int lev) -{ - auto dx_fp = CellSize(lev); - - if (Efield_fp_fft[lev][0]->local_size() == 1) - //Only one FFT patch on this MPI - { - for (MFIter mfi(*Efield_fp_fft[lev][0]); mfi.isValid(); ++mfi) - { - warpx_fft_dataplan_init(&nox_fft, &noy_fft, &noz_fft, - (*dataptr_fp_fft[lev])[mfi].data, &FFTData::N, - dx_fp.data(), &dt[lev], &fftw_plan_measure, &WarpX::do_nodal ); - } - } - else if (Efield_fp_fft[lev][0]->local_size() == 0) - // No FFT patch on this MPI rank (may happen with FFTW) - // Still need to call the MPI-FFT initialization routines - { - nullfftdata.reset(new FFTData()); - warpx_fft_dataplan_init(&nox_fft, &noy_fft, &noz_fft, - nullfftdata->data, &FFTData::N, - dx_fp.data(), &dt[lev], &fftw_plan_measure, - &WarpX::do_nodal ); - } - else - { - // Multiple FFT patches on this MPI rank - amrex::Abort("WarpX::InitFFTDataPlan: TODO"); - } - - if (lev > 0) - { - amrex::Abort("WarpX::InitFFTDataPlan: TODO"); - } -} - -void -WarpX::FreeFFT (int lev) -{ - nullfftdata.reset(); - - warpx_fft_nullify(); - - if (comm_fft[lev] != MPI_COMM_NULL) { - MPI_Comm_free(&comm_fft[lev]); - } - comm_fft[lev] = MPI_COMM_NULL; -} - -void -WarpX::PushPSATD_hybridFFT (int lev, amrex::Real /* dt */) -{ - WARPX_PROFILE_VAR_NS("WarpXFFT::CopyDualGrid", blp_copy); - WARPX_PROFILE_VAR_NS("PICSAR::FftPushEB", blp_push_eb); - - auto period_fp = geom[lev].periodicity(); - - WARPX_PROFILE_VAR_START(blp_copy); - Efield_fp_fft[lev][0]->ParallelCopy(*Efield_fp[lev][0], 0, 0, 1, Efield_fp[lev][0]->nGrow(), 0, period_fp); - Efield_fp_fft[lev][1]->ParallelCopy(*Efield_fp[lev][1], 0, 0, 1, Efield_fp[lev][1]->nGrow(), 0, period_fp); - Efield_fp_fft[lev][2]->ParallelCopy(*Efield_fp[lev][2], 0, 0, 1, Efield_fp[lev][2]->nGrow(), 0, period_fp); - Bfield_fp_fft[lev][0]->ParallelCopy(*Bfield_fp[lev][0], 0, 0, 1, Bfield_fp[lev][0]->nGrow(), 0, period_fp); - Bfield_fp_fft[lev][1]->ParallelCopy(*Bfield_fp[lev][1], 0, 0, 1, Bfield_fp[lev][1]->nGrow(), 0, period_fp); - Bfield_fp_fft[lev][2]->ParallelCopy(*Bfield_fp[lev][2], 0, 0, 1, Bfield_fp[lev][2]->nGrow(), 0, period_fp); - current_fp_fft[lev][0]->ParallelCopy(*current_fp[lev][0], 0, 0, 1, current_fp[lev][0]->nGrow(), 0, period_fp); - current_fp_fft[lev][1]->ParallelCopy(*current_fp[lev][1], 0, 0, 1, current_fp[lev][1]->nGrow(), 0, period_fp); - current_fp_fft[lev][2]->ParallelCopy(*current_fp[lev][2], 0, 0, 1, current_fp[lev][2]->nGrow(), 0, period_fp); - rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 2, rho_fp[lev]->nGrow(), 0, period_fp); - WARPX_PROFILE_VAR_STOP(blp_copy); - - WARPX_PROFILE_VAR_START(blp_push_eb); - if (Efield_fp_fft[lev][0]->local_size() == 1) - //Only one FFT patch on this MPI - { - for (MFIter mfi(*Efield_fp_fft[lev][0]); mfi.isValid(); ++mfi) - { - warpx_fft_push_eb(WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][0])[mfi]), - WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][1])[mfi]), - WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][2])[mfi]), - WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][0])[mfi]), - WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][1])[mfi]), - WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][2])[mfi]), - WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][0])[mfi]), - WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][1])[mfi]), - WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][2])[mfi]), - WARPX_TO_FORTRAN_N_ANYD((*rho_fp_fft[lev])[mfi],0), - WARPX_TO_FORTRAN_N_ANYD((*rho_fp_fft[lev])[mfi],1)); - } - } - else if (Efield_fp_fft[lev][0]->local_size() == 0) - // No FFT patch on this MPI rank - // Still need to call the MPI-FFT routine. - { - FArrayBox fab(Box(IntVect::TheZeroVector(), IntVect::TheUnitVector())); - warpx_fft_push_eb(WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab), - WARPX_TO_FORTRAN_ANYD(fab)); - } - else - // Multiple FFT patches on this MPI rank - { - amrex::Abort("WarpX::PushPSATD: TODO"); - } - WARPX_PROFILE_VAR_STOP(blp_push_eb); - - WARPX_PROFILE_VAR_START(blp_copy); - CopyDataFromFFTToValid(*Efield_fp[lev][0], *Efield_fp_fft[lev][0], ba_valid_fp_fft[lev], geom[lev]); - CopyDataFromFFTToValid(*Efield_fp[lev][1], *Efield_fp_fft[lev][1], ba_valid_fp_fft[lev], geom[lev]); - CopyDataFromFFTToValid(*Efield_fp[lev][2], *Efield_fp_fft[lev][2], ba_valid_fp_fft[lev], geom[lev]); - CopyDataFromFFTToValid(*Bfield_fp[lev][0], *Bfield_fp_fft[lev][0], ba_valid_fp_fft[lev], geom[lev]); - CopyDataFromFFTToValid(*Bfield_fp[lev][1], *Bfield_fp_fft[lev][1], ba_valid_fp_fft[lev], geom[lev]); - CopyDataFromFFTToValid(*Bfield_fp[lev][2], *Bfield_fp_fft[lev][2], ba_valid_fp_fft[lev], geom[lev]); - WARPX_PROFILE_VAR_STOP(blp_copy); - - if (lev > 0) - { - amrex::Abort("WarpX::PushPSATD: TODO"); - } -} - -#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 b/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 deleted file mode 100644 index 96bfb6111..000000000 --- a/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 +++ /dev/null @@ -1,326 +0,0 @@ -! Copyright 2019 Andrew Myers, Maxence Thevenet, Remi Lehe -! Weiqun Zhang -! -! This file is part of WarpX. -! -! License: BSD-3-Clause-LBNL - - -module warpx_fft_module - use amrex_error_module, only : amrex_error, amrex_abort - use amrex_fort_module, only : amrex_real - use iso_c_binding - implicit none - - include 'fftw3-mpi.f03' - - private - public :: warpx_fft_mpi_init, warpx_fft_domain_decomp, warpx_fft_dataplan_init, warpx_fft_nullify, & - warpx_fft_push_eb - -contains - -!> @brief -!! Set the communicator of the PICSAR module to the one that is passed in argument - subroutine warpx_fft_mpi_init (fcomm) bind(c,name='warpx_fft_mpi_init') - use shared_data, only : comm, rank, nproc - integer, intent(in), value :: fcomm - - integer :: ierr, lnproc, lrank - - comm = fcomm - - call mpi_comm_size(comm, lnproc, ierr) - nproc = lnproc - - call mpi_comm_rank(comm, lrank, ierr) - rank = lrank - -#ifdef _OPENMP - ierr = fftw_init_threads() - if (ierr.eq.0) call amrex_error("fftw_init_threads failed") -#endif - call fftw_mpi_init() -#ifdef _OPENMP - call dfftw_init_threads(ierr) - if (ierr.eq.0) call amrex_error("dfftw_init_threads failed") -#endif - end subroutine warpx_fft_mpi_init - -!> @brief -!! Ask FFTW to do domain decomposition. -! -! This is always a 1d domain decomposition along z ; it is typically -! done on the *FFT sub-groups*, not the all domain - subroutine warpx_fft_domain_decomp (warpx_local_nz, warpx_local_z0, global_lo, global_hi) & - bind(c,name='warpx_fft_domain_decomp') - use picsar_precision, only : idp - use shared_data, only : comm, & - nx_global, ny_global, nz_global, & ! size of global FFT - nx, ny, nz ! size of local subdomains - use mpi_fftw3, only : local_nz, local_z0, fftw_mpi_local_size_3d, alloc_local - - integer, intent(out) :: warpx_local_nz, warpx_local_z0 - integer, dimension(3), intent(in) :: global_lo, global_hi - - nx_global = INT(global_hi(1)-global_lo(1)+1,idp) - ny_global = INT(global_hi(2)-global_lo(2)+1,idp) - nz_global = INT(global_hi(3)-global_lo(3)+1,idp) - - alloc_local = fftw_mpi_local_size_3d( & - INT(nz_global,C_INTPTR_T), & - INT(ny_global,C_INTPTR_T), & - INT(nx_global,C_INTPTR_T)/2+1, & - comm, local_nz, local_z0) - - nx = nx_global - ny = ny_global - nz = local_nz - - warpx_local_nz = local_nz - warpx_local_z0 = local_z0 - end subroutine warpx_fft_domain_decomp - - -!> @brief -!! Set all the flags and metadata of the PICSAR FFT module. -!! Allocate the auxiliary arrays of `fft_data` -!! -!! Note: fft_data is a stuct containing 22 pointers to arrays -!! 1-11: padded arrays in real space ; 12-22 arrays for the fields in Fourier space - subroutine warpx_fft_dataplan_init (nox, noy, noz, fft_data, ndata, dx_wrpx, dt_wrpx, fftw_measure, do_nodal) & - bind(c,name='warpx_fft_dataplan_init') - USE picsar_precision, only: idp - use shared_data, only : c_dim, p3dfft_flag, fftw_plan_measure, & - fftw_with_mpi, fftw_threads_ok, fftw_hybrid, fftw_mpi_transpose, & - nx, ny, nz, & ! size of local subdomains - nkx, nky, nkz, & ! size of local ffts - iz_min_r, iz_max_r, iy_min_r, iy_max_r, ix_min_r, ix_max_r, & ! loop bounds - dx, dy, dz - use fields, only : nxguards, nyguards, nzguards, & ! size of guard regions - ex_r, ey_r, ez_r, bx_r, by_r, bz_r, & - jx_r, jy_r, jz_r, rho_r, rhoold_r, & - exf, eyf, ezf, bxf, byf, bzf, & - jxf, jyf, jzf, rhof, rhooldf, & - l_spectral, l_staggered, norderx, nordery, norderz - use mpi_fftw3, only : alloc_local - use omp_lib, only: omp_get_max_threads - USE gpstd_solver, only: init_gpstd - USE fourier_psaotd, only: init_plans_fourier_mpi - use params, only : dt - - integer, intent(in) :: nox, noy, noz, ndata - integer, intent(in) :: fftw_measure - integer, intent(in) :: do_nodal - type(c_ptr), intent(inout) :: fft_data(ndata) - real(c_double), intent(in) :: dx_wrpx(3), dt_wrpx - - integer(idp) :: nopenmp - integer :: nx_padded - integer, dimension(3) :: shp - integer(kind=c_size_t) :: sz - - ! No need to distinguish physical and guard cells for the global FFT; - ! only nx+2*nxguards counts. Thus we declare 0 guard cells for simplicity - nxguards = 0_idp - nyguards = 0_idp - nzguards = 0_idp - - ! For the calculation of the modified [k] vectors - if (do_nodal == 0) then - l_staggered = .TRUE. - else - l_staggered = .FALSE. - endif - norderx = int(nox, idp) - nordery = int(noy, idp) - norderz = int(noz, idp) - ! Define parameters of FFT plans - c_dim = INT(AMREX_SPACEDIM,idp) ! Dimensionality of the simulation (2d/3d) - fftw_with_mpi = .TRUE. ! Activate MPI FFTW - fftw_hybrid = .FALSE. ! FFT per MPI subgroup (instead of global) - fftw_mpi_transpose = .FALSE. ! Do not transpose the data - fftw_plan_measure = (fftw_measure .ne. 0) - p3dfft_flag = .FALSE. - l_spectral = .TRUE. ! Activate spectral Solver, using FFT -#ifdef _OPENMP - fftw_threads_ok = .TRUE. - nopenmp = OMP_GET_MAX_THREADS() -#else - fftw_threads_ok = .FALSE. - nopenmp = 1 -#endif - - ! Allocate padded arrays for MPI FFTW - nx_padded = 2*(nx/2 + 1) - shp = [nx_padded, int(ny), int(nz)] - sz = 2*alloc_local - fft_data(1) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(1), ex_r, shp) - fft_data(2) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(2), ey_r, shp) - fft_data(3) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(3), ez_r, shp) - fft_data(4) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(4), bx_r, shp) - fft_data(5) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(5), by_r, shp) - fft_data(6) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(6), bz_r, shp) - fft_data(7) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(7), jx_r, shp) - fft_data(8) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(8), jy_r, shp) - fft_data(9) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(9), jz_r, shp) - fft_data(10) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(10), rho_r, shp) - fft_data(11) = fftw_alloc_real(sz) - call c_f_pointer(fft_data(11), rhoold_r, shp) - - ! Set array bounds when copying ex to ex_r in PICSAR - ix_min_r = 1; ix_max_r = nx - iy_min_r = 1; iy_max_r = ny - iz_min_r = 1; iz_max_r = nz - ! Allocate Fourier space fields of the same size - nkx = nx/2 + 1 - nky = ny - nkz = nz - shp = [int(nkx), int(nky), int(nkz)] - sz = alloc_local - fft_data(12) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(12), exf, shp) - fft_data(13) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(13), eyf, shp) - fft_data(14) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(14), ezf, shp) - fft_data(15) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(15), bxf, shp) - fft_data(16) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(16), byf, shp) - fft_data(17) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(17), bzf, shp) - fft_data(18) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(18), jxf, shp) - fft_data(19) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(19), jyf, shp) - fft_data(20) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(20), jzf, shp) - fft_data(21) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(21), rhof, shp) - fft_data(22) = fftw_alloc_complex(sz) - call c_f_pointer(fft_data(22), rhooldf, shp) -!$acc enter data create (exf,eyf,ezf,bxf,byf,bzf,jxf,jyf,jzf,rhof,rhooldf) - if (ndata < 22) then - call amrex_abort("size of fft_data is too small") - end if - - dx = dx_wrpx(1) - dy = dx_wrpx(2) - dz = dx_wrpx(3) - dt = dt_wrpx - - ! Initialize the matrix blocks for the PSATD solver - CALL init_gpstd() - ! Initialize the plans for fftw with MPI - CALL init_plans_fourier_mpi(nopenmp) - - end subroutine warpx_fft_dataplan_init - - - subroutine warpx_fft_nullify () bind(c,name='warpx_fft_nullify') - use fields, only: ex_r, ey_r, ez_r, bx_r, by_r, bz_r, & - jx_r, jy_r, jz_r, rho_r, rhoold_r, & - exf, eyf, ezf, bxf, byf, bzf, & - jxf, jyf, jzf, rhof, rhooldf - use mpi_fftw3, only : plan_r2c_mpi, plan_c2r_mpi - nullify(ex_r) - nullify(ey_r) - nullify(ez_r) - nullify(bx_r) - nullify(by_r) - nullify(bz_r) - nullify(jx_r) - nullify(jy_r) - nullify(jz_r) - nullify(rho_r) - nullify(rhoold_r) - nullify(exf) - nullify(eyf) - nullify(ezf) - nullify(bxf) - nullify(byf) - nullify(bzf) - nullify(jxf) - nullify(jyf) - nullify(jzf) - nullify(rhof) - nullify(rhooldf) - call fftw_destroy_plan(plan_r2c_mpi) - call fftw_destroy_plan(plan_c2r_mpi) - call fftw_mpi_cleanup() - end subroutine warpx_fft_nullify - - - subroutine warpx_fft_push_eb ( & - ex_wrpx, exlo, exhi, & - ey_wrpx, eylo, eyhi, & - ez_wrpx, ezlo, ezhi, & - bx_wrpx, bxlo, bxhi, & - by_wrpx, bylo, byhi, & - bz_wrpx, bzlo, bzhi, & - jx_wrpx, jxlo, jxhi, & - jy_wrpx, jylo, jyhi, & - jz_wrpx, jzlo, jzhi, & - rhoold_wrpx, r1lo, r1hi, & - rho_wrpx, r2lo, r2hi) & - bind(c,name='warpx_fft_push_eb') - - use fields, only: ex, ey, ez, bx, by, bz, jx, jy, jz - use shared_data, only: rhoold, rho - use constants, only: num - - integer, dimension(3), intent(in) :: exlo, exhi, eylo, eyhi, ezlo, ezhi, bxlo, bxhi, & - bylo, byhi, bzlo, bzhi, jxlo, jxhi, jylo, jyhi, jzlo, jzhi, r1lo, r1hi, r2lo, r2hi - REAL(num), INTENT(INOUT), TARGET :: ex_wrpx(0:exhi(1)-exlo(1),0:exhi(2)-exlo(2),0:exhi(3)-exlo(3)) - REAL(num), INTENT(INOUT), TARGET :: ey_wrpx(0:eyhi(1)-eylo(1),0:eyhi(2)-eylo(2),0:eyhi(3)-eylo(3)) - REAL(num), INTENT(INOUT), TARGET :: ez_wrpx(0:ezhi(1)-ezlo(1),0:ezhi(2)-ezlo(2),0:ezhi(3)-ezlo(3)) - REAL(num), INTENT(INOUT), TARGET :: bx_wrpx(0:bxhi(1)-bxlo(1),0:bxhi(2)-bxlo(2),0:bxhi(3)-bxlo(3)) - REAL(num), INTENT(INOUT), TARGET :: by_wrpx(0:byhi(1)-bylo(1),0:byhi(2)-bylo(2),0:byhi(3)-bylo(3)) - REAL(num), INTENT(INOUT), TARGET :: bz_wrpx(0:bzhi(1)-bzlo(1),0:bzhi(2)-bzlo(2),0:bzhi(3)-bzlo(3)) - REAL(num), INTENT(INOUT), TARGET :: jx_wrpx(0:jxhi(1)-jxlo(1),0:jxhi(2)-jxlo(2),0:jxhi(3)-jxlo(3)) - REAL(num), INTENT(INOUT), TARGET :: jy_wrpx(0:jyhi(1)-jylo(1),0:jyhi(2)-jylo(2),0:jyhi(3)-jylo(3)) - REAL(num), INTENT(INOUT), TARGET :: jz_wrpx(0:jzhi(1)-jzlo(1),0:jzhi(2)-jzlo(2),0:jzhi(3)-jzlo(3)) - REAL(num), INTENT(INOUT), TARGET :: rhoold_wrpx(0:r1hi(1)-r1lo(1),0:r1hi(2)-r1lo(2),0:r1hi(3)-r1lo(3)) - REAL(num), INTENT(INOUT), TARGET :: rho_wrpx(0:r2hi(1)-r2lo(1),0:r2hi(2)-r2lo(2),0:r2hi(3)-r2lo(3)) - - ! Point the fields in the PICSAR modules to the fields provided by WarpX - ex => ex_wrpx - ey => ey_wrpx - ez => ez_wrpx - bx => bx_wrpx - by => by_wrpx - bz => bz_wrpx - jx => jx_wrpx - jy => jy_wrpx - jz => jz_wrpx - rho => rho_wrpx - rhoold => rhoold_wrpx - - ! Call the corresponding PICSAR function - CALL push_psatd_ebfield() - - ex => null() - ey => null() - ez => null() - bx => null() - by => null() - bz => null() - jx => null() - jy => null() - jz => null() - rho => null() - rhoold => null() - end subroutine warpx_fft_push_eb - -end module warpx_fft_module diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index c1e0c59d0..c53070827 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -37,7 +37,6 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // For local FFTs, boxes in spectral space start at 0 in // each direction and have the same number of points as the // (cell-centered) real space box - // TODO: this will be different for the hybrid FFT scheme Box realspace_bx = realspace_ba[i]; IntVect fft_size = realspace_bx.length(); // Because the spectral solver uses real-to-complex FFTs, we only @@ -115,7 +114,6 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, k[i] = (i-N)*dk; } } - // TODO: this will be different for the hybrid FFT scheme } return k_comp; } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6bf412b3c..09b1b7dbb 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -95,13 +95,7 @@ WarpX::PushPSATD (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent"); - if (fft_hybrid_mpi_decomposition){ -#ifdef WARPX_USE_PSATD_HYBRID - PushPSATD_hybridFFT(lev, a_dt); -#endif - } else { - PushPSATD_localFFT(lev, a_dt); - } + PushPSATD(lev, a_dt); // Evolve the fields in the PML boxes if (do_pml && pml[lev]->ok()) { @@ -111,7 +105,7 @@ WarpX::PushPSATD (amrex::Real a_dt) } void -WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) +WarpX::PushPSATD (int lev, amrex::Real /* dt */) { // Update the fields on the fine and coarse patch PushPSATDSinglePatch( *spectral_solver_fp[lev], diff --git a/Source/FortranInterface/Make.package b/Source/FortranInterface/Make.package deleted file mode 100644 index a1301fb7d..000000000 --- a/Source/FortranInterface/Make.package +++ /dev/null @@ -1 +0,0 @@ -VPATH_LOCATIONS += $(WARPX_HOME)/Source/FortranInterface diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H deleted file mode 100644 index a264780a3..000000000 --- a/Source/FortranInterface/WarpX_f.H +++ /dev/null @@ -1,62 +0,0 @@ -/* Copyright 2019 Andrew Myers, Aurore Blelly, Axel Huebl - * David Grote, Jean-Luc Vay, Ligia Diana Amorim - * Luca Fedeli, Mathieu Lobet, Maxence Thevenet - * Remi Lehe, Revathi Jambunathan, Weiqun Zhang - * - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef WARPX_F_H_ -#define WARPX_F_H_ - -#include <AMReX_BLFort.H> - -#ifdef __cplusplus - -#if AMREX_SPACEDIM==2 -#define WARPX_ARLIM_ANYD(x) std::array<int,3>{(x)[0], 0, (x)[1]}.data() -#else -#define WARPX_ARLIM_ANYD(x) x -#endif - -#define WARPX_TO_FORTRAN_BOX(x) WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect()) -#define WARPX_TO_FORTRAN_ANYD(x) (x).dataPtr(), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect()) -#define WARPX_TO_FORTRAN_N_ANYD(x,n) (x).dataPtr(n), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect()) - -#endif - -#ifdef __cplusplus -extern "C" -{ -#endif - -#ifdef WARPX_USE_PSATD - void warpx_fft_mpi_init (int fcomm); - void warpx_fft_domain_decomp (int* warpx_local_nz, int* warpx_local_z0, - const int* global_lo, const int* global_hi); - void warpx_fft_dataplan_init (const int* nox, const int* noy, const int* noz, - void* fft_data, const int* ndata, - const amrex_real* dx_w, const amrex_real* dt_w, - const int* fftw_plan_measure, const int* do_nodal ); - void warpx_fft_nullify (); - void warpx_fft_push_eb (amrex_real* ex_w, const int* exlo, const int* exhi, - amrex_real* ey_w, const int* eylo, const int* eyhi, - amrex_real* ez_w, const int* ezlo, const int* ezhi, - amrex_real* bx_w, const int* bxlo, const int* bxhi, - amrex_real* by_w, const int* bylo, const int* byhi, - amrex_real* bz_w, const int* bzlo, const int* bzhi, - amrex_real* jx_w, const int* jxlo, const int* jxhi, - amrex_real* jy_w, const int* jylo, const int* jyhi, - amrex_real* jz_w, const int* jzlo, const int* jzhi, - amrex_real* rhoold_w, const int* r1lo, const int* r1hi, - amrex_real* rho_w, const int* r2lo, const int* r2hi); - -#endif - -#ifdef __cplusplus -} -#endif - -#endif //WARPX_F_H_ diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 1f242f761..8f6dd9004 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -422,41 +422,6 @@ WarpX::InitLevelData (int lev, Real /*time*/) } } -#ifdef WARPX_USE_PSATD_HYBRID - -void -WarpX::InitLevelDataFFT (int lev, Real time) -{ - - Efield_fp_fft[lev][0]->setVal(0.0); - Efield_fp_fft[lev][1]->setVal(0.0); - Efield_fp_fft[lev][2]->setVal(0.0); - Bfield_fp_fft[lev][0]->setVal(0.0); - Bfield_fp_fft[lev][1]->setVal(0.0); - Bfield_fp_fft[lev][2]->setVal(0.0); - current_fp_fft[lev][0]->setVal(0.0); - current_fp_fft[lev][1]->setVal(0.0); - current_fp_fft[lev][2]->setVal(0.0); - rho_fp_fft[lev]->setVal(0.0); - - if (lev > 0) - { - Efield_cp_fft[lev][0]->setVal(0.0); - Efield_cp_fft[lev][1]->setVal(0.0); - Efield_cp_fft[lev][2]->setVal(0.0); - Bfield_cp_fft[lev][0]->setVal(0.0); - Bfield_cp_fft[lev][1]->setVal(0.0); - Bfield_cp_fft[lev][2]->setVal(0.0); - current_cp_fft[lev][0]->setVal(0.0); - current_cp_fft[lev][1]->setVal(0.0); - current_cp_fft[lev][2]->setVal(0.0); - rho_cp_fft[lev]->setVal(0.0); - } - -} - -#endif - void WarpX::InitializeExternalFieldsOnGridUsingParser ( MultiFab *mfx, MultiFab *mfy, MultiFab *mfz, diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 72dce3480..226dc7500 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -63,7 +63,6 @@ include $(WARPX_HOME)/Source/BoundaryConditions/Make.package include $(WARPX_HOME)/Source/Diagnostics/Make.package include $(WARPX_HOME)/Source/FieldSolver/Make.package include $(WARPX_HOME)/Source/Filter/Make.package -include $(WARPX_HOME)/Source/FortranInterface/Make.package include $(WARPX_HOME)/Source/Initialization/Make.package include $(WARPX_HOME)/Source/Laser/Make.package include $(WARPX_HOME)/Source/Parallelization/Make.package @@ -166,11 +165,6 @@ ifeq ($(USE_PSATD),TRUE) INCLUDE_LOCATIONS += $(FFTW_HOME)/include LIBRARY_LOCATIONS += $(FFTW_HOME)/lib endif - ifeq ($(USE_PSATD_PICSAR),TRUE) - # Compile with PICSAR's Hybrid-MPI PSATD - DEFINES += -DWARPX_USE_PSATD_HYBRID - DEFINES += -DFFTW # PICSAR uses it - endif else # Use cuFFT libraries += -lcufft diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H index 70948baa8..c2b9227d2 100644 --- a/Source/Parallelization/GuardCellManager.H +++ b/Source/Parallelization/GuardCellManager.H @@ -24,7 +24,6 @@ public: * \param do_fdtd_nci_corr bool, whether to use Godfrey NCI corrector * \param do_nodal bool, whether the field solver is nodal * \param do_moving_window bool, whether to use moving window - * \param do_fft_mpi_dec bool, whether to do parallel FFTs for PSATD * \param aux_is_nodal bool, true if the aux grid is nodal * \param moving_window_dir direction of moving window * \param nox order of current deposition @@ -40,7 +39,6 @@ public: const bool do_fdtd_nci_corr, const bool do_nodal, const bool do_moving_window, - const bool do_fft_mpi_dec, const bool aux_is_nodal, const int moving_window_dir, const int nox, diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp index f36c64214..c75b7426f 100644 --- a/Source/Parallelization/GuardCellManager.cpp +++ b/Source/Parallelization/GuardCellManager.cpp @@ -17,7 +17,6 @@ guardCellManager::Init( const bool do_fdtd_nci_corr, const bool do_nodal, const bool do_moving_window, - const bool do_fft_mpi_dec, const bool aux_is_nodal, const int moving_window_dir, const int nox, @@ -95,40 +94,38 @@ guardCellManager::Init( ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int)); #ifdef WARPX_USE_PSATD - if (do_fft_mpi_dec == false){ - // All boxes should have the same number of guard cells - // (to avoid temporary parallel copies) - // Thus take the max of the required number of guards for each field - // Also: the number of guard cell should be enough to contain - // the stencil of the FFT solver. Here, this number (`ngFFT`) - // is determined *empirically* to be the order of the solver - // for nodal, and half the order of the solver for staggered. - - int ngFFt_x = do_nodal ? nox_fft : nox_fft/2.; - int ngFFt_y = do_nodal ? noy_fft : noy_fft/2.; - int ngFFt_z = do_nodal ? noz_fft : noz_fft/2.; - - ParmParse pp("psatd"); - pp.query("nx_guard", ngFFt_x); - pp.query("ny_guard", ngFFt_y); - pp.query("nz_guard", ngFFt_z); - - IntVect ngFFT = IntVect(AMREX_D_DECL(ngFFt_x, ngFFt_y, ngFFt_z)); - - for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){ - int ng_required = ngFFT[i_dim]; - // Get the max - ng_required = std::max( ng_required, ng_alloc_EB[i_dim] ); - ng_required = std::max( ng_required, ng_alloc_J[i_dim] ); - ng_required = std::max( ng_required, ng_alloc_Rho[i_dim] ); - ng_required = std::max( ng_required, ng_alloc_F[i_dim] ); - // Set the guard cells to this max - ng_alloc_EB[i_dim] = ng_required; - ng_alloc_J[i_dim] = ng_required; - ng_alloc_F[i_dim] = ng_required; - ng_alloc_Rho[i_dim] = ng_required; - ng_alloc_F_int = ng_required; - } + // All boxes should have the same number of guard cells + // (to avoid temporary parallel copies) + // Thus take the max of the required number of guards for each field + // Also: the number of guard cell should be enough to contain + // the stencil of the FFT solver. Here, this number (`ngFFT`) + // is determined *empirically* to be the order of the solver + // for nodal, and half the order of the solver for staggered. + + int ngFFt_x = do_nodal ? nox_fft : nox_fft/2.; + int ngFFt_y = do_nodal ? noy_fft : noy_fft/2.; + int ngFFt_z = do_nodal ? noz_fft : noz_fft/2.; + + ParmParse pp("psatd"); + pp.query("nx_guard", ngFFt_x); + pp.query("ny_guard", ngFFt_y); + pp.query("nz_guard", ngFFt_z); + + IntVect ngFFT = IntVect(AMREX_D_DECL(ngFFt_x, ngFFt_y, ngFFt_z)); + + for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){ + int ng_required = ngFFT[i_dim]; + // Get the max + ng_required = std::max( ng_required, ng_alloc_EB[i_dim] ); + ng_required = std::max( ng_required, ng_alloc_J[i_dim] ); + ng_required = std::max( ng_required, ng_alloc_Rho[i_dim] ); + ng_required = std::max( ng_required, ng_alloc_F[i_dim] ); + // Set the guard cells to this max + ng_alloc_EB[i_dim] = ng_required; + ng_alloc_J[i_dim] = ng_required; + ng_alloc_F[i_dim] = ng_required; + ng_alloc_Rho[i_dim] = ng_required; + ng_alloc_F_int = ng_required; } ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int)); #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c63b94a27..a085c440d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -11,7 +11,6 @@ #include "PhysicalParticleContainer.H" #include "MultiParticleContainer.H" -#include "FortranInterface/WarpX_f.H" #include "WarpX.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXUtil.H" diff --git a/Source/WarpX.H b/Source/WarpX.H index d343528d7..1a344a785 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -31,9 +31,6 @@ # include "FieldSolver/SpectralSolver/SpectralSolver.H" # endif #endif -#ifdef WARPX_USE_PSATD_HYBRID -# include "FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H" -#endif #include "Parallelization/GuardCellManager.H" @@ -809,21 +806,6 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice; -#ifdef WARPX_USE_PSATD_HYBRID - // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) - // This includes data for the FFT guard cells (between FFT groups) - amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp_fft; - amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_fp_fft; - amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_fft; - amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_fp_fft; - - amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_cp_fft; - amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_cp_fft; - amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_cp_fft; - amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_cp_fft; -#endif - - bool fft_hybrid_mpi_decomposition = false; bool fft_periodic_single_box = false; int nox_fft = 16; int noy_fft = 16; @@ -833,9 +815,8 @@ private: private: void EvolvePSATD (int numsteps); void PushPSATD (amrex::Real dt); - void PushPSATD_localFFT (int lev, amrex::Real dt); + void PushPSATD (int lev, amrex::Real dt); - int ngroups_fft = 4; int fftw_plan_measure = 1; # ifdef WARPX_DIM_RZ @@ -849,33 +830,6 @@ private: amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> m_fdtd_solver_fp; amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> m_fdtd_solver_cp; -#ifdef WARPX_USE_PSATD_HYBRID -private: - amrex::Vector<std::unique_ptr<amrex::LayoutData<FFTData> > > dataptr_fp_fft; - amrex::Vector<std::unique_ptr<amrex::LayoutData<FFTData> > > dataptr_cp_fft; - - // Domain decomposition containing the valid part of the dual grids (i.e. without FFT guard cells) - amrex::Vector<amrex::BoxArray> ba_valid_fp_fft; - amrex::Vector<amrex::BoxArray> ba_valid_cp_fft; - - amrex::Vector<amrex::Box> domain_fp_fft; // "global" domain for the group this process belongs to - amrex::Vector<amrex::Box> domain_cp_fft; - - amrex::Vector<MPI_Comm> comm_fft; - amrex::Vector<int> color_fft; - - void AllocLevelDataFFT (int lev); - void InitLevelDataFFT (int lev, amrex::Real time); - - void InitFFTComm (int lev); - void FFTDomainDecomposition (int lev, amrex::BoxArray& ba_fft, amrex::DistributionMapping& dm_fft, - amrex::BoxArray& ba_valid, amrex::Box& domain_fft, - const amrex::Box& domain); - void InitFFTDataPlan (int lev); - void FreeFFT (int lev); - void PushPSATD_hybridFFT (int lev, amrex::Real dt); -#endif - #if defined(BL_USE_SENSEI_INSITU) amrex::AmrMeshInSituBridge *insitu_bridge; #endif diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 54deb59d4..1130a060f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -303,29 +303,6 @@ WarpX::WarpX () #endif m_fdtd_solver_fp.resize(nlevs_max); m_fdtd_solver_cp.resize(nlevs_max); -#ifdef WARPX_USE_PSATD_HYBRID - Efield_fp_fft.resize(nlevs_max); - Bfield_fp_fft.resize(nlevs_max); - current_fp_fft.resize(nlevs_max); - rho_fp_fft.resize(nlevs_max); - - Efield_cp_fft.resize(nlevs_max); - Bfield_cp_fft.resize(nlevs_max); - current_cp_fft.resize(nlevs_max); - rho_cp_fft.resize(nlevs_max); - - dataptr_fp_fft.resize(nlevs_max); - dataptr_cp_fft.resize(nlevs_max); - - ba_valid_fp_fft.resize(nlevs_max); - ba_valid_cp_fft.resize(nlevs_max); - - domain_fp_fft.resize(nlevs_max); - domain_cp_fft.resize(nlevs_max); - - comm_fft.resize(nlevs_max,MPI_COMM_NULL); - color_fft.resize(nlevs_max,-1); -#endif #ifdef BL_USE_SENSEI_INSITU insitu_bridge = nullptr; @@ -738,9 +715,7 @@ WarpX::ReadParameters () #ifdef WARPX_USE_PSATD { ParmParse pp("psatd"); - pp.query("hybrid_mpi_decomposition", fft_hybrid_mpi_decomposition); pp.query("periodic_single_box_fft", fft_periodic_single_box); - pp.query("ngroups_fft", ngroups_fft); pp.query("fftw_plan_measure", fftw_plan_measure); pp.query("nox", nox_fft); pp.query("noy", noy_fft); @@ -809,17 +784,6 @@ WarpX::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_grids, { AllocLevelData(lev, new_grids, new_dmap); InitLevelData(lev, time); - -#ifdef WARPX_USE_PSATD - if (fft_hybrid_mpi_decomposition){ -#ifdef WARPX_USE_PSATD_HYBRID - AllocLevelDataFFT(lev); - InitLevelDataFFT(lev, time); -#else - amrex::Abort("The option `psatd.fft_hybrid_mpi_decomposition` does not work on GPU."); -#endif - } -#endif } void @@ -855,30 +819,6 @@ WarpX::ClearLevel (int lev) rho_cp[lev].reset(); costs[lev].reset(); - - -#ifdef WARPX_USE_PSATD_HYBRID - for (int i = 0; i < 3; ++i) { - Efield_fp_fft[lev][i].reset(); - Bfield_fp_fft[lev][i].reset(); - current_fp_fft[lev][i].reset(); - - Efield_cp_fft[lev][i].reset(); - Bfield_cp_fft[lev][i].reset(); - current_cp_fft[lev][i].reset(); - } - - rho_fp_fft[lev].reset(); - rho_cp_fft[lev].reset(); - - dataptr_fp_fft[lev].reset(); - dataptr_cp_fft[lev].reset(); - - ba_valid_fp_fft[lev] = BoxArray(); - ba_valid_cp_fft[lev] = BoxArray(); - - FreeFFT(lev); -#endif } void @@ -892,7 +832,6 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d WarpX::use_fdtd_nci_corr, do_nodal, do_moving_window, - fft_hybrid_mpi_decomposition, aux_is_nodal, moving_window_dir, WarpX::nox, @@ -942,6 +881,8 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // // The fine patch // + std::array<Real,3> dx = CellSize(lev); + Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE+ngextra)); Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE+ngextra)); Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE+ngextra)); @@ -975,40 +916,35 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm { rho_fp[lev].reset(new MultiFab(amrex::convert(ba,rho_nodal_flag),dm,2*ncomps,ngRho)); } - if (fft_hybrid_mpi_decomposition == false){ - // Allocate and initialize the spectral solver - std::array<Real,3> dx = CellSize(lev); -#if (AMREX_SPACEDIM == 3) - RealVect dx_vect(dx[0], dx[1], dx[2]); -#elif (AMREX_SPACEDIM == 2) - RealVect dx_vect(dx[0], dx[2]); -#endif - - // Check whether the option periodic, single box is valid here - if (fft_periodic_single_box) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( geom[0].isAllPeriodic() && ba.size()==1 && lev==0, - "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box."); - } - // Get the cell-centered box - BoxArray realspace_ba = ba; // Copy box - realspace_ba.enclosedCells(); // Make it cell-centered - // Define spectral solver -#ifdef WARPX_DIM_RZ - realspace_ba.grow(1, ngE[1]); // add guard cells only in z - spectral_solver_fp[lev].reset( new SpectralSolverRZ( realspace_ba, dm, - n_rz_azimuthal_modes, noz_fft, do_nodal, dx_vect, dt[lev], lev ) ); -#else - if ( fft_periodic_single_box == false ) { - realspace_ba.grow(ngE); // add guard cells - } - bool const pml=false; - spectral_solver_fp[lev].reset( new SpectralSolver( realspace_ba, dm, - nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, dx_vect, dt[lev], - pml, fft_periodic_single_box ) ); -#endif + // Allocate and initialize the spectral solver +# if (AMREX_SPACEDIM == 3) + RealVect dx_vect(dx[0], dx[1], dx[2]); +# elif (AMREX_SPACEDIM == 2) + RealVect dx_vect(dx[0], dx[2]); +# endif + // Check whether the option periodic, single box is valid here + if (fft_periodic_single_box) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( geom[0].isAllPeriodic() && ba.size()==1 && lev==0, + "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box."); + } + // Get the cell-centered box + BoxArray realspace_ba = ba; // Copy box + realspace_ba.enclosedCells(); // Make it cell-centered + // Define spectral solver +# ifdef WARPX_DIM_RZ + realspace_ba.grow(1, ngE[1]); // add guard cells only in z + spectral_solver_fp[lev].reset( new SpectralSolverRZ( realspace_ba, dm, + n_rz_azimuthal_modes, noz_fft, do_nodal, dx_vect, dt[lev], lev ) ); +# else + if ( fft_periodic_single_box == false ) { + realspace_ba.grow(ngE); // add guard cells } + bool const pml=false; + spectral_solver_fp[lev].reset( new SpectralSolver( realspace_ba, dm, + nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, dx_vect, dt[lev], + pml, fft_periodic_single_box ) ); +# endif #endif - std::array<Real,3> const dx = CellSize(lev); m_fdtd_solver_fp[lev].reset( new FiniteDifferenceSolver(maxwell_fdtd_solver_id, dx, do_nodal) ); // @@ -1051,6 +987,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm { BoxArray cba = ba; cba.coarsen(refRatio(lev-1)); + std::array<Real,3> cdx = CellSize(lev-1); // Create the MultiFabs for B Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); @@ -1079,33 +1016,28 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm { rho_cp[lev].reset(new MultiFab(amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho)); } - if (fft_hybrid_mpi_decomposition == false){ - // Allocate and initialize the spectral solver - std::array<Real,3> cdx = CellSize(lev-1); -#if (AMREX_SPACEDIM == 3) - RealVect cdx_vect(cdx[0], cdx[1], cdx[2]); -#elif (AMREX_SPACEDIM == 2) - RealVect cdx_vect(cdx[0], cdx[2]); -#endif - // Get the cell-centered box, with guard cells - BoxArray realspace_ba = cba;// Copy box - realspace_ba.enclosedCells(); // Make it cell-centered - // Define spectral solver -#ifdef WARPX_DIM_RZ - realspace_ba.grow(1, ngE[1]); // add guard cells only in z - spectral_solver_cp[lev].reset( new SpectralSolverRZ( realspace_ba, dm, - n_rz_azimuthal_modes, noz_fft, do_nodal, cdx_vect, dt[lev], lev ) ); - -#else - realspace_ba.grow(ngE); // add guard cells - spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm, - nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, cdx_vect, dt[lev] ) ); -#endif - } + // Allocate and initialize the spectral solver +# if (AMREX_SPACEDIM == 3) + RealVect cdx_vect(cdx[0], cdx[1], cdx[2]); +# elif (AMREX_SPACEDIM == 2) + RealVect cdx_vect(cdx[0], cdx[2]); +# endif + // Get the cell-centered box, with guard cells + BoxArray realspace_ba = cba;// Copy box + realspace_ba.enclosedCells(); // Make it cell-centered + // Define spectral solver +# ifdef WARPX_DIM_RZ + realspace_ba.grow(1, ngE[1]); // add guard cells only in z + spectral_solver_cp[lev].reset( new SpectralSolverRZ( realspace_ba, dm, + n_rz_azimuthal_modes, noz_fft, do_nodal, cdx_vect, dt[lev], lev ) ); +# else + realspace_ba.grow(ngE); // add guard cells + spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm, + nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, cdx_vect, dt[lev] ) ); +# endif #endif - std::array<Real,3> cdx = CellSize(lev-1); - m_fdtd_solver_cp[lev].reset( - new FiniteDifferenceSolver( maxwell_fdtd_solver_id, cdx, do_nodal ) ); + m_fdtd_solver_cp[lev].reset( + new FiniteDifferenceSolver( maxwell_fdtd_solver_id, cdx, do_nodal ) ); } // |