diff options
| -rw-r--r-- | Source/Laser/LaserParticleContainer.H | 40 | ||||
| -rw-r--r-- | Source/Laser/LaserParticleContainer.cpp | 102 | ||||
| -rw-r--r-- | Source/Laser/LaserProfiles.H | 202 | ||||
| -rw-r--r-- | Source/Laser/LaserProfiles.cpp | 133 | ||||
| -rw-r--r-- | Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp | 45 | ||||
| -rw-r--r-- | Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp | 134 | ||||
| -rw-r--r-- | Source/Laser/LaserProfilesImpl/LaserProfileHarris.cpp | 79 | ||||
| -rw-r--r-- | Source/Laser/LaserProfilesImpl/Make.package | 6 | ||||
| -rw-r--r-- | Source/Laser/Make.package | 7 |
9 files changed, 500 insertions, 248 deletions
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 35eadf064..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. @@ -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 |