diff options
Diffstat (limited to 'Source')
31 files changed, 718 insertions, 616 deletions
diff --git a/Source/Diagnostics/BackTransformedDiagnostic.H b/Source/Diagnostics/BackTransformedDiagnostic.H index 9e24caa1b..94d557605 100644 --- a/Source/Diagnostics/BackTransformedDiagnostic.H +++ b/Source/Diagnostics/BackTransformedDiagnostic.H @@ -8,7 +8,6 @@ #include <AMReX_PlotFileUtil.H> #include <AMReX_ParallelDescriptor.H> #include <AMReX_Geometry.H> -#include <AMReX_CudaContainers.H> #include "MultiParticleContainer.H" #include "WarpXConst.H" diff --git a/Source/Evolve/WarpXEvolveES.cpp b/Source/Evolve/WarpXEvolveES.cpp index effd6ec96..7a57dfa80 100644 --- a/Source/Evolve/WarpXEvolveES.cpp +++ b/Source/Evolve/WarpXEvolveES.cpp @@ -179,30 +179,6 @@ void WarpX::zeroOutBoundary(amrex::MultiFab& input_data, bndry_data.FillBoundary(); } -void WarpX::sumFineToCrseNodal(const amrex::MultiFab& fine, - amrex::MultiFab& crse, - const amrex::Geometry& cgeom, - const amrex::IntVect& ratio) { - const BoxArray& fine_BA = fine.boxArray(); - const DistributionMapping& fine_dm = fine.DistributionMap(); - BoxArray coarsened_fine_BA = fine_BA; - coarsened_fine_BA.coarsen(ratio); - - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, 1, 0); - coarsened_fine_data.setVal(0.0); - - for (MFIter mfi(coarsened_fine_data); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - const Box& crse_box = coarsened_fine_data[mfi].box(); - const Box& fine_box = fine[mfi].box(); - WRPX_SUM_FINE_TO_CRSE_NODAL(bx.loVect(), bx.hiVect(), ratio.getVect(), - coarsened_fine_data[mfi].dataPtr(), crse_box.loVect(), crse_box.hiVect(), - fine[mfi].dataPtr(), fine_box.loVect(), fine_box.hiVect()); - } - - crse.copy(coarsened_fine_data, cgeom.periodicity(), FabArrayBase::ADD); -} - void WarpX::fixRHSForSolve(Vector<std::unique_ptr<MultiFab> >& rhs, const Vector<std::unique_ptr<FabArray<BaseFab<int> > > >& masks) const { diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index 5eb84b96a..1ecb1bff4 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -1,5 +1,4 @@ #include <AMReX_MultiFab.H> -#include <AMReX_CudaContainers.H> #ifndef WARPX_FILTER_H_ #define WARPX_FILTER_H_ diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index b51a83387..e7f02197f 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -16,11 +16,9 @@ #if (AMREX_SPACEDIM == 3) -#define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_3d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_3d #define WRPX_BUILD_MASK warpx_build_mask_3d #define WRPX_COMPUTE_E_NODAL warpx_compute_E_nodal_3d -#define WRPX_DEPOSIT_CIC warpx_deposit_cic_3d #define WRPX_INTERPOLATE_CIC warpx_interpolate_cic_3d #define WRPX_INTERPOLATE_CIC_TWO_LEVELS warpx_interpolate_cic_two_levels_3d #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d @@ -28,11 +26,9 @@ #elif (AMREX_SPACEDIM == 2) -#define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_2d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_2d #define WRPX_BUILD_MASK warpx_build_mask_2d #define WRPX_COMPUTE_E_NODAL warpx_compute_E_nodal_2d -#define WRPX_DEPOSIT_CIC warpx_deposit_cic_2d #define WRPX_INTERPOLATE_CIC warpx_interpolate_cic_2d #define WRPX_INTERPOLATE_CIC_TWO_LEVELS warpx_interpolate_cic_two_levels_2d #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_2d diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H index ba0d26fc4..f5819a908 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -47,6 +47,81 @@ private: amrex::Real m_ux_th, m_uy_th, m_uz_th; }; +// struct whose getMomentum returns momentum for 1 particle with relativistc +// drift velocity beta, from the Maxwell-Juttner distribution. Method is from +// Zenitani 2015 (Phys. Plasmas 22, 042116). +struct InjectorMomentumJuttner +{ + // Constructor whose inputs are: + // the temperature parameter theta, + // boost velocity/c beta, + // and boost direction dir respectively. + InjectorMomentumJuttner(amrex::Real t, amrex::Real b, int d) noexcept + : theta(t), beta(b), dir(d) + {} + + AMREX_GPU_HOST_DEVICE + amrex::XDim3 + getMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + { + // Sobol method for sampling MJ Speeds, + // from Zenitani 2015 (Phys. Plasmas 22, 042116). + amrex::Real x1, x2, gamma; + amrex::Real u [3]; + x1 = 0.; + gamma = 0.; + u[dir] = 0.; + // This condition is equation 10 in Zenitani, + // though x1 is defined differently. + while(u[dir]-gamma <= x1) + { + u[dir] = -theta* + std::log(amrex::Random()*amrex::Random()*amrex::Random()); + gamma = std::sqrt(1+std::pow(u[dir],2)); + x1 = theta*std::log(amrex::Random()); + } + // The following code samples a random unit vector + // and multiplies the result by speed u[dir]. + x1 = amrex::Random(); + x2 = amrex::Random(); + // Direction dir is an input parameter that sets the boost direction: + // 'x' -> d = 0, 'y' -> d = 1, 'z' -> d = 2. + u[(dir+1)%3] = 2*u[dir]*std::sqrt(x1*(1-x1))*std::sin(2*M_PI*x2); + u[(dir+2)%3] = 2*u[dir]*std::sqrt(x1*(1-x1))*std::cos(2*M_PI*x2); + // The value of dir is the boost direction to be transformed. + u[dir] = u[dir]*(2*x1-1); + x1 = amrex::Random(); + // The following condition is equtaion 32 in Zenitani, called + // The flipping method. It transforms the intergral: d3x' -> d3x + // where d3x' is the volume element for positions in the boosted frame. + // The particle positions and densities can be initialized in the + // simulation frame with this method. + // The flipping method can similarly transform any + // symmetric distribution from one reference frame to another moving at + // a relative velocity of beta. + // An equivalent alternative to this method native to WarpX + // would be to initialize the particle positions and densities in the + // frame moving at speed beta, and then perform a Lorentz transform + // on their positions and MJ sampled velocities to the simulation frame. + if(-beta*u[dir]/gamma>x1) + { + u[dir] = -u[dir]; + } + // This Lorentz transform is equation 17 in Zenitani. + // It transforms the integral d3u' -> d3u + // where d3u' is the volume element for momentum in the boosted frame. + u[dir] = 1/std::sqrt(1-pow(beta,2))*(u[dir]+gamma*beta); + // Note that if beta = 0 then the flipping method and Lorentz transform + // have no effect on the u[dir] direction. + return amrex::XDim3 {u[0],u[1],u[2]}; + } + +private: + int dir; + amrex::Real beta, theta; +}; + + // struct whose getMomentum returns momentum for 1 particle, for // radial expansion struct InjectorMomentumRadialExpansion @@ -124,6 +199,13 @@ struct InjectorMomentum object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th) { } + // This constructor stores a InjectorMomentumJuttner in union object. + InjectorMomentum (InjectorMomentumJuttner* t, + amrex::Real theta, amrex::Real beta, int dir) + : type(Type::juttner), + object(t, theta, beta, dir) + { } + // This constructor stores a InjectorMomentumCustom in union object. InjectorMomentum (InjectorMomentumCustom* t, std::string const& a_species_name) @@ -165,6 +247,10 @@ struct InjectorMomentum { return object.gaussian.getMomentum(x,y,z); } + case Type::juttner: + { + return object.juttner.getMomentum(x,y,z); + } case Type::constant: { return object.constant.getMomentum(x,y,z); @@ -186,7 +272,7 @@ struct InjectorMomentum } private: - enum struct Type { constant, custom, gaussian, radial_expansion, parser }; + enum struct Type { constant, custom, gaussian, juttner, radial_expansion, parser}; Type type; // An instance of union Object constructs and stores any one of @@ -204,6 +290,9 @@ private: amrex::Real a_uz_m, amrex::Real a_ux_th, amrex::Real a_uy_th, amrex::Real a_uz_th) noexcept : gaussian(a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th) {} + Object (InjectorMomentumJuttner*, + amrex::Real t, amrex::Real b, int dir) noexcept + : juttner(t,b,dir) {} Object (InjectorMomentumRadialExpansion*, amrex::Real u_over_r) noexcept : radial_expansion(u_over_r) {} @@ -215,6 +304,7 @@ private: InjectorMomentumConstant constant; InjectorMomentumCustom custom; InjectorMomentumGaussian gaussian; + InjectorMomentumJuttner juttner; InjectorMomentumRadialExpansion radial_expansion; InjectorMomentumParser parser; }; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index af04c6b31..32711cd31 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -269,6 +269,28 @@ void PlasmaInjector::parseMomentum (ParmParse& pp) // Construct InjectorMomentum with InjectorMomentumGaussian. inj_mom.reset(new InjectorMomentum((InjectorMomentumGaussian*)nullptr, ux_m, uy_m, uz_m, ux_th, uy_th, uz_th)); + } else if (mom_dist_s == "maxwell_juttner"){ + Real beta = 0.; + Real theta = 10.; + int dir = 0; + std::string direction = "x"; + pp.query("beta", beta); + pp.query("theta", theta); + pp.query("direction", direction); + if(direction == "x" || direction == "X"){ + dir = 0; + } else if (direction == "y" || direction == "Y"){ + dir = 1; + } else if (direction == "z" || direction == "Z"){ + dir = 2; + } else{ + std::stringstream stringstream; + stringstream << "Direction " << direction << " is not recognzied. Please enter x, y, or z."; + direction = stringstream.str(); + amrex::Abort(direction.c_str()); + } + // Construct InjectorMomentum with InjectorMomentumJuttner. + inj_mom.reset(new InjectorMomentum((InjectorMomentumJuttner*)nullptr, theta, beta, dir)); } else if (mom_dist_s == "radial_expansion") { Real u_over_r = 0.; pp.query("u_over_r", u_over_r); diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index e2a0743bc..63ace31fb 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -1,13 +1,14 @@ #ifndef WARPX_LaserParticleContainer_H_ #define WARPX_LaserParticleContainer_H_ -#include <limits> +#include <LaserProfiles.H> #include <WarpXParticleContainer.H> #include <WarpXConst.H> #include <WarpXParser.H> -enum class laser_t { Null, Gaussian, Harris, parse_field_function }; +#include <memory> +#include <limits> class LaserParticleContainer : public WarpXParticleContainer @@ -44,14 +45,6 @@ public: virtual void PostRestart () final; - void gaussian_laser_profile (const int np, amrex::Real const * AMREX_RESTRICT const Xp, - amrex::Real const * AMREX_RESTRICT const Yp, amrex::Real t, - amrex::Real * AMREX_RESTRICT const amplitude); - - void harris_laser_profile (const int np, amrex::Real const * AMREX_RESTRICT const Xp, - amrex::Real const * AMREX_RESTRICT const Yp, amrex::Real t, - amrex::Real * AMREX_RESTRICT const amplitude); - void calculate_laser_plane_coordinates (const int np, const int thread_num, amrex::Real * AMREX_RESTRICT const pplane_Xp, amrex::Real * AMREX_RESTRICT const pplane_Yp); @@ -68,17 +61,15 @@ protected: std::string laser_name; private: - // runtime paramters - laser_t profile = laser_t::Null; - amrex::Vector<amrex::Real> position; - amrex::Vector<amrex::Real> nvec; - amrex::Vector<amrex::Real> p_X; - amrex::Vector<amrex::Real> stc_direction; + amrex::Vector<amrex::Real> position; //! Coordinates of one of the point of the antenna + amrex::Vector<amrex::Real> nvec; //! Normal of the plane of the antenna + amrex::Vector<amrex::Real> p_X;// ! Polarization long pusher_algo = -1; amrex::Real e_max = std::numeric_limits<amrex::Real>::quiet_NaN(); amrex::Real wavelength = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real Z0_lab = 0; // Position of the antenna in the lab frame long min_particles_per_mode = 4; @@ -90,19 +81,6 @@ private: amrex::Real weight = std::numeric_limits<amrex::Real>::quiet_NaN(); amrex::Real mobility = std::numeric_limits<amrex::Real>::quiet_NaN(); - // Gaussian profile - amrex::Real profile_waist = std::numeric_limits<amrex::Real>::quiet_NaN(); - amrex::Real profile_duration = std::numeric_limits<amrex::Real>::quiet_NaN(); - amrex::Real profile_t_peak = std::numeric_limits<amrex::Real>::quiet_NaN(); - amrex::Real profile_focal_distance = std::numeric_limits<amrex::Real>::quiet_NaN(); - amrex::Real zeta = 0.; - amrex::Real beta = 0.; - amrex::Real phi2 = 0.; - amrex::Real theta_stc = 0.; - - // parse_field_function profile - WarpXParser parser; - std::string field_function; // laser particle domain amrex::RealBox laser_injection_box; @@ -118,6 +96,10 @@ private: void ContinuousInjection(const amrex::RealBox& injection_box) override; // Update position of the antenna void UpdateContinuousInjectionPosition(amrex::Real dt) override; + + // Unique (smart) pointer to the laser profile + std::unique_ptr<WarpXLaserProfiles::ILaserProfile> m_up_laser_profile; + }; #endif diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 9493672e0..a330200cc 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -1,4 +1,3 @@ - #include <limits> #include <cmath> #include <algorithm> @@ -11,6 +10,7 @@ #include <MultiParticleContainer.H> using namespace amrex; +using namespace WarpXLaserProfiles; namespace { @@ -34,70 +34,24 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, std::string laser_type_s; pp.get("profile", laser_type_s); std::transform(laser_type_s.begin(), laser_type_s.end(), laser_type_s.begin(), ::tolower); - if (laser_type_s == "gaussian") { - profile = laser_t::Gaussian; - } else if(laser_type_s == "harris") { - profile = laser_t::Harris; - } else if(laser_type_s == "parse_field_function") { - profile = laser_t::parse_field_function; - } else { - amrex::Abort("Unknown laser type"); - } // Parse the properties of the antenna pp.getarr("position", position); pp.getarr("direction", nvec); pp.getarr("polarization", p_X); + pp.query("pusher_algo", pusher_algo); pp.get("wavelength", wavelength); pp.get("e_max", e_max); pp.query("do_continuous_injection", do_continuous_injection); pp.query("min_particles_per_mode", min_particles_per_mode); - if ( profile == laser_t::Gaussian ) { - // Parse the properties of the Gaussian profile - pp.get("profile_waist", profile_waist); - pp.get("profile_duration", profile_duration); - pp.get("profile_t_peak", profile_t_peak); - pp.get("profile_focal_distance", profile_focal_distance); - stc_direction = p_X; - pp.queryarr("stc_direction", stc_direction); - pp.query("zeta", zeta); - pp.query("beta", beta); - pp.query("phi2", phi2); - } - - if ( profile == laser_t::Harris ) { - // Parse the properties of the Harris profile - pp.get("profile_waist", profile_waist); - pp.get("profile_duration", profile_duration); - pp.get("profile_focal_distance", profile_focal_distance); - } - - if ( profile == laser_t::parse_field_function ) { - // Parse the properties of the parse_field_function profile - pp.get("field_function(X,Y,t)", field_function); - parser.define(field_function); - parser.registerVariables({"X","Y","t"}); - - ParmParse ppc("my_constants"); - std::set<std::string> symbols = parser.symbols(); - symbols.erase("X"); - symbols.erase("Y"); - symbols.erase("t"); // after removing variables, we are left with constants - for (auto it = symbols.begin(); it != symbols.end(); ) { - Real v; - if (ppc.query(it->c_str(), v)) { - parser.setConstant(*it, v); - it = symbols.erase(it); - } else { - ++it; - } - } - for (auto const& s : symbols) { // make sure there no unknown symbols - amrex::Abort("Laser Profile: Unknown symbol "+s); - } + //Select laser profile + if(laser_profiles_dictionary.count(laser_type_s) == 0){ + amrex::Abort(std::string("Unknown laser type: ").append(laser_type_s)); } + m_up_laser_profile = laser_profiles_dictionary.at(laser_type_s)(); + //__________ // Plane normal Real s = 1.0_rt / std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); @@ -128,22 +82,6 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, p_Y = CrossProduct(nvec, p_X); // The second polarization vector - s = 1.0_rt / std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); - stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s }; - Real const dp2 = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp2) < 1.0e-14, - "stc_direction is not perpendicular to the laser plane vector"); - - // Get angle between p_X and stc_direction - // in 2d, stcs are in the simulation plane -#if AMREX_SPACEDIM == 3 - theta_stc = acos(stc_direction[0]*p_X[0] + - stc_direction[1]*p_X[1] + - stc_direction[2]*p_X[2]); -#else - theta_stc = 0.; -#endif - #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) u_X = p_X; u_Y = p_Y; @@ -189,6 +127,14 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, "warpx.boost_direction = z. TODO: all directions."); } } + + //Init laser profile + CommonLaserParameters common_params; + common_params.wavelength = wavelength; + common_params.e_max = e_max; + common_params.p_X = p_X; + common_params.nvec = nvec; + m_up_laser_profile->init(pp, ParmParse{"my_constants"}, common_params); } /* \brief Check if laser particles enter the box, and inject if necessary. @@ -454,7 +400,7 @@ LaserParticleContainer::Evolve (int lev, int const thread_num = 0; #endif - Cuda::ManagedDeviceVector<Real> plane_Xp, plane_Yp, amplitude_E; + Gpu::ManagedDeviceVector<Real> plane_Xp, plane_Yp, amplitude_E; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -503,21 +449,9 @@ LaserParticleContainer::Evolve (int lev, // Calculate the laser amplitude to be emitted, // at the position of the emission plane - if (profile == laser_t::Gaussian) { - gaussian_laser_profile(np, plane_Xp.dataPtr(), plane_Yp.dataPtr(), - t_lab, amplitude_E.dataPtr()); - } - - if (profile == laser_t::Harris) { - harris_laser_profile(np, plane_Xp.dataPtr(), plane_Yp.dataPtr(), - t_lab, amplitude_E.dataPtr()); - } - - if (profile == laser_t::parse_field_function) { - for (int i = 0; i < np; ++i) { - amplitude_E[i] = parser.eval(plane_Xp[i], plane_Yp[i], t); - } - } + m_up_laser_profile->fill_amplitude( + np, plane_Xp.dataPtr(), plane_Yp.dataPtr(), + t_lab, amplitude_E.dataPtr()); // Calculate the corresponding momentum and position for the particles update_laser_particle(np, uxp.dataPtr(), uyp.dataPtr(), diff --git a/Source/Laser/LaserProfiles.H b/Source/Laser/LaserProfiles.H new file mode 100644 index 000000000..528309492 --- /dev/null +++ b/Source/Laser/LaserProfiles.H @@ -0,0 +1,202 @@ +#ifndef WARPX_LaserProfiles_H_ +#define WARPX_LaserProfiles_H_ + +#include <AMReX_REAL.H> +#include <WarpXParser.H> +#include <AMReX_ParmParse.H> +#include <AMReX_Vector.H> + +#include <WarpXParser.H> + +#include <map> +#include <string> +#include <memory> +#include <functional> +#include <limits> + +namespace WarpXLaserProfiles { + +/** Common laser profile parameters + * + * Parameters for each laser profile as shared among all laser profile classes. + */ +struct CommonLaserParameters +{ + amrex::Real wavelength; //! central wavelength + amrex::Real e_max; //! maximum electric field at peak + amrex::Vector<amrex::Real> p_X;// ! Polarization + amrex::Vector<amrex::Real> nvec; //! Normal of the plane of the antenna +}; + + +/** Abstract interface for laser profile classes + * + * Each new laser profile should inherit this interface and implement two + * methods: init and fill_amplitude (described below). + * + * The declaration of a LaserProfile class should be placed in this file, + * while the implementation of the methods should be in a dedicated file in + * LaserProfilesImpl folder. LaserProfile classes should appear in + * laser_profiles_dictionary to be used by LaserParticleContainer. + */ +class ILaserProfile +{ +public: + /** Initialize Laser Profile + * + * Reads the section of the inputfile relative to the laser beam + * (e.g. <laser_name>.profile_t_peak, <laser_name>.profile_duration...) + * and the "my_constants" section. It also receives some common + * laser profile parameters. It uses these data to initialize the + * member variables of the laser profile class. + * + * @param[in] ppl should be amrex::ParmParse(laser_name) + * @param[in] ppc should be amrex::ParmParse("my_constants") + * @param[in] params common laser profile parameters + */ + virtual void + init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& ppc, + CommonLaserParameters params) = 0; + + /** Fill Electric Field Amplitude for each particle of the antenna. + * + * Xp, Yp and amplitude must be arrays with the same length + * + * @param[in] Xp X coordinate of the particles of the antenna + * @param[in] Yp Y coordinate of the particles of the antenna + * @param[in] t time (seconds) + * @param[out] amplitude of the electric field (V/m) + */ + virtual void + fill_amplitude ( + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) = 0; + + virtual ~ILaserProfile(){}; +}; + +/** + * Gaussian laser profile + */ +class GaussianLaserProfile : public ILaserProfile +{ + +public: + void + init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& ppc, + CommonLaserParameters params) override final; + + void + fill_amplitude ( + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) override final; + +private: + struct { + amrex::Real waist = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real duration = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real t_peak = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real focal_distance = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real zeta = 0; + amrex::Real beta = 0; + amrex::Real phi2 = 0; + + amrex::Vector<amrex::Real> stc_direction; //! Direction of the spatio-temporal couplings + amrex::Real theta_stc; //! Angle between polarization (p_X) and direction of spatiotemporal coupling (stc_direction) + } m_params; + + CommonLaserParameters m_common_params; +}; + +/** + * Harris laser profile + */ +class HarrisLaserProfile : public ILaserProfile +{ + +public: + void + init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& ppc, + CommonLaserParameters params) override final; + + void + fill_amplitude ( + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) override final; + +private: + struct { + amrex::Real waist = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real duration = std::numeric_limits<amrex::Real>::quiet_NaN(); + amrex::Real focal_distance = std::numeric_limits<amrex::Real>::quiet_NaN(); + } m_params; + + CommonLaserParameters m_common_params; +}; + +/** + * Laser profile defined by the used with an analytical expression + */ +class FieldFunctionLaserProfile : public ILaserProfile +{ + +public: + void + init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& ppc, + CommonLaserParameters params) override final; + + void + fill_amplitude ( + const int np, + amrex::Real const * AMREX_RESTRICT const Xp, + amrex::Real const * AMREX_RESTRICT const Yp, + amrex::Real t, + amrex::Real * AMREX_RESTRICT const amplitude) override final; + +private: + struct{ + std::string field_function; + } m_params; + + WarpXParser m_parser; +}; + +/** + * Maps laser profile names to lambdas returing unique pointers + * to the corresponding laser profile objects. + */ +const +std::map< +std::string, +std::function<std::unique_ptr<ILaserProfile>()> +> +laser_profiles_dictionary = +{ + {"gaussian", + [] () {return std::make_unique<GaussianLaserProfile>();} }, + {"harris", + [] () {return std::make_unique<HarrisLaserProfile>();} }, + {"parse_field_function", + [] () {return std::make_unique<FieldFunctionLaserProfile>();} } +}; + +} //WarpXLaserProfiles + +#endif //WARPX_LaserProfiles_H_ diff --git a/Source/Laser/LaserProfiles.cpp b/Source/Laser/LaserProfiles.cpp deleted file mode 100644 index 44411cedf..000000000 --- a/Source/Laser/LaserProfiles.cpp +++ /dev/null @@ -1,133 +0,0 @@ - -#include <WarpX_Complex.H> -#include <LaserParticleContainer.H> - -using namespace amrex; - -/* \brief compute field amplitude for a Gaussian laser, at particles' position - * - * Both Xp and Yp are given in laser plane coordinate. - * For each particle with position Xp and Yp, this routine computes the - * amplitude of the laser electric field, stored in array amplitude. - * - * \param np: number of laser particles - * \param Xp: pointer to first component of positions of laser particles - * \param Yp: pointer to second component of positions of laser particles - * \param t: Current physical time - * \param amplitude: pointer to array of field amplitude. - */ -void -LaserParticleContainer::gaussian_laser_profile ( - const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, - Real t, Real * AMREX_RESTRICT const amplitude) -{ - Complex I(0,1); - // Calculate a few factors which are independent of the macroparticle - const Real k0 = 2.*MathConst::pi/wavelength; - const Real inv_tau2 = 1. / (profile_duration * profile_duration); - const Real oscillation_phase = k0 * PhysConst::c * ( t - profile_t_peak ); - // The coefficients below contain info about Gouy phase, - // laser diffraction, and phase front curvature - const Complex diffract_factor = 1._rt + I * profile_focal_distance - * 2._rt/( k0 * profile_waist * profile_waist ); - const Complex inv_complex_waist_2 = 1._rt / ( profile_waist*profile_waist * diffract_factor ); - - // Time stretching due to STCs and phi2 complex envelope - // (1 if zeta=0, beta=0, phi2=0) - const Complex stretch_factor = 1._rt + 4._rt * - (zeta+beta*profile_focal_distance) * (zeta+beta*profile_focal_distance) - * (inv_tau2*inv_complex_waist_2) + - 2._rt *I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; - - // Amplitude and monochromatic oscillations - Complex prefactor = e_max * MathFunc::exp( I * oscillation_phase ); - - // Because diffract_factor is a complex, the code below takes into - // account the impact of the dimensionality on both the Gouy phase - // and the amplitude of the laser -#if (AMREX_SPACEDIM == 3) - prefactor = prefactor / diffract_factor; -#elif (AMREX_SPACEDIM == 2) - prefactor = prefactor / MathFunc::sqrt(diffract_factor); -#endif - - // Copy member variables to tmp copies for GPU runs. - Real tmp_profile_t_peak = profile_t_peak; - Real tmp_beta = beta; - Real tmp_zeta = zeta; - Real tmp_theta_stc = theta_stc; - Real tmp_profile_focal_distance = profile_focal_distance; - // Loop through the macroparticle to calculate the proper amplitude - amrex::ParallelFor( - np, - [=] AMREX_GPU_DEVICE (int i) { - const Complex stc_exponent = 1._rt / stretch_factor * inv_tau2 * - MathFunc::pow((t - tmp_profile_t_peak - - tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - - 2._rt *I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - *( tmp_zeta - tmp_beta*tmp_profile_focal_distance ) * inv_complex_waist_2),2); - // stcfactor = everything but complex transverse envelope - const Complex stcfactor = prefactor * MathFunc::exp( - stc_exponent ); - // Exp argument for transverse envelope - const Complex exp_argument = - ( Xp[i]*Xp[i] + Yp[i]*Yp[i] ) * inv_complex_waist_2; - // stcfactor + transverse envelope - amplitude[i] = ( stcfactor * MathFunc::exp( exp_argument ) ).real(); - } - ); -} - -/* \brief compute field amplitude for a Harris laser function, at particles' position - * - * Both Xp and Yp are given in laser plane coordinate. - * For each particle with position Xp and Yp, this routine computes the - * amplitude of the laser electric field, stored in array amplitude. - * - * \param np: number of laser particles - * \param Xp: pointer to first component of positions of laser particles - * \param Yp: pointer to second component of positions of laser particles - * \param t: Current physical time - * \param amplitude: pointer to array of field amplitude. - */ -void -LaserParticleContainer::harris_laser_profile ( - const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, - Real t, Real * AMREX_RESTRICT const amplitude) -{ - // This function uses the Harris function as the temporal profile of the pulse - const Real omega0 = 2.*MathConst::pi*PhysConst::c/wavelength; - const Real zR = MathConst::pi * profile_waist*profile_waist / wavelength; - const Real wz = profile_waist * - std::sqrt(1. + profile_focal_distance*profile_focal_distance/zR*zR); - const Real inv_wz_2 = 1./(wz*wz); - Real inv_Rz; - if (profile_focal_distance == 0.){ - inv_Rz = 0.; - } else { - inv_Rz = -profile_focal_distance / - ( profile_focal_distance*profile_focal_distance + zR*zR ); - } - const Real arg_env = 2.*MathConst::pi*t/profile_duration; - - // time envelope is given by the Harris function - Real time_envelope = 0.; - if (t < profile_duration) - time_envelope = 1./32. * (10. - 15.*std::cos(arg_env) + - 6.*std::cos(2.*arg_env) - - std::cos(3.*arg_env)); - - // Copy member variables to tmp copies for GPU runs. - Real tmp_e_max = e_max; - // Loop through the macroparticle to calculate the proper amplitude - amrex::ParallelFor( - np, - [=] AMREX_GPU_DEVICE (int i) { - const Real space_envelope = - std::exp(- ( Xp[i]*Xp[i] + Yp[i]*Yp[i] ) * inv_wz_2); - const Real arg_osc = omega0*t - omega0/PhysConst::c* - (Xp[i]*Xp[i] + Yp[i]*Yp[i]) * inv_Rz / 2.; - const Real oscillations = std::cos(arg_osc); - amplitude[i] = tmp_e_max * time_envelope * - space_envelope * oscillations; - } - ); -} diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp new file mode 100644 index 000000000..3c9d7379a --- /dev/null +++ b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp @@ -0,0 +1,45 @@ +#include <LaserProfiles.H> + +#include <WarpX_Complex.H> + +using namespace amrex; +using namespace WarpXLaserProfiles; + +void +FieldFunctionLaserProfile::init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& ppc, + CommonLaserParameters params) +{ + // Parse the properties of the parse_field_function profile + ppl.get("field_function(X,Y,t)", m_params.field_function); + m_parser.define(m_params.field_function); + m_parser.registerVariables({"X","Y","t"}); + + std::set<std::string> symbols = m_parser.symbols(); + symbols.erase("X"); + symbols.erase("Y"); + symbols.erase("t"); // after removing variables, we are left with constants + for (auto it = symbols.begin(); it != symbols.end(); ) { + Real v; + if (ppc.query(it->c_str(), v)) { + m_parser.setConstant(*it, v); + it = symbols.erase(it); + } else { + ++it; + } + } + for (auto const& s : symbols) { // make sure there no unknown symbols + amrex::Abort("Laser Profile: Unknown symbol "+s); + } +} + +void +FieldFunctionLaserProfile::fill_amplitude ( + const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, + Real t, Real * AMREX_RESTRICT const amplitude) +{ + for (int i = 0; i < np; ++i) { + amplitude[i] = m_parser.eval(Xp[i], Yp[i], t); + } +} diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp new file mode 100644 index 000000000..a0b5dd855 --- /dev/null +++ b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp @@ -0,0 +1,134 @@ +#include <LaserProfiles.H> + +#include <WarpX_Complex.H> +#include <WarpXConst.H> + +#include <cmath> + +using namespace amrex; +using namespace WarpXLaserProfiles; + +void +GaussianLaserProfile::init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& /* ppc */, + CommonLaserParameters params) +{ + //Copy common params + m_common_params = params; + + // Parse the properties of the Gaussian profile + ppl.get("profile_waist", m_params.waist); + ppl.get("profile_duration", m_params.duration); + ppl.get("profile_t_peak", m_params.t_peak); + ppl.get("profile_focal_distance", m_params.focal_distance); + ppl.query("zeta", m_params.zeta); + ppl.query("beta", m_params.beta); + ppl.query("phi2", m_params.phi2); + + m_params.stc_direction = m_common_params.p_X; + ppl.queryarr("stc_direction", m_params.stc_direction); + auto const s = 1.0_rt / std::sqrt( + m_params.stc_direction[0]*m_params.stc_direction[0] + + m_params.stc_direction[1]*m_params.stc_direction[1] + + m_params.stc_direction[2]*m_params.stc_direction[2]); + m_params.stc_direction = { + m_params.stc_direction[0]*s, + m_params.stc_direction[1]*s, + m_params.stc_direction[2]*s }; + auto const dp2 = + std::inner_product( + m_common_params.nvec.begin(), + m_common_params.nvec.end(), + m_params.stc_direction.begin(), 0.0); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp2) < 1.0e-14, + "stc_direction is not perpendicular to the laser plane vector"); + + // Get angle between p_X and stc_direction + // in 2d, stcs are in the simulation plane +#if AMREX_SPACEDIM == 3 + m_params.theta_stc = acos( + m_params.stc_direction[0]*m_common_params.p_X[0] + + m_params.stc_direction[1]*m_common_params.p_X[1] + + m_params.stc_direction[2]*m_common_params.p_X[2]); +#else + m_params.theta_stc = 0.; +#endif + +} + +/* \brief compute field amplitude for a Gaussian laser, at particles' position + * + * Both Xp and Yp are given in laser plane coordinate. + * For each particle with position Xp and Yp, this routine computes the + * amplitude of the laser electric field, stored in array amplitude. + * + * \param np: number of laser particles + * \param Xp: pointer to first component of positions of laser particles + * \param Yp: pointer to second component of positions of laser particles + * \param t: Current physical time + * \param amplitude: pointer to array of field amplitude. + */ +void +GaussianLaserProfile::fill_amplitude ( + const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, + Real t, Real * AMREX_RESTRICT const amplitude) +{ + Complex I(0,1); + // Calculate a few factors which are independent of the macroparticle + const Real k0 = 2.*MathConst::pi/m_common_params.wavelength; + const Real inv_tau2 = 1._rt /(m_params.duration * m_params.duration); + const Real oscillation_phase = k0 * PhysConst::c * ( t - m_params.t_peak ); + // The coefficients below contain info about Gouy phase, + // laser diffraction, and phase front curvature + const Complex diffract_factor = + 1._rt + I * m_params.focal_distance * 2._rt/ + ( k0 * m_params.waist * m_params.waist ); + const Complex inv_complex_waist_2 = + 1._rt /(m_params.waist*m_params.waist * diffract_factor ); + + // Time stretching due to STCs and phi2 complex envelope + // (1 if zeta=0, beta=0, phi2=0) + const Complex stretch_factor = 1._rt + 4._rt * + (m_params.zeta+m_params.beta*m_params.focal_distance) + * (m_params.zeta+m_params.beta*m_params.focal_distance) + * (inv_tau2*inv_complex_waist_2) + 2._rt *I * (m_params.phi2 + - m_params.beta*m_params.beta*k0*m_params.focal_distance) * inv_tau2; + + // Amplitude and monochromatic oscillations + Complex prefactor = + m_common_params.e_max * MathFunc::exp( I * oscillation_phase ); + + // Because diffract_factor is a complex, the code below takes into + // account the impact of the dimensionality on both the Gouy phase + // and the amplitude of the laser +#if (AMREX_SPACEDIM == 3) + prefactor = prefactor / diffract_factor; +#elif (AMREX_SPACEDIM == 2) + prefactor = prefactor / MathFunc::sqrt(diffract_factor); +#endif + + // Copy member variables to tmp copies for GPU runs. + auto const tmp_profile_t_peak = m_params.t_peak; + auto const tmp_beta = m_params.beta; + auto const tmp_zeta = m_params.zeta; + auto const tmp_theta_stc = m_params.theta_stc; + auto const tmp_profile_focal_distance = m_params.focal_distance; + // Loop through the macroparticle to calculate the proper amplitude + amrex::ParallelFor( + np, + [=] AMREX_GPU_DEVICE (int i) { + const Complex stc_exponent = 1._rt / stretch_factor * inv_tau2 * + MathFunc::pow((t - tmp_profile_t_peak - + tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - + 2._rt *I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) + *( tmp_zeta - tmp_beta*tmp_profile_focal_distance ) * inv_complex_waist_2),2); + // stcfactor = everything but complex transverse envelope + const Complex stcfactor = prefactor * MathFunc::exp( - stc_exponent ); + // Exp argument for transverse envelope + const Complex exp_argument = - ( Xp[i]*Xp[i] + Yp[i]*Yp[i] ) * inv_complex_waist_2; + // stcfactor + transverse envelope + amplitude[i] = ( stcfactor * MathFunc::exp( exp_argument ) ).real(); + } + ); +} diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp new file mode 100644 index 000000000..55374c5ea --- /dev/null +++ b/Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp @@ -0,0 +1,79 @@ +#include <LaserProfiles.H> + +#include <WarpX_Complex.H> +#include <WarpXConst.H> + +using namespace amrex; +using namespace WarpXLaserProfiles; + +void +HarrisLaserProfile::init ( + const amrex::ParmParse& ppl, + const amrex::ParmParse& /* ppc */, + CommonLaserParameters params) +{ + // Parse the properties of the Harris profile + ppl.get("profile_waist", m_params.waist); + ppl.get("profile_duration", m_params.duration); + ppl.get("profile_focal_distance", m_params.focal_distance); + //Copy common params + m_common_params = params; +} + +/* \brief compute field amplitude for a Harris laser function, at particles' position + * + * Both Xp and Yp are given in laser plane coordinate. + * For each particle with position Xp and Yp, this routine computes the + * amplitude of the laser electric field, stored in array amplitude. + * + * \param np: number of laser particles + * \param Xp: pointer to first component of positions of laser particles + * \param Yp: pointer to second component of positions of laser particles + * \param t: Current physical time + * \param amplitude: pointer to array of field amplitude. + */ +void +HarrisLaserProfile::fill_amplitude ( + const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, + Real t, Real * AMREX_RESTRICT const amplitude) +{ + // This function uses the Harris function as the temporal profile of the pulse + const Real omega0 = + 2._rt*MathConst::pi*PhysConst::c/m_common_params.wavelength; + const Real zR = MathConst::pi * m_params.waist*m_params.waist + / m_common_params.wavelength; + const Real wz = m_params.waist * + std::sqrt(1._rt + m_params.focal_distance*m_params.focal_distance/(zR*zR)); + const Real inv_wz_2 = 1._rt/(wz*wz); + Real inv_Rz; + if (m_params.focal_distance == 0.){ + inv_Rz = 0.; + } else { + inv_Rz = -m_params.focal_distance / + ( m_params.focal_distance*m_params.focal_distance + zR*zR ); + } + const Real arg_env = 2._rt*MathConst::pi*t/m_params.duration; + + // time envelope is given by the Harris function + Real time_envelope = 0.; + if (t < m_params.duration) + time_envelope = 1._rt/32._rt * (10._rt - 15._rt*std::cos(arg_env) + + 6._rt*std::cos(2._rt*arg_env) - + std::cos(3._rt*arg_env)); + + // Copy member variables to tmp copies for GPU runs. + const auto tmp_e_max = m_common_params.e_max; + // Loop through the macroparticle to calculate the proper amplitude + amrex::ParallelFor( + np, + [=] AMREX_GPU_DEVICE (int i) { + const Real space_envelope = + std::exp(- ( Xp[i]*Xp[i] + Yp[i]*Yp[i] ) * inv_wz_2); + const Real arg_osc = omega0*t - omega0/PhysConst::c* + (Xp[i]*Xp[i] + Yp[i]*Yp[i]) * inv_Rz / 2._rt; + const Real oscillations = std::cos(arg_osc); + amplitude[i] = tmp_e_max * time_envelope * + space_envelope * oscillations; + } + ); +} diff --git a/Source/Laser/LaserProfilesImpl/Make.package b/Source/Laser/LaserProfilesImpl/Make.package new file mode 100644 index 000000000..32284c4e4 --- /dev/null +++ b/Source/Laser/LaserProfilesImpl/Make.package @@ -0,0 +1,6 @@ +CEXE_sources += LaserProfileGaussian.cpp +CEXE_sources += LaserProfileHarris.cpp +CEXE_sources += LaserProfileFieldFunction.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Laser/LaserProfilesImpl +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Laser/LaserProfilesImpl diff --git a/Source/Laser/Make.package b/Source/Laser/Make.package index 0babbb8e4..f2de8797d 100644 --- a/Source/Laser/Make.package +++ b/Source/Laser/Make.package @@ -1,6 +1,9 @@ -CEXE_sources += LaserParticleContainer.cpp -CEXE_sources += LaserProfiles.cpp CEXE_headers += LaserParticleContainer.H +CEXE_sources += LaserParticleContainer.cpp +CEXE_headers += LaserProfiles.H +CEXE_headers += LaserAntennaParams.H + +include $(WARPX_HOME)/Source/Laser/LaserProfilesImpl/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Laser VPATH_LOCATIONS += $(WARPX_HOME)/Source/Laser diff --git a/Source/Parallelization/WarpXComm.H b/Source/Parallelization/WarpXComm.H new file mode 100644 index 000000000..81f00088e --- /dev/null +++ b/Source/Parallelization/WarpXComm.H @@ -0,0 +1,33 @@ +#ifndef WARPX_PARALLELIZATION_COMM_H_ +#define WARPX_PARALLELIZATION_COMM_H_ + +#include <AMReX_MultiFab.H> + +/** \brief Fills the values of the current on the coarse patch by + * averaging the values of the current of the fine patch (on the same level). + * Also fills the guards of the coarse patch. + * + * \param[in] fine fine patches to interpolate from + * \param[out] coarse coarse patches to interpolate to + * \param[in] refinement_ratio integer ratio between the two + */ +void +interpolateCurrentFineToCoarse ( + std::array< amrex::MultiFab const *, 3 > const & fine, + std::array< amrex::MultiFab *, 3 > const & coarse, + int const refinement_ratio); + +/** \brief Fills the values of the charge density on the coarse patch by + * averaging the values of the charge density of the fine patch (on the same level). + * + * \param[in] fine fine patches to interpolate from + * \param[out] coarse coarse patches to interpolate to + * \param[in] refinement_ratio integer ratio between the two + */ +void +interpolateDensityFineToCoarse ( + const amrex::MultiFab& fine, + amrex::MultiFab& coarse, + int const refinement_ratio); + +#endif // WARPX_PARALLELIZATION_COMM_H_ diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index f7e9086f9..3cc0eb563 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -1,3 +1,4 @@ +#include <WarpXComm.H> #include <WarpXComm_K.H> #include <WarpX.H> #include <WarpX_f.H> @@ -546,9 +547,9 @@ WarpX::SyncCurrent () } void -WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, - std::array< amrex::MultiFab *, 3 > const & coarse, - int const refinement_ratio) +interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, + std::array< amrex::MultiFab *, 3 > const & coarse, + int const refinement_ratio) { BL_PROFILE("interpolateCurrentFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); @@ -606,7 +607,7 @@ WarpX::SyncRho () } void -WarpX::interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) +interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio) { BL_PROFILE("interpolateDensityFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index a6c124ddd..6d34fa0ef 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -1,4 +1,3 @@ -F90EXE_sources += deposit_cic.F90 F90EXE_sources += interpolate_cic.F90 F90EXE_sources += push_particles_ES.F90 diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index fca22daa2..55768c6fc 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -321,11 +321,7 @@ void MultiParticleContainer::Redistribute () { for (auto& pc : allcontainers) { - if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { - pc->RedistributeCPU(); - } else { - pc->Redistribute(); - } + pc->Redistribute(); } } @@ -333,11 +329,7 @@ void MultiParticleContainer::RedistributeLocal (const int num_ghost) { for (auto& pc : allcontainers) { - if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { - pc->RedistributeCPU(0, 0, 0, num_ghost); - } else { - pc->Redistribute(0, 0, 0, num_ghost); - } + pc->Redistribute(0, 0, 0, num_ghost); } } diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 9132cf4a9..851726d57 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -41,9 +41,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; // Don't push momenta for photons diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 7e52b52e1..fd47ac8a0 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -51,9 +51,9 @@ void PhotonParticleContainer::InitData() void PhotonParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<ParticleReal>& xp, - Cuda::ManagedDeviceVector<ParticleReal>& yp, - Cuda::ManagedDeviceVector<ParticleReal>& zp, + Gpu::ManagedDeviceVector<ParticleReal>& xp, + Gpu::ManagedDeviceVector<ParticleReal>& yp, + Gpu::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 17a504719..174488eb9 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -94,9 +94,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full); virtual void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 51690d659..938c80de5 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -62,16 +62,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp "Radiation reaction can be enabled only if Boris pusher is used"); //_____________________________ -#ifdef AMREX_USE_GPU - Print()<<"\n-----------------------------------------------------\n"; - Print()<<"WARNING: field ionization on GPU uses RedistributeCPU\n"; - Print()<<"-----------------------------------------------------\n\n"; - //AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - //do_field_ionization == 0, - //"Field ionization does not work on GPU so far, because the current " - //"version of Redistribute in AMReX does not work with runtime parameters"); -#endif - #ifdef WARPX_QED //Add real component if QED is enabled pp.query("do_qed", m_do_qed); @@ -1406,7 +1396,7 @@ PhysicalParticleContainer::SplitParticles(int lev) { auto& mypc = WarpX::GetInstance().GetPartContainer(); auto& pctmp_split = mypc.GetPCtmp(); - Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; + Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; RealVector psplit_x, psplit_y, psplit_z, psplit_w; RealVector psplit_ux, psplit_uy, psplit_uz; long np_split_to_add = 0; @@ -1564,9 +1554,9 @@ PhysicalParticleContainer::SplitParticles(int lev) void PhysicalParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<ParticleReal>& xp, - Cuda::ManagedDeviceVector<ParticleReal>& yp, - Cuda::ManagedDeviceVector<ParticleReal>& zp, + Gpu::ManagedDeviceVector<ParticleReal>& xp, + Gpu::ManagedDeviceVector<ParticleReal>& yp, + Gpu::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index 3abbb4afe..a2473c5ad 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -44,9 +44,9 @@ public: DtType a_dt_type=DtType::Full) override; virtual void PushPX(WarpXParIter& pti, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& xp, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& yp, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& zp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& xp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& yp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; virtual void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 507206041..bee71fba1 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -78,7 +78,7 @@ RigidInjectedParticleContainer::RemapParticles() // Note that the particles are already in the boosted frame. // This value is saved to advance the particles not injected yet - Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; + Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -138,7 +138,7 @@ RigidInjectedParticleContainer::BoostandRemapParticles() #pragma omp parallel #endif { - Cuda::ManagedDeviceVector<ParticleReal> xp, yp, zp; + Gpu::ManagedDeviceVector<ParticleReal> xp, yp, zp; for (WarpXParIter pti(*this, 0); pti.isValid(); ++pti) { @@ -209,9 +209,9 @@ RigidInjectedParticleContainer::BoostandRemapParticles() void RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, - Cuda::ManagedDeviceVector<ParticleReal>& xp, - Cuda::ManagedDeviceVector<ParticleReal>& yp, - Cuda::ManagedDeviceVector<ParticleReal>& zp, + Gpu::ManagedDeviceVector<ParticleReal>& xp, + Gpu::ManagedDeviceVector<ParticleReal>& yp, + Gpu::ManagedDeviceVector<ParticleReal>& zp, Real dt, DtType a_dt_type) { @@ -222,7 +222,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& uzp = attribs[PIdx::uz]; // Save the position and momenta, making copies - Cuda::ManagedDeviceVector<ParticleReal> xp_save, yp_save, zp_save; + Gpu::ManagedDeviceVector<ParticleReal> xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; ParticleReal* const AMREX_RESTRICT x = xp.dataPtr(); diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index 28a1b73b7..80eaaf9cb 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -2,7 +2,6 @@ #define WARPX_PARTICLES_SORTING_SORTINGUTILS_H_ #include <WarpXParticleContainer.H> -#include <AMReX_CudaContainers.H> #include <AMReX_Gpu.H> #ifdef AMREX_USE_GPU #include <thrust/partition.h> @@ -43,7 +42,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin, // On GPU: Use thrust int const* AMREX_RESTRICT predicate_ptr = predicate.dataPtr(); ForwardIterator const sep = thrust::stable_partition( - thrust::cuda::par(amrex::Cuda::The_ThrustCachedAllocator()), + thrust::cuda::par(amrex::Gpu::The_ThrustCachedAllocator()), index_begin, index_end, [predicate_ptr] AMREX_GPU_DEVICE (long i) { return predicate_ptr[i]; } ); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index dbd913c5b..7f86930b2 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -73,12 +73,12 @@ public: WarpXParIter (ContainerType& pc, int level); #if (AMREX_SPACEDIM == 2) - void GetPosition (amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& x, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& y, - amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& z) const; - void SetPosition (const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& x, - const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& y, - const amrex::Cuda::ManagedDeviceVector<amrex::ParticleReal>& z); + void GetPosition (amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& x, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& y, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& z) const; + void SetPosition (const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& x, + const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& y, + const amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal>& z); #endif const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { return GetStructOfArrays().GetRealData(); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index d5520ba06..d22de00e0 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -4,6 +4,7 @@ #include <MultiParticleContainer.H> #include <WarpXParticleContainer.H> #include <AMReX_AmrParGDB.H> +#include <WarpXComm.H> #include <WarpX_f.H> #include <WarpX.H> #include <WarpXAlgorithmSelection.H> @@ -502,66 +503,54 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, void WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, bool local) { + // Loop over the refinement levels + int const finest_level = rho.size() - 1; + for (int lev = 0; lev < finest_level; ++lev) { - int num_levels = rho.size(); - int finest_level = num_levels - 1; - - // each level deposits it's own particles - const int ng = rho[0]->nGrow(); - for (int lev = 0; lev < num_levels; ++lev) { - - rho[lev]->setVal(0.0, ng); - - const auto& gm = m_gdb->Geom(lev); - const auto& ba = m_gdb->ParticleBoxArray(lev); - - const Real* dx = gm.CellSize(); - const Real* plo = gm.ProbLo(); - BoxArray nba = ba; - nba.surroundingNodes(); - - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = nba[pti]; - + // Loop over particle tiles and deposit charge on each level +#ifdef _OPENMP + #pragma omp parallel + { + int thread_num = omp_get_thread_num(); +#else + int thread_num = 0; +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const long np = pti.numParticles(); auto& wp = pti.GetAttribs(PIdx::w); - const auto& particles = pti.GetArrayOfStructs(); - int nstride = particles.dataShape().first; - const long np = pti.numParticles(); - FArrayBox& rhofab = (*rho[lev])[pti]; + pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - WRPX_DEPOSIT_CIC(particles.dataPtr(), nstride, np, - wp.dataPtr(), &this->charge, - rhofab.dataPtr(), box.loVect(), box.hiVect(), - plo, dx, &ng); - } + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } - if (!local) rho[lev]->SumBoundary(gm.periodicity()); + DepositCharge(pti, wp, ion_lev, rho[lev].get(), 0, 0, np, thread_num, lev, lev); + } +#ifdef _OPENMP + } +#endif + // Exchange guard cells + if (!local) rho[lev]->SumBoundary( m_gdb->Geom(lev).periodicity() ); } - // now we average down fine to crse - std::unique_ptr<MultiFab> crse; + // Now that the charge has been deposited at each level, + // we average down from fine to crse for (int lev = finest_level - 1; lev >= 0; --lev) { - const BoxArray& fine_BA = rho[lev+1]->boxArray(); const DistributionMapping& fine_dm = rho[lev+1]->DistributionMap(); - BoxArray coarsened_fine_BA = fine_BA; + BoxArray coarsened_fine_BA = rho[lev+1]->boxArray(); coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); coarsened_fine_data.setVal(0.0); - IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME + int const refinement_ratio = 2; - for (MFIter mfi(coarsened_fine_data); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); - const Box& crse_box = coarsened_fine_data[mfi].box(); - const Box& fine_box = (*rho[lev+1])[mfi].box(); - WRPX_SUM_FINE_TO_CRSE_NODAL(bx.loVect(), bx.hiVect(), ratio.getVect(), - coarsened_fine_data[mfi].dataPtr(), crse_box.loVect(), crse_box.hiVect(), - (*rho[lev+1])[mfi].dataPtr(), fine_box.loVect(), fine_box.hiVect()); - } - - rho[lev]->copy(coarsened_fine_data, m_gdb->Geom(lev).periodicity(), FabArrayBase::ADD); + interpolateDensityFineToCoarse( *rho[lev+1], coarsened_fine_data, refinement_ratio ); + rho[lev]->ParallelAdd( coarsened_fine_data, m_gdb->Geom(lev).periodicity() ); } } @@ -840,5 +829,3 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, if (pld.m_lev == lev-1){ } } - - diff --git a/Source/Particles/deposit_cic.F90 b/Source/Particles/deposit_cic.F90 deleted file mode 100644 index 2831ce96b..000000000 --- a/Source/Particles/deposit_cic.F90 +++ /dev/null @@ -1,134 +0,0 @@ -module warpx_ES_deposit_cic - - use iso_c_binding - use amrex_fort_module, only : amrex_real, amrex_particle_real - - implicit none - -contains - -! This routine computes the charge density due to the particles using cloud-in-cell -! deposition. The Fab rho is assumed to be node-centered. -! -! Arguments: -! particles : a pointer to the particle array-of-structs -! ns : the stride length of particle struct (the size of the struct in number of reals) -! np : the number of particles -! weights : the particle weights -! charge : the charge of this particle species -! rho : a Fab that will contain the charge density on exit -! lo : a pointer to the lo corner of this valid box for rho, in index space -! hi : a pointer to the hi corner of this valid box for rho, in index space -! plo : the real position of the left-hand corner of the problem domain -! dx : the mesh spacing -! ng : the number of ghost cells in rho -! - subroutine warpx_deposit_cic_3d(particles, ns, np, & - weights, charge, rho, lo, hi, plo, dx, & - ng) & - bind(c,name='warpx_deposit_cic_3d') - integer, value, intent(in) :: ns, np - real(amrex_particle_real), intent(in) :: particles(ns,np) - real(amrex_particle_real), intent(in) :: weights(np) - real(amrex_real), intent(in) :: charge - integer, intent(in) :: lo(3) - integer, intent(in) :: hi(3) - integer, intent(in) :: ng - real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng, lo(3)-ng:hi(3)+ng) - real(amrex_real), intent(in) :: plo(3) - real(amrex_real), intent(in) :: dx(3) - - integer i, j, k, n - real(amrex_real) wx_lo, wy_lo, wz_lo, wx_hi, wy_hi, wz_hi - real(amrex_real) lx, ly, lz - real(amrex_real) inv_dx(3) - real(amrex_real) qp, inv_vol - - inv_dx = 1.0d0/dx - - inv_vol = inv_dx(1) * inv_dx(2) * inv_dx(3) - - do n = 1, np - - qp = weights(n) * charge * inv_vol - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - lz = (particles(3, n) - plo(3))*inv_dx(3) - - i = floor(lx) - j = floor(ly) - k = floor(lz) - - wx_hi = lx - i - wy_hi = ly - j - wz_hi = lz - k - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - wz_lo = 1.0d0 - wz_hi - - rho(i, j, k ) = rho(i, j, k ) + wx_lo*wy_lo*wz_lo*qp - rho(i, j, k+1) = rho(i, j, k+1) + wx_lo*wy_lo*wz_hi*qp - rho(i, j+1, k ) = rho(i, j+1, k ) + wx_lo*wy_hi*wz_lo*qp - rho(i, j+1, k+1) = rho(i, j+1, k+1) + wx_lo*wy_hi*wz_hi*qp - rho(i+1, j, k ) = rho(i+1, j, k ) + wx_hi*wy_lo*wz_lo*qp - rho(i+1, j, k+1) = rho(i+1, j, k+1) + wx_hi*wy_lo*wz_hi*qp - rho(i+1, j+1, k ) = rho(i+1, j+1, k ) + wx_hi*wy_hi*wz_lo*qp - rho(i+1, j+1, k+1) = rho(i+1, j+1, k+1) + wx_hi*wy_hi*wz_hi*qp - - end do - - end subroutine warpx_deposit_cic_3d - - subroutine warpx_deposit_cic_2d(particles, ns, np, & - weights, charge, rho, lo, hi, plo, dx, & - ng) & - bind(c,name='warpx_deposit_cic_2d') - integer, value, intent(in) :: ns, np - real(amrex_particle_real), intent(in) :: particles(ns,np) - real(amrex_particle_real), intent(in) :: weights(np) - real(amrex_real), intent(in) :: charge - integer, intent(in) :: lo(2) - integer, intent(in) :: hi(2) - integer, intent(in) :: ng - real(amrex_real), intent(inout) :: rho(lo(1)-ng:hi(1)+ng, lo(2)-ng:hi(2)+ng) - real(amrex_real), intent(in) :: plo(2) - real(amrex_real), intent(in) :: dx(2) - - integer i, j, n - real(amrex_real) wx_lo, wy_lo, wx_hi, wy_hi - real(amrex_real) lx, ly - real(amrex_real) inv_dx(2) - real(amrex_real) qp, inv_vol - - inv_dx = 1.0d0/dx - - inv_vol = inv_dx(1) * inv_dx(2) - - do n = 1, np - - qp = weights(n) * charge * inv_vol - - lx = (particles(1, n) - plo(1))*inv_dx(1) - ly = (particles(2, n) - plo(2))*inv_dx(2) - - i = floor(lx) - j = floor(ly) - - wx_hi = lx - i - wy_hi = ly - j - - wx_lo = 1.0d0 - wx_hi - wy_lo = 1.0d0 - wy_hi - - rho(i, j ) = rho(i, j ) + wx_lo*wy_lo*qp - rho(i, j+1) = rho(i, j+1) + wx_lo*wy_hi*qp - rho(i+1, j ) = rho(i+1, j ) + wx_hi*wy_lo*qp - rho(i+1, j+1) = rho(i+1, j+1) + wx_hi*wy_hi*qp - - end do - - end subroutine warpx_deposit_cic_2d - -end module warpx_ES_deposit_cic diff --git a/Source/Utils/utils_ES.F90 b/Source/Utils/utils_ES.F90 index ce143bb94..baadeb7af 100644 --- a/Source/Utils/utils_ES.F90 +++ b/Source/Utils/utils_ES.F90 @@ -7,79 +7,6 @@ module warpx_ES_utils contains - subroutine warpx_sum_fine_to_crse_nodal_3d (lo, hi, lrat, crse, clo, chi, fine, flo, fhi) & - bind(c, name="warpx_sum_fine_to_crse_nodal_3d") - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: clo(3), chi(3) - integer, intent(in) :: flo(3), fhi(3) - integer, intent(in) :: lrat(3) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3)) - real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3)) - - integer :: i, j, k, ii, jj, kk - - do k = lo(3), hi(3) - kk = k * lrat(3) - do j = lo(2), hi(2) - jj = j * lrat(2) - do i = lo(1), hi(1) - ii = i * lrat(1) - crse(i,j,k) = fine(ii,jj,kk) + & -! These six fine nodes are shared by two coarse nodes... - 0.5d0 * (fine(ii-1,jj,kk) + fine(ii+1,jj,kk) + & - fine(ii,jj-1,kk) + fine(ii,jj+1,kk) + & - fine(ii,jj,kk-1) + fine(ii,jj,kk+1)) + & -! ... these twelve are shared by four... - 0.25d0 * (fine(ii,jj-1,kk-1) + fine(ii,jj+1,kk-1) + & - fine(ii,jj-1,kk+1) + fine(ii,jj+1,kk+1) + & - fine(ii-1,jj,kk-1) + fine(ii+1,jj,kk-1) + & - fine(ii-1,jj,kk+1) + fine(ii+1,jj,kk+1) + & - fine(ii-1,jj-1,kk) + fine(ii+1,jj-1,kk) + & - fine(ii-1,jj+1,kk) + fine(ii+1,jj+1,kk)) + & -! ... and these eight are shared by eight - 0.125d0 * (fine(ii-1,jj-1,kk-1) + fine(ii-1,jj-1,kk+1) + & - fine(ii-1,jj+1,kk-1) + fine(ii-1,jj+1,kk+1) + & - fine(ii+1,jj-1,kk-1) + fine(ii+1,jj-1,kk+1) + & - fine(ii+1,jj+1,kk-1) + fine(ii+1,jj+1,kk+1)) -! ... note that we have 27 nodes in total... - crse(i,j,k) = crse(i,j,k) / 8.d0 - end do - end do - end do - - end subroutine warpx_sum_fine_to_crse_nodal_3d - - subroutine warpx_sum_fine_to_crse_nodal_2d (lo, hi, lrat, crse, clo, chi, fine, flo, fhi) & - bind(c, name="warpx_sum_fine_to_crse_nodal_2d") - - integer, intent(in) :: lo(2), hi(2) - integer, intent(in) :: clo(2), chi(2) - integer, intent(in) :: flo(2), fhi(2) - integer, intent(in) :: lrat(2) - real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2)) - real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2)) - - integer :: i, j, ii, jj - - do j = lo(2), hi(2) - jj = j * lrat(2) - do i = lo(1), hi(1) - ii = i * lrat(1) - crse(i,j) = fine(ii,jj) + & -! These four fine nodes are shared by two coarse nodes... - 0.5d0 * (fine(ii-1,jj) + fine(ii+1,jj) + & - fine(ii,jj-1) + fine(ii,jj+1)) + & -! ... and these four are shared by four... - 0.25d0 * (fine(ii-1,jj-1) + fine(ii-1,jj+1) + & - fine(ii-1,jj+1) + fine(ii+1,jj+1)) -! ... note that we have 9 nodes in total... - crse(i,j) = crse(i,j) / 4.d0 - end do - end do - - end subroutine warpx_sum_fine_to_crse_nodal_2d - subroutine warpx_zero_out_bndry_3d (lo, hi, input_data, bndry_data, mask) & bind(c,name='warpx_zero_out_bndry_3d') diff --git a/Source/WarpX.H b/Source/WarpX.H index 9050196d6..04308ddee 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -393,9 +393,6 @@ private: void zeroOutBoundary(amrex::MultiFab& input_data, amrex::MultiFab& bndry_data, const amrex::FabArray<amrex::BaseFab<int> >& mask) const; - void sumFineToCrseNodal(const amrex::MultiFab& fine, amrex::MultiFab& crse, - const amrex::Geometry& cgeom, const amrex::IntVect& ratio); - void fixRHSForSolve(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhs, const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks) const ; @@ -441,29 +438,6 @@ private: std::array<std::unique_ptr<amrex::MultiFab>, 3> getInterpolatedB(int lev) const; - /** \brief Fills the values of the current on the coarse patch by - * averaging the values of the current of the fine patch (on the same level). - * Also fills the guards of the coarse patch. - * - * \param[in] fine fine patches to interpolate from - * \param[out] coarse coarse patches to interpolate to - * \param[in] refinement_ratio integer ratio between the two - */ - void interpolateCurrentFineToCoarse (std::array< amrex::MultiFab const *, 3 > const & fine, - std::array< amrex::MultiFab *, 3 > const & coarse, - int const refinement_ratio); - - /** \brief Fills the values of the charge density on the coarse patch by - * averaging the values of the charge density of the fine patch (on the same level). - * - * \param[in] fine fine patches to interpolate from - * \param[out] coarse coarse patches to interpolate to - * \param[in] refinement_ratio integer ratio between the two - */ - void interpolateDensityFineToCoarse (const amrex::MultiFab& fine, - amrex::MultiFab& coarse, - int const refinement_ratio); - void ExchangeWithPmlB (int lev); void ExchangeWithPmlE (int lev); void ExchangeWithPmlF (int lev); |