diff options
Diffstat (limited to 'Source')
24 files changed, 627 insertions, 162 deletions
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 60b0b5fa3..e98561be1 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -1,10 +1,10 @@ - #include <cmath> #include <limits> #include <WarpX.H> #include <WarpXConst.H> #include <WarpX_f.H> +#include <WarpXUtil.H> #ifdef WARPX_USE_PY #include <WarpX_py.H> #endif @@ -38,7 +38,7 @@ WarpX::EvolveEM (int numsteps) { Real walltime_beg_step = amrex::second(); - // Start loop on time steps + // Start loop on time steps amrex::Print() << "\nSTEP " << step+1 << " starts ...\n"; #ifdef WARPX_USE_PY if (warpx_py_beforestep) warpx_py_beforestep(); @@ -53,16 +53,16 @@ WarpX::EvolveEM (int numsteps) if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); - // Reset the costs to 0 - for (int lev = 0; lev <= finest_level; ++lev) { - costs[lev]->setVal(0.0); - } + // Reset the costs to 0 + for (int lev = 0; lev <= finest_level; ++lev) { + costs[lev]->setVal(0.0); + } } for (int lev = 0; lev <= finest_level; ++lev) { - // Perform running average of the costs - // (Giving more importance to most recent costs) - (*costs[lev].get()).mult( (1. - 2./load_balance_int) ); + // Perform running average of the costs + // (Giving more importance to most recent costs) + (*costs[lev].get()).mult( (1. - 2./load_balance_int) ); } } @@ -81,8 +81,8 @@ WarpX::EvolveEM (int numsteps) } is_synchronized = false; } else { - // Beyond one step, we have E^{n} and B^{n}. - // Particles have p^{n-1/2} and x^{n}. + // Beyond one step, we have E^{n} and B^{n}. + // Particles have p^{n-1/2} and x^{n}. FillBoundaryE(); FillBoundaryB(); UpdateAuxilaryData(); @@ -97,6 +97,10 @@ WarpX::EvolveEM (int numsteps) amrex::Abort("Unsupported do_subcycling type"); } + if (num_mirrors>0){ + applyMirrors(cur_time); + } + #ifdef WARPX_USE_PY if (warpx_py_beforeEsolve) warpx_py_beforeEsolve(); #endif @@ -105,8 +109,10 @@ WarpX::EvolveEM (int numsteps) UpdateAuxilaryData(); for (int lev = 0; lev <= finest_level; ++lev) { mypc->PushP(lev, 0.5*dt[lev], - *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], - *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); + *Efield_aux[lev][0],*Efield_aux[lev][1], + *Efield_aux[lev][2], + *Bfield_aux[lev][0],*Bfield_aux[lev][1], + *Bfield_aux[lev][2]); } is_synchronized = true; } @@ -118,18 +124,18 @@ WarpX::EvolveEM (int numsteps) ++istep[lev]; } - cur_time += dt[0]; + cur_time += dt[0]; bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && - (insitu_int > 0) && ((step+1) % insitu_int == 0); + (insitu_int > 0) && ((step+1) % insitu_int == 0); bool move_j = is_synchronized || to_make_plot || do_insitu; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. - int num_moved = MoveWindow(move_j); + int num_moved = MoveWindow(move_j); if (max_level == 0) { int num_redistribute_ghost = num_moved + 1; @@ -139,11 +145,11 @@ WarpX::EvolveEM (int numsteps) mypc->Redistribute(); } - bool to_sort = (sort_int > 0) && ((step+1) % sort_int == 0); - if (to_sort) { - amrex::Print() << "re-sorting particles \n"; - mypc->SortParticlesByCell(); - } + bool to_sort = (sort_int > 0) && ((step+1) % sort_int == 0); + if (to_sort) { + amrex::Print() << "re-sorting particles \n"; + mypc->SortParticlesByCell(); + } amrex::Print()<< "STEP " << step+1 << " ends." << " TIME = " << cur_time << " DT = " << dt[0] << "\n"; @@ -153,10 +159,10 @@ WarpX::EvolveEM (int numsteps) << " s; This step = " << walltime_end_step-walltime_beg_step << " s; Avg. per step = " << walltime/(step+1) << " s\n"; - // sync up time - for (int i = 0; i <= max_level; ++i) { - t_new[i] = cur_time; - } + // sync up time + for (int i = 0; i <= max_level; ++i) { + t_new[i] = cur_time; + } if (do_boosted_frame_diagnostic) { std::unique_ptr<MultiFab> cell_centered_data = nullptr; @@ -166,7 +172,7 @@ WarpX::EvolveEM (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - if (to_make_plot || do_insitu) + if (to_make_plot || do_insitu) { FillBoundaryE(); FillBoundaryB(); @@ -178,34 +184,34 @@ WarpX::EvolveEM (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } - last_plot_file_step = step+1; - last_insitu_step = step+1; + last_plot_file_step = step+1; + last_insitu_step = step+1; - if (to_make_plot) - WritePlotFile(); - - if (do_insitu) - UpdateInSitu(); - } + if (to_make_plot) + WritePlotFile(); + if (do_insitu) + UpdateInSitu(); + } - if (check_int > 0 && (step+1) % check_int == 0) { - last_check_file_step = step+1; - WriteCheckPointFile(); - } + if (check_int > 0 && (step+1) % check_int == 0) { + last_check_file_step = step+1; + WriteCheckPointFile(); + } - if (cur_time >= stop_time - 1.e-3*dt[0]) { - max_time_reached = true; - break; - } + if (cur_time >= stop_time - 1.e-3*dt[0]) { + max_time_reached = true; + break; + } #ifdef WARPX_USE_PY if (warpx_py_afterstep) warpx_py_afterstep(); #endif - // End loop on time steps + // End loop on time steps } - bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step && (max_time_reached || istep[0] >= max_step); + bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step + && (max_time_reached || istep[0] >= max_step); bool do_insitu = (insitu_start >= istep[0]) && (insitu_int > 0) && (istep[0] > last_insitu_step) && (max_time_reached || istep[0] >= max_step); @@ -218,8 +224,10 @@ WarpX::EvolveEM (int numsteps) for (int lev = 0; lev <= finest_level; ++lev) { mypc->FieldGather(lev, - *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], - *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); + *Efield_aux[lev][0],*Efield_aux[lev][1], + *Efield_aux[lev][2], + *Bfield_aux[lev][0],*Bfield_aux[lev][1], + *Bfield_aux[lev][2]); } if (write_plot_file) @@ -229,8 +237,9 @@ WarpX::EvolveEM (int numsteps) UpdateInSitu(); } - if (check_int > 0 && istep[0] > last_check_file_step && (max_time_reached || istep[0] >= max_step)) { - WriteCheckPointFile(); + if (check_int > 0 && istep[0] > last_check_file_step && + (max_time_reached || istep[0] >= max_step)) { + WriteCheckPointFile(); } if (do_boosted_frame_diagnostic) { @@ -462,23 +471,23 @@ WarpX::ComputeDt () Real deltat = 0.; if (maxwell_fdtd_solver_id == 0) { - // CFL time step Yee solver + // CFL time step Yee solver #ifdef WARPX_RZ - // Derived semi-analytically by R. Lehe - deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); + // Derived semi-analytically by R. Lehe + deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); #else - deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), - + 1./(dx[1]*dx[1]), - + 1./(dx[2]*dx[2]))) * PhysConst::c ); + deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), + + 1./(dx[1]*dx[1]), + + 1./(dx[2]*dx[2]))) * PhysConst::c ); #endif } else { - // CFL time step CKC solver + // CFL time step CKC solver #if (BL_SPACEDIM == 3) - const Real delta = std::min(dx[0],std::min(dx[1],dx[2])); + const Real delta = std::min(dx[0],std::min(dx[1],dx[2])); #elif (BL_SPACEDIM == 2) - const Real delta = std::min(dx[0],dx[1]); + const Real delta = std::min(dx[0],dx[1]); #endif - deltat = cfl*delta/PhysConst::c; + deltat = cfl*delta/PhysConst::c; } dt.resize(0); dt.resize(max_level+1,deltat); @@ -493,3 +502,61 @@ WarpX::ComputeDt () dt[0] = const_dt; } } + +/* \brief Apply perfect mirror condition inside the box (not at a boundary). + * In practice, set all fields to 0 on a section of the simulation domain + * (as for a perfect conductor with a given thickness). + * The mirror normal direction has to be parallel to the z axis. + */ +void +WarpX::applyMirrors(Real time){ + // Loop over the mirrors + for(int i_mirror=0; i_mirror<num_mirrors; ++i_mirror){ + // Get mirror properties (lower and upper z bounds) + Real z_min = mirror_z[i_mirror]; + Real z_max_tmp = z_min + mirror_z_width[i_mirror]; + // Boost quantities for boosted frame simulations + if (gamma_boost>1){ + z_min = z_min/gamma_boost - PhysConst::c*beta_boost*time; + z_max_tmp = z_max_tmp/gamma_boost - PhysConst::c*beta_boost*time; + } + // Loop over levels + for(int lev=0; lev<=finest_level; lev++){ + // Make sure that the mirror contains at least + // mirror_z_npoints[i_mirror] cells + Real dz = WarpX::CellSize(lev)[2]; + Real z_max = std::max(z_max_tmp, + z_min+mirror_z_npoints[i_mirror]*dz); + // Get fine patch field MultiFabs + MultiFab& Ex = *Efield_fp[lev][0].get(); + MultiFab& Ey = *Efield_fp[lev][1].get(); + MultiFab& Ez = *Efield_fp[lev][2].get(); + MultiFab& Bx = *Bfield_fp[lev][0].get(); + MultiFab& By = *Bfield_fp[lev][1].get(); + MultiFab& Bz = *Bfield_fp[lev][2].get(); + // Set each field to zero between z_min and z_max + NullifyMF(Ex, lev, z_min, z_max); + NullifyMF(Ey, lev, z_min, z_max); + NullifyMF(Ez, lev, z_min, z_max); + NullifyMF(Bx, lev, z_min, z_max); + NullifyMF(By, lev, z_min, z_max); + NullifyMF(Bz, lev, z_min, z_max); + if (lev>0){ + // Get coarse patch field MultiFabs + MultiFab& Ex = *Efield_cp[lev][0].get(); + MultiFab& Ey = *Efield_cp[lev][1].get(); + MultiFab& Ez = *Efield_cp[lev][2].get(); + MultiFab& Bx = *Bfield_cp[lev][0].get(); + MultiFab& By = *Bfield_cp[lev][1].get(); + MultiFab& Bz = *Bfield_cp[lev][2].get(); + // Set each field to zero between z_min and z_max + NullifyMF(Ex, lev, z_min, z_max); + NullifyMF(Ey, lev, z_min, z_max); + NullifyMF(Ez, lev, z_min, z_max); + NullifyMF(Bx, lev, z_min, z_max); + NullifyMF(By, lev, z_min, z_max); + NullifyMF(Bz, lev, z_min, z_max); + } + } + } +} diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index 50127914d..b0ee54095 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,11 +1,12 @@ CEXE_headers += WarpX_Complex.H CEXE_headers += SpectralSolver.H +CEXE_sources += SpectralSolver.cpp CEXE_headers += SpectralFieldData.H CEXE_sources += SpectralFieldData.cpp -CEXE_headers += PsatdAlgorithm.H -CEXE_sources += PsatdAlgorithm.cpp CEXE_headers += SpectralKSpace.H CEXE_sources += SpectralKSpace.cpp +include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H deleted file mode 100644 index acefcc466..000000000 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef WARPX_PSATD_ALGORITHM_H_ -#define WARPX_PSATD_ALGORITHM_H_ - -#include <SpectralKSpace.H> -#include <SpectralFieldData.H> - -/* \brief Class that updates the field in spectral space - * and stores the coefficients of the corresponding update equation. - */ -class PsatdAlgorithm -{ - using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; - - public: - PsatdAlgorithm(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const bool nodal, const amrex::Real dt); - PsatdAlgorithm() = default; // Default constructor - PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; - void pushSpectralFields(SpectralFieldData& f) const; - - private: - // Modified finite-order vectors - KVectorComponent modified_kx_vec, modified_kz_vec; -#if (AMREX_SPACEDIM==3) - KVectorComponent modified_ky_vec; -#endif - SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; -}; - -#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package new file mode 100644 index 000000000..c62c21f44 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += SpectralBaseAlgorithm.H +CEXE_headers += PsatdAlgorithm.H +CEXE_sources += PsatdAlgorithm.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H new file mode 100644 index 000000000..0487e5226 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -0,0 +1,24 @@ +#ifndef WARPX_PSATD_ALGORITHM_H_ +#define WARPX_PSATD_ALGORITHM_H_ + +#include <SpectralBaseAlgorithm.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithm : public SpectralBaseAlgorithm +{ + public: + PsatdAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::Real dt); + // Redefine update equation from base class + virtual void pushSpectralFields(SpectralFieldData& f) const override final; + + private: + SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; +}; + +#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index ada7506c3..37892d35a 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -9,14 +9,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const Real dt) -// Compute and assign the modified k vectors -: modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)), -#if (AMREX_SPACEDIM==3) - modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)), - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal)) -#else - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) -#endif + // Initialize members of base class + : SpectralBaseAlgorithm( spectral_kspace, dm, + norder_x, norder_y, norder_z, nodal ) { const BoxArray& ba = spectral_kspace.spectralspace_ba; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H new file mode 100644 index 000000000..602eb2473 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -0,0 +1,51 @@ +#ifndef WARPX_SPECTRAL_BASE_ALGORITHM_H_ +#define WARPX_SPECTRAL_BASE_ALGORITHM_H_ + +#include <SpectralKSpace.H> +#include <SpectralFieldData.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + * + * `SpectralBaseAlgorithm` is only a base class and cannot be used directly. + * Instead use its subclasses, which implement the specific field update + * equations for a given spectral algorithm. + */ +class SpectralBaseAlgorithm +{ + public: + // Member function that updates the fields in spectral space ; + // meant to be overridden in subclasses + virtual void pushSpectralFields(SpectralFieldData& f) const = 0; + // The destructor should also be a virtual function, so that + // a pointer to subclass of `SpectraBaseAlgorithm` actually + // calls the subclass's destructor. + virtual ~SpectralBaseAlgorithm() {}; + + protected: // Meant to be used in the subclasses + + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; + + // Constructor + SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal) + // Compute and assign the modified k vectors + : modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)), +#if (AMREX_SPACEDIM==3) + modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)), + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal)) +#else + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) +#endif + {}; + + // Modified finite-order vectors + KVectorComponent modified_kx_vec, modified_kz_vec; +#if (AMREX_SPACEDIM==3) + KVectorComponent modified_ky_vec; +#endif +}; + +#endif // WARPX_SPECTRAL_BASE_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 7444452af..d4019a9a3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -1,8 +1,7 @@ #ifndef WARPX_SPECTRAL_SOLVER_H_ #define WARPX_SPECTRAL_SOLVER_H_ -#include <SpectralKSpace.H> -#include <PsatdAlgorithm.H> +#include <SpectralBaseAlgorithm.H> #include <SpectralFieldData.H> /* \brief Top-level class for the electromagnetic spectral solver @@ -14,28 +13,17 @@ class SpectralSolver { public: - // Inline definition of the member functions of `SpectralSolver` + // Inline definition of the member functions of `SpectralSolver`, + // except the constructor (see `SpectralSolver.cpp`) // The body of these functions is short, since the work is done in the // underlying classes `SpectralFieldData` and `PsatdAlgorithm` - /* \brief Initialize the spectral solver */ + // Constructor SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { - // Initialize all structures using the same distribution mapping dm - - // - Initialize k space object (Contains info about the size of - // the spectral space corresponding to each box in `realspace_ba`, - // as well as the value of the corresponding k coordinates) - const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); - // - Initialize the algorithm (coefficients) over this space - algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, - norder_z, nodal, dt ); - // - Initialize arrays for fields in Fourier space + FFT plans - field_data = SpectralFieldData( realspace_ba, k_space, dm ); - }; + const amrex::RealVect dx, const amrex::Real dt ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally @@ -59,14 +47,20 @@ class SpectralSolver /* \brief Update the fields in spectral space, over one timestep */ void pushSpectralFields(){ BL_PROFILE("SpectralSolver::pushSpectralFields"); - algorithm.pushSpectralFields( field_data ); + // Virtual function: the actual function used here depends + // on the sub-class of `SpectralBaseAlgorithm` that was + // initialized in the constructor of `SpectralSolver` + algorithm->pushSpectralFields( field_data ); }; private: SpectralFieldData field_data; // Store field in spectral space // and perform the Fourier transforms - PsatdAlgorithm algorithm; // Contains Psatd coefficients - // and field update equation + std::unique_ptr<SpectralBaseAlgorithm> algorithm; + // Defines field update equation in spectral space, + // and the associated coefficients. + // SpectralBaseAlgorithm is a base class ; this pointer is meant + // to point an instance of a *sub-class* defining a specific algorithm }; #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp new file mode 100644 index 000000000..c21c3cfb1 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -0,0 +1,35 @@ +#include <SpectralKSpace.H> +#include <SpectralSolver.H> +#include <PsatdAlgorithm.H> + +/* \brief Initialize the spectral Maxwell solver + * + * This function selects the spectral algorithm to be used, allocates the + * corresponding coefficients for the discretized field update equation, + * and prepares the structures that store the fields in spectral space. + */ +SpectralSolver::SpectralSolver( + const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::RealVect dx, const amrex::Real dt ) { + + // Initialize all structures using the same distribution mapping dm + + // - Initialize k space object (Contains info about the size of + // the spectral space corresponding to each box in `realspace_ba`, + // as well as the value of the corresponding k coordinates) + const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); + + // - Select the algorithm depending on the input parameters + // Initialize the corresponding coefficients over k space + // TODO: Add more algorithms + selection depending on input parameters + // For the moment, this only uses the standard PsatdAlgorithm + algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + + // - Initialize arrays for fields in spectral space + FFT plans + field_data = SpectralFieldData( realspace_ba, k_space, dm ); + +}; diff --git a/Source/FortranInterface/WarpX_f.F90 b/Source/FortranInterface/WarpX_f.F90 index ccd6dd48a..d7f7a8053 100644 --- a/Source/FortranInterface/WarpX_f.F90 +++ b/Source/FortranInterface/WarpX_f.F90 @@ -186,6 +186,41 @@ contains end subroutine warpx_compute_dive_2d + subroutine warpx_compute_dive_rz (lo, hi, dive, dlo, dhi, & + Ex, xlo, xhi, Ey, ylo, yhi, Ez, zlo, zhi, dx, rmin) & + bind(c, name='warpx_compute_dive_rz') + integer, intent(in) :: lo(2),hi(2),dlo(2),dhi(2),xlo(2),xhi(2),ylo(2),yhi(2),zlo(2),zhi(2) + real(amrex_real), intent(in) :: dx(3), rmin + real(amrex_real), intent(in ) :: Ex (xlo(1):xhi(1),xlo(2):xhi(2)) + real(amrex_real), intent(in ) :: Ey (ylo(1):yhi(1),ylo(2):yhi(2)) + real(amrex_real), intent(in ) :: Ez (zlo(1):zhi(1),zlo(2):zhi(2)) + real(amrex_real), intent(inout) :: dive(dlo(1):dhi(1),dlo(2):dhi(2)) + + integer :: i,k + real(amrex_real) :: dxinv(3) + real(amrex_real) :: ru, rd + + dxinv = 1.d0/dx + + do k = lo(2), hi(2) + do i = lo(1), hi(1) + if (i == 0 .and. rmin == 0.) then + ! the bulk equation diverges on axis + ! (due to the 1/r terms). The following expressions regularize + ! these divergences. + dive(i,k) = 4.*dxinv(1) * Ex(i,k) & + + dxinv(3) * (Ez(i,k) - Ez(i,k-1)) + else + ru = 1.d0 + 0.5d0/(rmin*dxinv(1) + i) + rd = 1.d0 - 0.5d0/(rmin*dxinv(1) + i) + dive(i,k) = dxinv(1) * (ru*Ex(i,k) - rd*Ex(i-1,k)) & + + dxinv(3) * (Ez(i,k) - Ez(i,k-1)) + end if + end do + end do + end subroutine warpx_compute_dive_rz + + subroutine warpx_sync_current_2d (lo, hi, crse, clo, chi, fine, flo, fhi, dir) & bind(c, name='warpx_sync_current_2d') integer, intent(in) :: lo(2), hi(2), flo(2), fhi(2), clo(2), chi(2), dir diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index f246f9f54..3d92b7651 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -48,7 +48,6 @@ #elif (AMREX_SPACEDIM == 2) #define WRPX_COMPUTE_DIVB warpx_compute_divb_2d -#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d #define WRPX_SYNC_CURRENT warpx_sync_current_2d #define WRPX_SYNC_RHO warpx_sync_rho_2d @@ -74,6 +73,12 @@ #define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs #define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_2d +#ifdef WARPX_RZ +#define WRPX_COMPUTE_DIVE warpx_compute_dive_rz +#else +#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d +#endif + #endif #ifdef __cplusplus @@ -102,6 +107,12 @@ extern "C" const long* nox, const long* noy,const long* noz, const long* lvect, const long* charge_depo_algo); + // Charge deposition finalize for RZ + void warpx_charge_deposition_rz_volume_scaling( + amrex::Real* rho, const long* rho_ng, const int* rho_ntot, + const amrex::Real* rmin, + const amrex::Real* dr); + // Current deposition void warpx_current_deposition( amrex::Real* jx, const long* jx_ng, const int* jx_ntot, @@ -362,7 +373,11 @@ extern "C" const BL_FORT_FAB_ARG_ANYD(ex), const BL_FORT_FAB_ARG_ANYD(ey), const BL_FORT_FAB_ARG_ANYD(ez), - const amrex::Real* dx); + const amrex::Real* dx +#ifdef WARPX_RZ + ,const amrex::Real* rmin +#endif + ); void WRPX_PUSH_PML_BVEC(const int* xlo, const int* xhi, const int* ylo, const int* yhi, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 151342ff5..6d6a2e411 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -18,7 +18,8 @@ #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2drz_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz -#define WRPX_PXR_RZ_VOLUME_SCALING apply_rz_volume_scaling +#define WRPX_PXR_RZ_VOLUME_SCALING_RHO apply_rz_volume_scaling_rho +#define WRPX_PXR_RZ_VOLUME_SCALING_J apply_rz_volume_scaling_j #define WRPX_PXR_PUSH_BVEC pxrpush_emrz_bvec #define WRPX_PXR_PUSH_EVEC pxrpush_emrz_evec @@ -227,9 +228,15 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! Dimension 2 #elif (AMREX_SPACEDIM==2) +#ifdef WARPX_RZ + logical(pxr_logical) :: l_2drz = .TRUE._c_long +#else + logical(pxr_logical) :: l_2drz = .FALSE._c_long +#endif + CALL pxr_depose_rho_n_2dxz(rho,np,xp,yp,zp,w,q,xmin,zmin,dx,dz,nx,nz,& nxguard,nzguard,nox,noz, & - .TRUE._c_long, .FALSE._c_long, .FALSE._c_long, 0_c_long) + .TRUE._c_long, .FALSE._c_long, l_2drz, 0_c_long) #endif @@ -238,6 +245,44 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! _________________________________________________________________ !> !> @brief + !> Applies the inverse volume scaling for RZ charge deposition + !> + !> @details + !> The scaling is done for both single mode (FDTD) and + !> multi mode (spectral) (todo) + ! + !> @param[inout] rho charge array + !> @param[in] rmin tile grid minimum radius + !> @param[in] dr radial space discretization steps + !> @param[in] nx,ny,nz number of cells + !> @param[in] nxguard,nyguard,nzguard number of guard cells + !> + subroutine warpx_charge_deposition_rz_volume_scaling(rho,rho_ng,rho_ntot,rmin,dr) & + bind(C, name="warpx_charge_deposition_rz_volume_scaling") + + integer, intent(in) :: rho_ntot(AMREX_SPACEDIM) + integer(c_long), intent(in) :: rho_ng + real(amrex_real), intent(IN OUT):: rho(*) + real(amrex_real), intent(IN) :: rmin, dr + + integer(c_long) :: type_rz_depose = 1 + + ! Compute the number of valid cells and guard cells + integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) + rho_nvalid = rho_ntot - 2*rho_ng + rho_nguards = rho_ng + +#ifdef WARPX_RZ + CALL WRPX_PXR_RZ_VOLUME_SCALING_RHO( & + rho,rho_nguards,rho_nvalid, & + rmin,dr,type_rz_depose) +#endif + + end subroutine warpx_charge_deposition_rz_volume_scaling + + ! _________________________________________________________________ + !> + !> @brief !> Main subroutine for the current deposition !> !> @details @@ -353,7 +398,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jz_nguards = jz_ng #ifdef WARPX_RZ - CALL WRPX_PXR_RZ_VOLUME_SCALING( & + CALL WRPX_PXR_RZ_VOLUME_SCALING_J( & jx,jx_nguards,jx_nvalid, & jy,jy_nguards,jy_nvalid, & jz,jz_nguards,jz_nvalid, & @@ -413,7 +458,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n END SELECT !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) #elif (AMREX_SPACEDIM == 2) CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) @@ -498,7 +543,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) #elif (AMREX_SPACEDIM == 2) CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index b825924c6..edcf402c9 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -1,4 +1,5 @@ CEXE_sources += CustomDensityProb.cpp +CEXE_sources += PlasmaProfiles.cpp CEXE_sources += WarpXInitData.cpp CEXE_sources += CustomMomentumProb.cpp CEXE_sources += PlasmaInjector.cpp diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index f6b76e8b6..cd5802a55 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -10,6 +10,8 @@ #include "AMReX_ParmParse.H" #include "AMReX_Utility.H" +enum class predefined_profile_flag { Null, parabolic_channel }; + /// /// PlasmaDensityProfile describes how the charge density /// is set in particle initialization. Subclasses must define a @@ -59,6 +61,24 @@ private: }; /// +/// This describes predefined density distributions. +/// +class PredefinedDensityProfile : public PlasmaDensityProfile +{ +public: + PredefinedDensityProfile(const std::string& species_name); + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const override; + amrex::Real ParabolicChannel(amrex::Real x, + amrex::Real y, + amrex::Real z) const; +private: + predefined_profile_flag which_profile = predefined_profile_flag::Null; + amrex::Vector<amrex::Real> params; +}; + +/// /// This describes a density function parsed in the input file. /// class ParseDensityProfile : public PlasmaDensityProfile diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 0da9318de..52b5d8b61 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -70,6 +70,17 @@ CustomDensityProfile::CustomDensityProfile(const std::string& species_name) pp.getarr("custom_profile_params", params); } +PredefinedDensityProfile::PredefinedDensityProfile(const std::string& species_name) +{ + ParmParse pp(species_name); + std::string which_profile_s; + pp.getarr("predefined_profile_params", params); + pp.query("predefined_profile_name", which_profile_s); + if (which_profile_s == "parabolic_channel"){ + which_profile = predefined_profile_flag::parabolic_channel; + } +} + ParseDensityProfile::ParseDensityProfile(std::string parse_density_function) : _parse_density_function(parse_density_function) { @@ -333,6 +344,8 @@ void PlasmaInjector::parseDensity(ParmParse pp){ rho_prof.reset(new ConstantDensityProfile(density)); } else if (rho_prof_s == "custom") { rho_prof.reset(new CustomDensityProfile(species_name)); + } else if (rho_prof_s == "predefined") { + rho_prof.reset(new PredefinedDensityProfile(species_name)); } else if (rho_prof_s == "parse_density_function") { pp.get("density_function(x,y,z)", str_density_function); rho_prof.reset(new ParseDensityProfile(str_density_function)); diff --git a/Source/Initialization/PlasmaProfiles.cpp b/Source/Initialization/PlasmaProfiles.cpp new file mode 100644 index 000000000..d9d207f7e --- /dev/null +++ b/Source/Initialization/PlasmaProfiles.cpp @@ -0,0 +1,41 @@ +#include <PlasmaInjector.H> +#include <cmath> +#include <iostream> +#include <WarpXConst.H> + +using namespace amrex; + +Real PredefinedDensityProfile::getDensity(Real x, Real y, Real z) const { + Real n; + if ( which_profile == predefined_profile_flag::parabolic_channel ) { + n = ParabolicChannel(x,y,z); + } + return n; +} + +/// +/// plateau between linear upramp and downramp, and parab transverse profile +/// +Real PredefinedDensityProfile::ParabolicChannel(Real x, Real y, Real z) const { + // params = [z_start ramp_up plateau ramp_down rc n0] + Real z_start = params[0]; + Real ramp_up = params[1]; + Real plateau = params[2]; + Real ramp_down = params[3]; + Real rc = params[4]; + Real n0 = params[5]; + Real n; + Real kp = PhysConst::q_e/PhysConst::c*sqrt( n0/(PhysConst::m_e*PhysConst::ep0) ); + + if ((z-z_start)>=0 and (z-z_start)<ramp_up ) { + n = (z-z_start)/ramp_up; + } else if ((z-z_start)>=ramp_up and (z-z_start)<ramp_up+plateau ) { + n = 1; + } else if ((z-z_start)>=ramp_up+plateau and (z-z_start)<ramp_up+plateau+ramp_down) { + n = 1-((z-z_start)-ramp_up-plateau)/ramp_down; + } else { + n = 0; + } + n *= n0*(1+4*(x*x+y*y)/(kp*kp*std::pow(rc,4))); + return n; +} diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package index cbb1b5234..3d1fcf1da 100644 --- a/Source/Parallelization/Make.package +++ b/Source/Parallelization/Make.package @@ -1,5 +1,6 @@ CEXE_sources += WarpXComm.cpp CEXE_sources += WarpXRegrid.cpp +CEXE_headers += WarpXSumGuardCells.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parallelization diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 28971eb0c..348b31c8b 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1,5 +1,6 @@ #include <WarpX.H> #include <WarpX_f.H> +#include <WarpXSumGuardCells.H> #include <AMReX_FillPatchUtil_F.H> @@ -355,7 +356,8 @@ WarpX::SyncCurrent () { BL_PROFILE("SyncCurrent()"); - // Restrict fine patch current onto the coarse patch, before fine patch SumBoundary + // Restrict fine patch current onto the coarse patch, before + // summing the guard cells of the fine patch for (int lev = 1; lev <= finest_level; ++lev) { current_cp[lev][0]->setVal(0.0); @@ -391,8 +393,9 @@ WarpX::SyncCurrent () bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); // Then swap j_fp and current_fp std::swap(j_fp[lev][idim], current_fp[lev][idim]); - // At this point, current_fp may have false values close to the - // edges of each FAB. This will be solved with a SumBoundary later. + // At this point, current_fp may have false values close to the + // edges of each FAB. This will be solved with when summing + // the guard cells later. // j_fp contains the exact MultiFab current_fp before this step. } } @@ -427,11 +430,11 @@ WarpX::SyncCurrent () { const auto& period = Geom(lev).periodicity(); // When using a bilinear filter with many passes, current_fp has - // temporarily more ghost cells here, so that its value inside + // temporarily more ghost cells here, so that its value inside // the domain is correct at the end of this stage. - current_fp[lev][0]->SumBoundary(period); - current_fp[lev][1]->SumBoundary(period); - current_fp[lev][2]->SumBoundary(period); + WarpXSumGuardCells(*(current_fp[lev][0]),period); + WarpXSumGuardCells(*(current_fp[lev][1]),period); + WarpXSumGuardCells(*(current_fp[lev][2]),period); } // Add fine level's coarse patch to coarse level's fine patch @@ -461,9 +464,9 @@ WarpX::SyncCurrent () for (int lev = 1; lev <= finest_level; ++lev) { const auto& cperiod = Geom(lev-1).periodicity(); - current_cp[lev][0]->SumBoundary(cperiod); - current_cp[lev][1]->SumBoundary(cperiod); - current_cp[lev][2]->SumBoundary(cperiod); + WarpXSumGuardCells(*(current_cp[lev][0]),cperiod); + WarpXSumGuardCells(*(current_cp[lev][1]),cperiod); + WarpXSumGuardCells(*(current_cp[lev][2]),cperiod); } if (WarpX::use_filter) { @@ -476,7 +479,7 @@ WarpX::SyncCurrent () std::swap(j_fp[lev][idim], current_fp[lev][idim]); // Then copy the interior of j_fp to current_fp. MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, 1, 0); - // current_fp has right number of ghost cells and + // current_fp has right number of ghost cells and // correct filtered values here. } } @@ -556,7 +559,8 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, { if (!rhof[0]) return; - // Restrict fine patch onto the coarse patch, before fine patch SumBoundary + // Restrict fine patch onto the coarse patch, + // before summing the guard cells of the fine patch for (int lev = 1; lev <= finest_level; ++lev) { rhoc[lev]->setVal(0.0); @@ -607,7 +611,7 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, for (int lev = 0; lev <= finest_level; ++lev) { const auto& period = Geom(lev).periodicity(); - rhof[lev]->SumBoundary(period); + WarpXSumGuardCells( *(rhof[lev]), period, 0, rhof[lev]->nComp() ); } // Add fine level's coarse patch to coarse level's fine patch @@ -631,7 +635,7 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { const auto& cperiod = Geom(lev-1).periodicity(); - rhoc[lev]->SumBoundary(cperiod); + WarpXSumGuardCells( *(rhoc[lev]), cperiod, 0, rhoc[lev]->nComp() ); } if (WarpX::use_filter) { @@ -729,10 +733,9 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *j[idim]); - jf.SumBoundary(period); - MultiFab::Copy(*j[idim], jf, 0, 0, 1, 0); + WarpXSumGuardCells(*(j[idim]), jf, period); } else { - j[idim]->SumBoundary(period); + WarpXSumGuardCells(*(j[idim]), period); } } } @@ -780,8 +783,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) MultiFab::Add(jfb, jfc, 0, 0, 1, ng); mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); - jfc.SumBoundary(period); - MultiFab::Copy(*current_cp[lev+1][idim], jfc, 0, 0, 1, 0); + WarpXSumGuardCells(*current_cp[lev+1][idim], jfc, period); } else if (use_filter) // but no buffer { @@ -792,8 +794,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_cp[lev+1][idim]->DistributionMap(), 1, ng); bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]); mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); - jf.SumBoundary(period); - MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0); + WarpXSumGuardCells(*current_cp[lev+1][idim], jf, period); } else if (current_buf[lev+1][idim]) // but no filter { @@ -803,14 +804,14 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, 1, current_buf[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); - current_cp[lev+1][idim]->SumBoundary(period); + WarpXSumGuardCells(*(current_cp[lev+1][idim]), period); } else // no filter, no buffer { mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); - current_cp[lev+1][idim]->SumBoundary(period); + WarpXSumGuardCells(*(current_cp[lev+1][idim]), period); } MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } @@ -840,10 +841,9 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i ng += bilinear_filter.stencil_length_each_dir-1; MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp); - rf.SumBoundary(period); - MultiFab::Copy(*r, rf, 0, icomp, ncomp, 0); + WarpXSumGuardCells(*r, rf, period, icomp, ncomp ); } else { - r->SumBoundary(icomp, ncomp, period); + WarpXSumGuardCells(*r, period, icomp, ncomp); } } @@ -885,9 +885,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) MultiFab::Add(rhofb, rhofc, 0, 0, ncomp, ng); mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); - - rhofc.SumBoundary(period); - MultiFab::Copy(*rho_cp[lev+1], rhofc, 0, 0, ncomp, 0); + WarpXSumGuardCells( *rho_cp[lev+1], rhofc, period, icomp, ncomp ); } else if (use_filter) // but no buffer { @@ -896,8 +894,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp); mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); - rf.SumBoundary(0, ncomp, period); - MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); + WarpXSumGuardCells( *rho_cp[lev+1], rf, period, icomp, ncomp ); } else if (charge_buf[lev+1]) // but no filter { @@ -908,14 +905,14 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) ncomp, charge_buf[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); - rho_cp[lev+1]->SumBoundary(icomp, ncomp, period); + WarpXSumGuardCells(*(rho_cp[lev+1]), period, icomp, ncomp); } else // no filter, no buffer { mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); - rho_cp[lev+1]->SumBoundary(icomp, ncomp, period); + WarpXSumGuardCells(*(rho_cp[lev+1]), period, icomp, ncomp); } ApplyFilterandSumBoundaryRho(lev, PatchType::fine, icomp, ncomp); MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0); diff --git a/Source/Parallelization/WarpXSumGuardCells.H b/Source/Parallelization/WarpXSumGuardCells.H new file mode 100644 index 000000000..24ad1b80f --- /dev/null +++ b/Source/Parallelization/WarpXSumGuardCells.H @@ -0,0 +1,61 @@ +#ifndef WARPX_SUM_GUARD_CELLS_H_ +#define WARPX_SUM_GUARD_CELLS_H_ + +#include <AMReX_MultiFab.H> + +/* \brief Sum the values of `mf`, where the different boxes overlap + * (i.e. in the guard cells) + * + * This is typically called for the sources of the Maxwell equations (J/rho) + * after deposition from the macroparticles. + * + * - When WarpX is compiled with a finite-difference scheme: this only + * updates the *valid* cells of `mf` + * - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this + * updates both the *valid* cells and *guard* cells. (This is because a + * spectral solver requires the value of the sources over a large stencil.) + */ +void +WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period, + const int icomp=0, const int ncomp=1){ +#ifdef WARPX_USE_PSATD + // Update both valid cells and guard cells + const amrex::IntVect n_updated_guards = mf.nGrowVect(); +#else + // Update only the valid cells + const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector(); +#endif + mf.SumBoundary(icomp, ncomp, n_updated_guards, period); +} + +/* \brief Sum the values of `src` where the different boxes overlap + * (i.e. in the guard cells) and copy them into `dst` + * + * This is typically called for the sources of the Maxwell equations (J/rho) + * after deposition from the macroparticles + filtering. + * + * - When WarpX is compiled with a finite-difference scheme: this only + * updates the *valid* cells of `dst` + * - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this + * updates both the *valid* cells and *guard* cells. (This is because a + * spectral solver requires the value of the sources over a large stencil.) + * + * Note: `i_comp` is the component where the results will be stored in `dst`; + * The component from which we copy in `src` is always 0. + */ +void +WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src, + const amrex::Periodicity& period, + const int icomp=0, const int ncomp=1){ +#ifdef WARPX_USE_PSATD + // Update both valid cells and guard cells + const amrex::IntVect n_updated_guards = dst.nGrowVect(); +#else + // Update only the valid cells + const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector(); +#endif + src.SumBoundary(icomp, ncomp, n_updated_guards, period); + amrex::Copy( dst, src, 0, icomp, ncomp, n_updated_guards ); +} + +#endif // WARPX_SUM_GUARD_CELLS_H_ diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 695faaa62..80f3882a0 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -641,6 +641,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); #ifndef AMREX_USE_GPU @@ -696,6 +701,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &cxyzmin_tile[0], &cdx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); #ifndef AMREX_USE_GPU @@ -852,6 +862,12 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + long ngRho = WarpX::nox; + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif #ifdef _OPENMP rhofab.atomicAdd(local_rho); diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index a973c2283..39fded8e8 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -1,7 +1,11 @@ #include <AMReX_REAL.H> #include <AMReX_Vector.H> +#include <AMReX_MultiFab.H> void ReadBoostedFrameParameters(amrex::Real& gamma_boost, amrex::Real& beta_boost, amrex::Vector<int>& boost_direction); void ConvertLabParamsToBoost(); + +void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, + amrex::Real zmax); diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4a884330a..4ec7ebb51 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -3,6 +3,7 @@ #include <WarpXUtil.H> #include <WarpXConst.H> #include <AMReX_ParmParse.H> +#include <WarpX.H> using namespace amrex; @@ -95,3 +96,44 @@ void ConvertLabParamsToBoost() pp_wpx.addarr("fine_tag_hi", fine_tag_hi); } } + +/* \brief Function that sets the value of MultiFab MF to zero for z between + * zmin and zmax. + */ +void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax){ + BL_PROFILE("WarpX::NullifyMF()"); +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for(amrex::MFIter mfi(mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi){ + const amrex::Box& bx = mfi.tilebox(); + // Get box lower and upper physical z bound, and dz + const amrex::Real zmin_box = WarpX::LowerCorner(bx, lev)[2]; + const amrex::Real zmax_box = WarpX::UpperCorner(bx, lev)[2]; + amrex::Real dz = WarpX::CellSize(lev)[2]; + // Get box lower index in the z direction +#if (AMREX_SPACEDIM==3) + const int lo_ind = bx.loVect()[2]; +#else + const int lo_ind = bx.loVect()[1]; +#endif + // Check if box intersect with [zmin, zmax] + if ( (zmin>zmin_box && zmin<=zmax_box) || + (zmax>zmin_box && zmax<=zmax_box) ){ + Array4<Real> arr = mf[mfi].array(); + // Set field to 0 between zmin and zmax + ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept{ +#if (AMREX_SPACEDIM==3) + const Real z_gridpoint = zmin_box+(k-lo_ind)*dz; +#else + const Real z_gridpoint = zmin_box+(j-lo_ind)*dz; +#endif + if ( (z_gridpoint >= zmin) && (z_gridpoint < zmax) ) { + arr(i,j,k) = 0.; + }; + } + ); + } + } +} diff --git a/Source/WarpX.H b/Source/WarpX.H index 583820646..6cec8e5be 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -144,6 +144,13 @@ public: static amrex::IntVect filter_npass_each_dir; BilinearFilter bilinear_filter; + static int num_mirrors; + amrex::Vector<amrex::Real> mirror_z; + amrex::Vector<amrex::Real> mirror_z_width; + amrex::Vector<int> mirror_z_npoints; + + void applyMirrors(amrex::Real time); + void ComputeDt (); int MoveWindow (bool move_j); void UpdatePlasmaInjectionPosition (amrex::Real dt); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 532858556..47ead98df 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -50,6 +50,8 @@ bool WarpX::use_filter = false; bool WarpX::serialize_ics = false; bool WarpX::refine_plasma = false; +int WarpX::num_mirrors = 0; + int WarpX::sort_int = -1; bool WarpX::do_boosted_frame_diagnostic = false; @@ -355,6 +357,16 @@ WarpX::ReadParameters () filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; #endif + pp.query("num_mirrors", num_mirrors); + if (num_mirrors>0){ + mirror_z.resize(num_mirrors); + pp.getarr("mirror_z", mirror_z, 0, num_mirrors); + mirror_z_width.resize(num_mirrors); + pp.getarr("mirror_z_width", mirror_z_width, 0, num_mirrors); + mirror_z_npoints.resize(num_mirrors); + pp.getarr("mirror_z_npoints", mirror_z_npoints, 0, num_mirrors); + } + pp.query("serialize_ics", serialize_ics); pp.query("refine_plasma", refine_plasma); pp.query("do_dive_cleaning", do_dive_cleaning); @@ -975,12 +987,19 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, for (MFIter mfi(divE, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); +#ifdef WARPX_RZ + const Real xmin = bx.smallEnd(0)*dx[0]; +#endif WRPX_COMPUTE_DIVE(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_N_ANYD(divE[mfi],dcomp), BL_TO_FORTRAN_ANYD((*E[0])[mfi]), BL_TO_FORTRAN_ANYD((*E[1])[mfi]), BL_TO_FORTRAN_ANYD((*E[2])[mfi]), - dx.data()); + dx.data() +#ifdef WARPX_RZ + ,&xmin +#endif + ); } } @@ -995,12 +1014,19 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, for (MFIter mfi(divE, true); mfi.isValid(); ++mfi) { Box bx = mfi.growntilebox(ngrow); +#ifdef WARPX_RZ + const Real xmin = bx.smallEnd(0)*dx[0]; +#endif WRPX_COMPUTE_DIVE(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_N_ANYD(divE[mfi],dcomp), BL_TO_FORTRAN_ANYD((*E[0])[mfi]), BL_TO_FORTRAN_ANYD((*E[1])[mfi]), BL_TO_FORTRAN_ANYD((*E[2])[mfi]), - dx.data()); + dx.data() +#ifdef WARPX_RZ + ,&xmin +#endif + ); } } |