diff options
Diffstat (limited to 'Source')
31 files changed, 829 insertions, 181 deletions
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 1c5dc0104..02c21529b 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) ); } } @@ -85,8 +85,8 @@ WarpX::EvolveEM (int numsteps) amrex::Print() << " in evolve after pushP \n"; } 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(); @@ -104,6 +104,10 @@ WarpX::EvolveEM (int numsteps) } amrex::Print() << " in evolve after onestep no sub \n"; + if (num_mirrors>0){ + applyMirrors(cur_time); + } + #ifdef WARPX_USE_PY if (warpx_py_beforeEsolve) warpx_py_beforeEsolve(); #endif @@ -112,8 +116,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; } @@ -125,18 +131,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; @@ -146,11 +152,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"; @@ -160,10 +166,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; @@ -173,7 +179,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(); @@ -185,34 +191,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); @@ -225,8 +231,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) @@ -236,8 +244,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) { @@ -265,8 +274,9 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_particlescraper) warpx_py_particlescraper(); if (warpx_py_beforedeposition) warpx_py_beforedeposition(); #endif - + amrex::Print() << " before push \n"; PushParticlesandDepose(cur_time); + amrex::Print() << " after push \n"; #ifdef WARPX_USE_PY if (warpx_py_afterdeposition) warpx_py_afterdeposition(); @@ -473,23 +483,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); @@ -504,3 +514,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/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/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 7cf31b97b..34743525e 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -3,11 +3,12 @@ #include <SpectralKSpace.H> #include <SpectralFieldData.H> +#include <SpectralBaseAlgorithm.H> /* \brief Class that updates the field in spectral space * and stores the coefficients of the corresponding update equation. */ -class PsatdAlgorithm +class PsatdAlgorithm : public SpectralBaseAlgorithm { using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 9e0bbfe06..8dd2a830f 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/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 5dd94d0c7..516e368ab 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -84,11 +84,19 @@ private: std::string field_function; // laser particle domain - amrex::RealBox prob_domain; - + amrex::RealBox laser_injection_box; + // Theoretical position of the antenna. Used if do_continuous_injection=1. + // Track the position of the antenna until it enters the simulation domain. + amrex::Vector<amrex::Real> updated_position; + void ComputeSpacing (int lev, amrex::Real& Sx, amrex::Real& Sy) const; void ComputeWeightMobility (amrex::Real Sx, amrex::Real Sy); void InitData (int lev); + // Inject the laser antenna during the simulation, if it started + // outside of the simulation domain and enters it. + void ContinuousInjection(const amrex::RealBox& injection_box) override; + // Update position of the antenna + void UpdateContinuousInjectionPosition(amrex::Real dt) override; }; #endif diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 2b56c3cfd..8a2590e5e 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -49,6 +49,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, pp.query("pusher_algo", pusher_algo); pp.get("wavelength", wavelength); pp.get("e_max", e_max); + pp.query("do_continuous_injection", do_continuous_injection); if ( profile == laser_t::Gaussian ) { // Parse the properties of the Gaussian profile @@ -148,16 +149,94 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, u_Y = {0., 1., 0.}; #endif - prob_domain = Geometry::ProbDomain(); + laser_injection_box= Geometry::ProbDomain(); { Vector<Real> lo, hi; if (pp.queryarr("prob_lo", lo)) { - prob_domain.setLo(lo); + laser_injection_box.setLo(lo); } if (pp.queryarr("prob_hi", hi)) { - prob_domain.setHi(hi); + laser_injection_box.setHi(hi); } } + + if (do_continuous_injection){ + // If laser antenna initially outside of the box, store its theoretical + // position in z_antenna_th + updated_position = position; + + // Sanity checks + int dir = WarpX::moving_window_dir; + std::vector<Real> windir(3, 0.0); +#if (AMREX_SPACEDIM==2) + windir[2*dir] = 1.0; +#else + windir[dir] = 1.0; +#endif + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (nvec[0]-windir[0]) + (nvec[1]-windir[1]) + (nvec[2]-windir[2]) + < 1.e-12, "do_continous_injection for laser particle only works" + + " if moving window direction and laser propagation direction are the same"); + if ( WarpX::gamma_boost>1 ){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + + (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + + (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, + "do_continous_injection for laser particle only works if " + + "warpx.boost_direction = z. TODO: all directions."); + } + } +} + +/* \brief Check if laser particles enter the box, and inject if necessary. + * \param injection_box: a RealBox where particles should be injected. + */ +void +LaserParticleContainer::ContinuousInjection (const RealBox& injection_box) +{ + // Input parameter injection_box contains small box where injection + // should occur. + // So far, LaserParticleContainer::laser_injection_box contains the + // outdated full problem domain at t=0. + + // Convert updated_position to Real* to use RealBox::contains(). +#if (AMREX_SPACEDIM == 3) + const Real* p_pos = updated_position.dataPtr(); +#else + const Real p_pos[2] = {updated_position[0], updated_position[2]}; +#endif + if ( injection_box.contains(p_pos) ){ + // Update laser_injection_box with current value + laser_injection_box = injection_box; + // Inject laser particles. LaserParticleContainer::InitData + // is called only once, when the antenna enters the simulation + // domain. + InitData(); + } +} + +/* \brief update position of the antenna if running in boosted frame. + * \param dt time step (level 0). + * The up-to-date antenna position is stored in updated_position. + */ +void +LaserParticleContainer::UpdateContinuousInjectionPosition(Real dt) +{ + int dir = WarpX::moving_window_dir; + if (do_continuous_injection and (WarpX::gamma_boost > 1)){ + // In boosted-frame simulations, the antenna has moved since the last + // call to this function, and injection position needs to be updated +#if ( AMREX_SPACEDIM == 3 ) + updated_position[dir] -= WarpX::beta_boost * + WarpX::boost_direction[dir] * PhysConst::c * dt; +#elif ( AMREX_SPACEDIM == 2 ) + // In 2D, dir=0 corresponds to x and dir=1 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for both 2D and 3D simulations. + updated_position[2*dir] -= WarpX::beta_boost * + WarpX::boost_direction[2*dir] * PhysConst::c * dt; +#endif + } } void @@ -175,6 +254,13 @@ LaserParticleContainer::InitData (int lev) ComputeSpacing(lev, S_X, S_Y); ComputeWeightMobility(S_X, S_Y); + // LaserParticleContainer::position contains the initial position of the + // laser antenna. In the boosted frame, the antenna is moving. + // Update its position with updated_position. + if (do_continuous_injection){ + position = updated_position; + } + auto Transform = [&](int i, int j) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) return { position[0] + (S_X*(i+0.5))*u_X[0] + (S_Y*(j+0.5))*u_Y[0], @@ -210,8 +296,8 @@ LaserParticleContainer::InitData (int lev) plane_hi[1] = std::max(plane_hi[1], j); }; - const Real* prob_lo = prob_domain.lo(); - const Real* prob_hi = prob_domain.hi(); + const Real* prob_lo = laser_injection_box.lo(); + const Real* prob_hi = laser_injection_box.hi(); #if (AMREX_SPACEDIM == 3) compute_min_max(prob_lo[0], prob_lo[1], prob_lo[2]); compute_min_max(prob_hi[0], prob_lo[1], prob_lo[2]); @@ -272,7 +358,7 @@ LaserParticleContainer::InitData (int lev) #else const Real x[2] = {pos[0], pos[2]}; #endif - if (prob_domain.contains(x)) + if (laser_injection_box.contains(x)) { for (int k = 0; k<2; ++k) { particle_x.push_back(pos[0]); 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 23d97133b..e9d814a94 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> @@ -360,7 +361,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); @@ -396,8 +398,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. } } @@ -432,11 +435,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 @@ -466,9 +469,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) { @@ -481,7 +484,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. } } @@ -561,7 +564,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); @@ -612,7 +616,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 @@ -636,7 +640,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) { @@ -734,10 +738,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); } } } @@ -785,8 +788,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 { @@ -797,8 +799,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 { @@ -808,14 +809,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); } @@ -845,10 +846,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); } } @@ -890,9 +890,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 { @@ -901,8 +899,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 { @@ -913,14 +910,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/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index f3fd522a9..0c5e49c04 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -169,7 +169,15 @@ public: const amrex::Real z_old, const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, amrex::Vector<WarpXParticleContainer::DiagnosticParticleData>& parts) const; - + + // Inject particles during the simulation (for particles entering the + // simulation domain after some iterations, due to flowing plasma and/or + // moving window). + void ContinuousInjection(const amrex::RealBox& injection_box) const; + // Update injection position for continuously-injected species. + void UpdateContinuousInjectionPosition(amrex::Real dt) const; + int doContinuousInjection() const; + // // Parameters for the Cherenkov corrector in the FDTD solver. // Both stencils are calculated ar runtime. diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index a4df1f83a..440906348 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -414,3 +414,48 @@ MultiParticleContainer } } } + +/* \brief Continuous injection for particles initially outside of the domain. + * \param injection_box: Domain where new particles should be injected. + * Loop over all WarpXParticleContainer in MultiParticleContainer and + * calls virtual function ContinuousInjection. + */ +void +MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const +{ + for (int i=0; i<nspecies+nlasers; i++){ + auto& pc = allcontainers[i]; + if (pc->do_continuous_injection){ + pc->ContinuousInjection(injection_box); + } + } +} + +/* \brief Update position of continuous injection parameters. + * \param dt: simulation time step (level 0) + * All classes inherited from WarpXParticleContainer do not have + * a position to update (PhysicalParticleContainer does not do anything). + */ +void +MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const +{ + for (int i=0; i<nspecies+nlasers; i++){ + auto& pc = allcontainers[i]; + if (pc->do_continuous_injection){ + pc->UpdateContinuousInjectionPosition(dt); + } + } +} + +int +MultiParticleContainer::doContinuousInjection() const +{ + int warpx_do_continuous_injection = 0; + for (int i=0; i<nspecies+nlasers; i++){ + auto& pc = allcontainers[i]; + if (pc->do_continuous_injection){ + warpx_do_continuous_injection = 1; + } + } + return warpx_do_continuous_injection; +} diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 362683879..4f365768b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -101,8 +101,6 @@ public: const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; - bool injected = false; - protected: std::string species_name; @@ -122,6 +120,9 @@ protected: int GetRefineFac(const amrex::Real x, const amrex::Real y, const amrex::Real z); std::unique_ptr<amrex::IArrayBox> m_refined_injection_mask = nullptr; + // Inject particles during the whole simulation + void ContinuousInjection(const amrex::RealBox& injection_box) override; + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 17e6d98d9..a1d4e1319 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -24,7 +24,7 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -81,6 +81,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_backward_propagation", do_backward_propagation); pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); + pp.query("do_continuous_injection", do_continuous_injection); } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) @@ -361,7 +362,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -602,7 +603,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -1331,7 +1332,9 @@ PhysicalParticleContainer::Evolve (int lev, pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); + amrex::Print() << " before deposit chage \n"; if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); + amrex::Print() << " after deposit chage \n"; if (! do_not_push) { @@ -2004,3 +2007,14 @@ int PhysicalParticleContainer::GetRefineFac(const Real x, const Real y, const Re return ref_fac; } + +/* \brief Inject particles during the simulation + * \param injection_box: domain where particles should be injected. + */ +void +PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) +{ + // Inject plasma on level 0. Paticles will be redistributed. + const int lev=0; + AddPlasma(lev, injection_box); +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 275554cd8..6fa02b824 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -183,7 +183,18 @@ public: int thread_num, int lev, amrex::Real dt ); - + + // If particles start outside of the domain, ContinuousInjection + // makes sure that they are initialized when they enter the domain, and + // NOT before. Virtual function, overriden by derived classes. + // Current status: + // PhysicalParticleContainer: implemented. + // LaserParticleContainer: implemented. + // RigidInjectedParticleContainer: not implemented. + virtual void ContinuousInjection(const amrex::RealBox& injection_box) {} + // Update optional sub-class-specific injection location. + virtual void UpdateContinuousInjectionPosition(amrex::Real dt) {} + /// /// This returns the total charge for all the particles in this ParticleContainer. /// This is needed when solving Poisson's equation with periodic boundary conditions. @@ -248,6 +259,12 @@ protected: static int do_not_push; + // Whether to allow particles outside of the simulation domain to be + // initialized when they enter the domain. + // This is currently required because continuous injection does not + // support all features allowed by direct injection. + int do_continuous_injection = 0; + amrex::Vector<amrex::FArrayBox> local_rho; amrex::Vector<amrex::FArrayBox> local_jx; amrex::Vector<amrex::FArrayBox> local_jy; @@ -258,6 +275,7 @@ protected: private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; + }; #endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 94b58ca07..89e21df1c 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -608,6 +608,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array<Real, 3>& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU + amrex::Print() << " before icomp data ptr " << icomp << "\n"; data_ptr = (*rhomf)[pti].dataPtr(icomp); auto rholen = (*rhomf)[pti].length(); #else @@ -630,6 +631,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const long nz = rholen[1]-1-2*ngRho; #endif BL_PROFILE_VAR_START(blp_pxr_chd); + amrex::Print() << " before warpxcharge deposition \n"; warpx_charge_deposition(data_ptr, &np_current, m_xp[thread_num].dataPtr(), m_yp[thread_num].dataPtr(), @@ -641,6 +643,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 +703,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 +864,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/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 05e171f22..18d89951d 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -9,7 +9,7 @@ WarpX::UpdatePlasmaInjectionPosition (Real dt) { int dir = moving_window_dir; // Continuously inject plasma in new cells (by default only on level 0) - if (WarpX::do_plasma_injection and (WarpX::gamma_boost > 1)){ + if (WarpX::warpx_do_continuous_injection and (WarpX::gamma_boost > 1)){ // In boosted-frame simulations, the plasma has moved since the last // call to this function, and injection position needs to be updated current_injection_position -= WarpX::beta_boost * @@ -33,7 +33,16 @@ WarpX::MoveWindow (bool move_j) // and of the plasma injection moving_window_x += moving_window_v * dt[0]; int dir = moving_window_dir; + + // Update warpx.current_injection_position + // PhysicalParticleContainer uses this injection position UpdatePlasmaInjectionPosition( dt[0] ); + if (WarpX::warpx_do_continuous_injection){ + // Update injection position for WarpXParticleContainer in mypc. + // Nothing to do for PhysicalParticleContainers + // For LaserParticleContainer, need to update the antenna position. + mypc->UpdateContinuousInjectionPosition( dt[0] ); + } // compute the number of cells to shift on the base level Real new_lo[AMREX_SPACEDIM]; @@ -134,7 +143,7 @@ WarpX::MoveWindow (bool move_j) } // Continuously inject plasma in new cells (by default only on level 0) - if (WarpX::do_plasma_injection) { + if (WarpX::warpx_do_continuous_injection) { const int lev = 0; @@ -162,15 +171,11 @@ WarpX::MoveWindow (bool move_j) particleBox.setLo( dir, new_injection_position ); particleBox.setHi( dir, current_injection_position ); } - // Perform the injection of new particles in particleBox + if (particleBox.ok() and (current_injection_position != new_injection_position)){ - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast<PhysicalParticleContainer&>(pc); - ppc.AddPlasma(lev, particleBox); - } - // Update the injection position + // Performs continuous injection of all WarpXParticleContainer + // in mypc. + mypc->ContinuousInjection(particleBox); current_injection_position = new_injection_position; } } 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..44e84033f 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); @@ -469,10 +476,10 @@ private: amrex::Real current_injection_position = 0; // Plasma injection parameters - int do_plasma_injection = 0; + int warpx_do_continuous_injection = 0; int num_injected_species = -1; amrex::Vector<int> injected_plasma_species; - + int do_electrostatic = 0; int n_buffer = 4; amrex::Real const_dt = 0.5e-11; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 5814d6829..83788d5b9 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; @@ -143,13 +145,14 @@ WarpX::WarpX () // Particle Container mypc = std::unique_ptr<MultiParticleContainer> (new MultiParticleContainer(this)); - - if (do_plasma_injection) { - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast<PhysicalParticleContainer&>(pc); - ppc.injected = true; + warpx_do_continuous_injection = mypc->doContinuousInjection(); + if (warpx_do_continuous_injection){ + if (moving_window_v >= 0){ + // Inject particles continuously from the right end of the box + current_injection_position = geom[0].ProbHi(moving_window_dir); + } else { + // Inject particles continuously from the left end of the box + current_injection_position = geom[0].ProbLo(moving_window_dir); } } @@ -298,21 +301,6 @@ WarpX::ReadParameters () moving_window_v *= PhysConst::c; } - pp.query("do_plasma_injection", do_plasma_injection); - if (do_plasma_injection) { - pp.get("num_injected_species", num_injected_species); - injected_plasma_species.resize(num_injected_species); - pp.getarr("injected_plasma_species", injected_plasma_species, - 0, num_injected_species); - if (moving_window_v >= 0){ - // Inject particles continuously from the right end of the box - current_injection_position = geom[0].ProbHi(moving_window_dir); - } else { - // Inject particles continuously from the left end of the box - current_injection_position = geom[0].ProbLo(moving_window_dir); - } - } - pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); if (do_boosted_frame_diagnostic) { @@ -355,6 +343,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); @@ -977,12 +975,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 + ); } } @@ -997,12 +1002,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 + ); } } |