From 29c2a53bf13f8d1e93a4197ac8f914527336a193 Mon Sep 17 00:00:00 2001 From: "Eya D." <81635404+EyaDammak@users.noreply.github.com> Date: Fri, 17 Jan 2025 13:14:53 +0100 Subject: [PATCH] Laser ionization (#1190) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: AlexanderSinn Co-authored-by: Alexander Sinn <64009254+AlexanderSinn@users.noreply.github.com> Co-authored-by: Maxence Thévenet --- CMakeLists.txt | 6 + docs/source/run/parameters.rst | 3 + .../analysis_laser_ionization.py | 59 ++++ .../laser_ionization/inputs_laser_ionization | 62 ++++ src/Hipace.cpp | 5 +- src/laser/MultiLaser.H | 3 + src/particles/particles_utils/FieldGather.H | 74 +++++ src/particles/plasma/MultiPlasma.H | 9 + src/particles/plasma/MultiPlasma.cpp | 19 +- .../plasma/PlasmaParticleContainer.H | 37 ++- .../plasma/PlasmaParticleContainer.cpp | 278 ++++++++++++++++-- .../plasma/PlasmaParticleContainerInit.cpp | 30 +- .../laser_ionization.1Rank.json | 24 ++ tests/checksum/reset_all_benchmarks.sh | 12 + tests/laser_ionization.1Rank.sh | 48 +++ 15 files changed, 614 insertions(+), 55 deletions(-) create mode 100755 examples/laser_ionization/analysis_laser_ionization.py create mode 100644 examples/laser_ionization/inputs_laser_ionization create mode 100644 tests/checksum/benchmarks_json/laser_ionization.1Rank.json create mode 100755 tests/laser_ionization.1Rank.sh diff --git a/CMakeLists.txt b/CMakeLists.txt index dc4604f2f3..aa02e26781 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -342,6 +342,12 @@ if(BUILD_TESTING) WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} ) + add_test(NAME laser_ionization.1Rank + COMMAND bash ${HiPACE_SOURCE_DIR}/tests/laser_ionization.1Rank.sh + $ ${HiPACE_SOURCE_DIR} ${HiPACE_COMPUTE} + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} + ) + if (NOT HiPACE_COMPUTE STREQUAL CUDA) # These tests only run on CPU diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index df01af4681..1324b844aa 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -446,6 +446,9 @@ When both are specified, the per-species value is used. * ``.can_ionize`` (`bool`) optional (default `0`) Whether this plasma can ionize. Can also be set to 1 by specifying ``.ionization_product``. +* ``.can_laser_ionize`` (`bool`) optional (default `.can_ionize`) + Whether this plasma can be ionized by a laser. + * ``.initial_ion_level`` (`int`) optional (default `-1`) The initial ionization state of the plasma. `0` for neutral gasses. If set, the plasma charge gets multiplied by this number. If the plasma species is not ionizable, diff --git a/examples/laser_ionization/analysis_laser_ionization.py b/examples/laser_ionization/analysis_laser_ionization.py new file mode 100755 index 0000000000..46b562a946 --- /dev/null +++ b/examples/laser_ionization/analysis_laser_ionization.py @@ -0,0 +1,59 @@ +#! /usr/bin/env python3 + +# +# This file is part of HiPACE++. +# +# Authors: EyaDammak + +import argparse +import numpy as np +import math +from openpmd_viewer import OpenPMDTimeSeries +import statistics +from scipy.constants import e as qe +parser = argparse.ArgumentParser( + description='Script to analyze the equality of two simulations') +parser.add_argument('--first', + dest='first', + required=True, + help='Path to the directory containing output files') +parser.add_argument('--second', + dest='second', + required=True, + help='Path to the directory containing output files') +args = parser.parse_args() + +ts_linear = OpenPMDTimeSeries(args.first) +ts_circular = OpenPMDTimeSeries(args.second) + +a0_linear = 0.00885126 +a0_circular = 0.00787934 + +nc = 1.75e27 +n0 = nc / 10000 + +iteration = 0 + +rho_elec_linear, _ = ts_linear.get_field(field='rho_elec', coord='z', iteration=iteration) +rho_elec_mean_linear = np.mean(rho_elec_linear, axis=(1, 2)) +rho_average_linear = statistics.mean(rho_elec_mean_linear[0:10]) +fraction_linear = -rho_average_linear / qe / n0 + +rho_elec_circular, _ = ts_circular.get_field(field='rho_elec', coord='z', iteration=iteration) +rho_elec_mean_circular = np.mean(rho_elec_circular, axis=(1, 2)) +rho_average_circular = statistics.mean(rho_elec_mean_circular[0:10]) #average over a thickness in the ionized region +fraction_circular = -rho_average_circular / qe / n0 + +fraction_warpx_linear = 0.41014984 # result from WarpX simulation +fraction_warpx_circular = 0.502250841 # result from WarpX simulation + +relative_diff_linear = np.abs( ( fraction_linear - fraction_warpx_linear ) / fraction_warpx_linear ) +relative_diff_circular = np.abs( ( fraction_circular - fraction_warpx_circular ) / fraction_warpx_circular ) + +tolerance = 0.25 +print(f"fraction_warpx_linear = {fraction_warpx_linear}") +print(f"fraction_hipace_linear = {fraction_linear}") +print(f"fraction_warpx_circular = {fraction_warpx_circular}") +print(f"fraction_hipace_circular = {fraction_circular}") + +assert ( (relative_diff_linear < tolerance) and (relative_diff_circular < tolerance) ), 'Test laser_ionization did not pass' diff --git a/examples/laser_ionization/inputs_laser_ionization b/examples/laser_ionization/inputs_laser_ionization new file mode 100644 index 0000000000..74a49ed066 --- /dev/null +++ b/examples/laser_ionization/inputs_laser_ionization @@ -0,0 +1,62 @@ +max_step = 0 +hipace.dt = 0 +hipace.verbose = 3 + +amr.n_cell = 32 32 50 + + +my_constants.kp_inv = 10.e-6 + +my_constants.nc = 1.75e27 +my_constants.n0 = nc/10000 + +my_constants.a0 = 0.00885126 + +my_constants.w0_um = 1.e8 +my_constants.tau_fs = 30/1.17741 +my_constants.lambda0_nm = 800 + +hipace.file_prefix = laser_ionization.1Rank +hipace.do_tiling = 0 +hipace.deposit_rho_individual = 1 + +geometry.prob_lo = -10.*kp_inv -10.*kp_inv -15.*kp_inv #box shifted +geometry.prob_hi = 10.*kp_inv 10.*kp_inv 5.*kp_inv + +lasers.names = laser +lasers.lambda0 = lambda0_nm * 1.e-9 + + +laser.a0 = a0 +laser.position_mean = 0. 0. 0 +laser.w0 = w0_um*1.e-6 +laser.tau = tau_fs*1.e-15 +laser.focal_distance = 50.e-6 + +plasmas.names = elec ion + +elec.density(x,y,z) = n0 +elec.ppc = 0 0 +elec.u_mean = 0.0 0.0 0.0 +elec.element = electron +elec.neutralize_background = false + +ion.density(x,y,z) = n0 +ion.ppc = 1 1 +ion.u_mean = 0.0 0.0 0.0 +ion.element = H +ion.mass_Da = 1.008 +ion.initial_ion_level = 0 +ion.can_laser_ionize = 1 +ion.ionization_product = elec + +amr.max_level = 0 + +diagnostic.output_period = 1 + +hipace.depos_order_xy = 0 + +boundary.field = Dirichlet +boundary.particle = Periodic + +diagnostic.diag_type = xyz diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 3934ca30bb..20700c02b9 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -690,11 +690,14 @@ Hipace::SolveOneSlice (int islice, int step) // copy fields (and laser) to diagnostic array FillFieldDiagnostics(current_N_level, islice); - // plasma ionization + // plasma field ionization for (int lev=0; lev +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void doLaserGatherShapeN (const amrex::Real xp, + const amrex::Real yp, + amrex::GpuComplex& Ap, + amrex::GpuComplex& ADxp, + amrex::GpuComplex& ADzetap, + Array3 const& laser_arr, + const amrex::Real dx_inv, + const amrex::Real dy_inv, + const amrex::Real dzeta_inv, + const amrex::Real x_pos_offset, + const amrex::Real y_pos_offset) +{ + using namespace amrex::literals; + + // x,y direction + const amrex::Real x = (xp-x_pos_offset)*dx_inv; + const amrex::Real y = (yp-y_pos_offset)*dy_inv; + + // --- Compute shape factors + // x direction + // j_cell leftmost cell in x that the particle touches. sx_cell shape factor along x + amrex::Real sx_cell[depos_order_xy + 1]; + const int i_cell = compute_shape_factor(sx_cell, x); + + // y direction + amrex::Real sy_cell[depos_order_xy + 1]; + const int j_cell = compute_shape_factor(sy_cell, y); + + // Gather field Aabssq, AabssqDxp, and AabssqDxp on particle from field on grid slice_arr + // the derivative is calculated on the fly + for (int iy=0; iy<=depos_order_xy; iy++){ + for (int ix=0; ix<=depos_order_xy; ix++){ + using namespace WhichLaserSlice; + const amrex::Real s_cell = sx_cell[ix]*sy_cell[iy]; + const int i = i_cell+ix; + const int j = j_cell+iy; + Ap += amrex::GpuComplex { + s_cell * laser_arr(i, j, n00j00_r), + s_cell * laser_arr(i, j, n00j00_i) + }; + ADxp += amrex::GpuComplex { + s_cell * 0.5_rt * dx_inv * + (laser_arr(i + 1, j, n00j00_r) - laser_arr(i - 1, j, n00j00_r)), + s_cell * 0.5_rt * dx_inv * + (laser_arr(i + 1, j, n00j00_i) - laser_arr(i - 1, j, n00j00_i)) + }; + ADzetap += amrex::GpuComplex { + s_cell * dzeta_inv * + (laser_arr(i, j, n00jp1_r) - laser_arr(i, j, n00j00_r)), + s_cell * dzeta_inv * + (laser_arr(i, j, n00jp1_i) - laser_arr(i, j, n00j00_i)) + }; + } + } +} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void doGatherEz (const amrex::Real xp, const amrex::Real yp, diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index ea8de84748..e8e5843134 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -98,6 +98,15 @@ public: */ void DoFieldIonization (const int lev, const amrex::Geometry& geom, const Fields& fields); + /** Calculates Ionization Probability from laser and makes new Plasma Particles + * + * \param[in] islice zeta slice index + * \param[in] laser_geom Geometry of the laser + * \param[in] laser the laser class + */ + void DoLaserIonization (const int islice, const amrex::Geometry& laser_geom, + const MultiLaser& laser); + /** \brief whether any plasma species uses a neutralizing background, e.g. no ion motion */ bool AnySpeciesNeutralizeBackground () const; diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index 982fa4e063..3e4799a5f5 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -46,16 +46,16 @@ MultiPlasma::InitData (amrex::Vector slice_ba, plasma.InitData(gm[0]); if(plasma.m_can_ionize) { - PlasmaParticleContainer* plasma_product = nullptr; for (int i=0; i m_density_func; /**< Density function for the plasma */ amrex::Real m_min_density {0.}; /**< minimal density at which particles are injected */ @@ -176,11 +189,19 @@ public: amrex::Real m_temperature_in_ev {0.}; /**< Temperature of the plasma in eV */ /** whether to add a neutralizing background of immobile particles of opposite charge */ bool m_neutralize_background = true; + + // general: + amrex::Real m_mass = 0; /**< mass of each particle of this species */ amrex::Real m_charge = 0; /**< charge of each particle of this species, per Ion level */ - int m_init_ion_lev = -1; /**< initial Ion level of each particle */ int m_n_subcycles = 1; /**< number of subcycles in the plasma particle push */ - bool m_can_ionize = false; /**< whether this plasma can ionize */ + + // ionization: + + bool m_can_ionize = false; /**< whether this plasma can ionize from anything */ + bool m_can_field_ionize = false; /**< whether this plasma can ionize from the field */ + bool m_can_laser_ionize = false; /**< whether this plasma can ionize from a laser */ + int m_init_ion_lev = -1; /**< initial Ion level of each particle */ std::string m_product_name = ""; /**< name of Ionization product plasma */ PlasmaParticleContainer* m_product_pc = nullptr; /**< Ionization product plasma */ /** to calculate Ionization probability with ADK formula */ @@ -189,10 +210,18 @@ public: amrex::Gpu::DeviceVector m_adk_exp_prefactor; /** to calculate Ionization probability with ADK formula */ amrex::Gpu::DeviceVector m_adk_power; + /** to calculate laser Ionization probability with ADK formula */ + amrex::Gpu::DeviceVector m_laser_adk_prefactor; + + // plasma sorting/reordering: + /** After how many slices the particles are reordered. 0: off */ int m_reorder_period = 0; /** 2D reordering index type. 0: cell, 1: node, 2: both */ amrex::IntVect m_reorder_idx_type = {0, 0, 0}; + + // in-situ diagnostic: + /** How often the insitu plasma diagnostics should be computed and written * Default is 0, meaning no output */ int m_insitu_period {0}; diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 6c5e03f66a..31fa883973 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -3,7 +3,7 @@ * This file is part of HiPACE++. * * Authors: AlexanderSinn, Andrew Myers, MaxThevenet, Severin Diederichs - * Weiqun Zhang, Angel Ferran Pousa + * Weiqun Zhang, Angel Ferran Pousa, EyaDammak * License: BSD-3-Clause-LBNL */ #include "Hipace.H" @@ -65,9 +65,14 @@ PlasmaParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_mass != 0, "The plasma particle mass must be specified"); bool ion_lev_specified = queryWithParser(pp, "initial_ion_level", m_init_ion_lev); - m_can_ionize = pp.contains("ionization_product"); + m_can_field_ionize = pp.contains("ionization_product"); + + queryWithParser(pp, "can_ionize", m_can_field_ionize); + m_can_laser_ionize = false; + queryWithParser(pp, "can_laser_ionize", m_can_laser_ionize); + + m_can_ionize = m_can_field_ionize || m_can_laser_ionize; - queryWithParser(pp, "can_ionize", m_can_ionize); if(m_can_ionize) { m_neutralize_background = false; // change default AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_init_ion_lev >= 0, @@ -265,7 +270,7 @@ IonizationModule (const int lev, const Fields& fields, const amrex::Real background_density_SI) { - if (!m_can_ionize) return; + if (!m_can_field_ionize) return; HIPACE_PROFILE("PlasmaParticleContainer::IonizationModule()"); using namespace amrex::literals; @@ -302,7 +307,7 @@ IonizationModule (const int lev, auto& soa_ion = ptile_ion.GetStructOfArrays(); // For momenta and weights const amrex::Real clightsq = 1.0_rt / ( phys_const.c * phys_const.c ); - // calcuation of E0 in SI units for denormalization + // Calculation of E0 in SI units for denormalization const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * PhysConstSI::q_e*PhysConstSI::q_e / (PhysConstSI::ep0 * PhysConstSI::m_e) ); @@ -329,6 +334,11 @@ IonizationModule (const int lev, long num_ions = ptile_ion.numParticles(); + + // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time + // and calculates ionization probability. If ionization occurs, it increments + // `p_num_new_electrons` to calculate the number of ionized electrons. + // It also constructs a mask with 1 boolean per macro-ion: 1 if ionized, 0 otherwise. amrex::AnyCTO( amrex::TypeList< amrex::CompileTimeOptions<0, 1, 2, 3> @@ -336,7 +346,8 @@ IonizationModule (const int lev, Hipace::m_depos_order_xy }, [&] (auto cto_func) { - amrex::ParallelForRNG(num_ions, cto_func); + amrex::ParallelForRNG(num_ions, cto_func); // enables the use of amrex::Random within the loop + }, [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine, auto depos_order_xy) { @@ -344,11 +355,11 @@ IonizationModule (const int lev, if (amrex::ConstParticleIDWrapper(idcpup[ip]) < 0 || amrex::ConstParticleCPUWrapper(idcpup[ip]) != lev) return; - // avoid temp slice + // Avoid temp slice const amrex::Real xp = x_prev[ip]; const amrex::Real yp = y_prev[ip]; - // define field at particle position reals + // Define field at particle position reals amrex::ParticleReal ExmByp = 0., EypBxp = 0., Ezp = 0.; amrex::ParticleReal Bxp = 0., Byp = 0., Bzp = 0.; @@ -377,7 +388,7 @@ IonizationModule (const int lev, { ion_lev[ip] += 1; p_ion_mask[ip] = 1; - amrex::Gpu::Atomic::Add( p_num_new_electrons, 1u ); + amrex::Gpu::Atomic::Add( p_num_new_electrons, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` } }); amrex::Gpu::streamSynchronize(); @@ -385,12 +396,12 @@ IonizationModule (const int lev, if (num_new_electrons.dataValue() == 0) continue; if(Hipace::m_verbose >= 3) { - amrex::Print() << "Number of ionized Plasma Particles: " + amrex::Print() << "Number of ionized Plasma Particles (field): " << num_new_electrons.dataValue() << "\n"; } - // resize electron particle tile + // Resize electron particle tile const auto old_size = ptile_elec.numParticles(); const auto new_size = old_size + num_new_electrons.dataValue(); ptile_elec.resize(new_size); @@ -406,15 +417,16 @@ IonizationModule (const int lev, amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); + // This kernel adds the new ionized electrons to the Plasma Particle Container amrex::ParallelFor(num_ions, [=] AMREX_GPU_DEVICE (long ip) { if(p_ion_mask[ip] != 0) { - const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); + const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` const long pidx = pid + old_size; // Copy ion data to new electron - amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // only for valid/invalid + amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // sets the ionized electron ID to 2 (valid/invalid) for the new electron amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = lev; // current level arrdata_elec[PlasmaIdx::x ][pidx] = arrdata_ion[PlasmaIdx::x ][ip]; arrdata_elec[PlasmaIdx::y ][pidx] = arrdata_ion[PlasmaIdx::y ][ip]; @@ -422,7 +434,213 @@ IonizationModule (const int lev, arrdata_elec[PlasmaIdx::w ][pidx] = arrdata_ion[PlasmaIdx::w ][ip]; arrdata_elec[PlasmaIdx::ux ][pidx] = 0._rt; arrdata_elec[PlasmaIdx::uy ][pidx] = 0._rt; - // later we could consider adding a finite temperature to the ionized electrons + // Later we could consider adding a finite temperature to the ionized electrons + arrdata_elec[PlasmaIdx::psi ][pidx] = 1._rt; + arrdata_elec[PlasmaIdx::x_prev ][pidx] = arrdata_ion[PlasmaIdx::x_prev][ip]; + arrdata_elec[PlasmaIdx::y_prev ][pidx] = arrdata_ion[PlasmaIdx::y_prev][ip]; + arrdata_elec[PlasmaIdx::ux_half_step ][pidx] = 0._rt; + arrdata_elec[PlasmaIdx::uy_half_step ][pidx] = 0._rt; + arrdata_elec[PlasmaIdx::psi_half_step][pidx] = 1._rt; +#ifdef HIPACE_USE_AB5_PUSH +#ifdef AMREX_USE_GPU +#pragma unroll +#endif + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + arrdata_elec[iforce][pidx] = 0._rt; + } +#endif + int_arrdata_elec[PlasmaIdx::ion_lev][pidx] = init_ion_lev; + } + }); + + // Synchronize before ion_mask and ip_elec go out of scope + amrex::Gpu::streamSynchronize(); + } +} + +void +PlasmaParticleContainer:: +LaserIonization (const int islice, + const amrex::Geometry& laser_geom, + const MultiLaser& laser, + const amrex::Real background_density_SI) +{ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( !m_can_laser_ionize || laser.UseLaser(), + "Error: LaserIonization requires the laser to be enabled in the current slice."); + if (!m_can_laser_ionize || !laser.UseLaser(islice)) return; + HIPACE_PROFILE("PlasmaParticleContainer::LaserIonization()"); + + using namespace amrex::literals; + using Complex = amrex::GpuComplex; + constexpr Complex I(0.,1.); + + const PhysConst phys_const = get_phys_const(); + + // Loop over particle boxes with both ion and electron Particle Containers at the same time + for (amrex::MFIter mfi_ion = MakeMFIter(0, DfltMfi); mfi_ion.isValid(); ++mfi_ion) + { + // Extract laser array + Array3 const laser_arr = laser.getSlices().const_array(mfi_ion); + + const CheckDomainBounds laser_bounds {laser_geom}; + + // Extract properties associated with physical size of the box + const amrex::Real dx_inv = laser_geom.InvCellSize(0); + const amrex::Real dy_inv = laser_geom.InvCellSize(1); + const amrex::Real dzeta_inv = laser_geom.InvCellSize(2); + + // Offset for converting positions to indexes + amrex::Real const x_pos_offset = GetPosOffset(0, laser_geom, laser_geom.Domain()); + amrex::Real const y_pos_offset = GetPosOffset(1, laser_geom, laser_geom.Domain()); + + auto& plevel_ion = GetParticles(0); + auto index = std::make_pair(mfi_ion.index(), mfi_ion.LocalTileIndex()); + if(plevel_ion.find(index) == plevel_ion.end()) continue; + auto& ptile_elec = m_product_pc->DefineAndReturnParticleTile(0, + mfi_ion.index(), mfi_ion.LocalTileIndex()); + auto& ptile_ion = plevel_ion.at(index); + + auto& soa_ion = ptile_ion.GetStructOfArrays(); // for momenta and weights + + const amrex::Real clightsq = 1.0_rt / ( phys_const.c * phys_const.c ); + // Calcuation of E0 in SI units for denormalization + const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * + PhysConstSI::q_e*PhysConstSI::q_e / + (PhysConstSI::ep0 * PhysConstSI::m_e) ); + const amrex::Real E0 = Hipace::m_normalized_units ? + wp * PhysConstSI::m_e * PhysConstSI::c / PhysConstSI::q_e : 1; + const amrex::Real lambda0 = laser.GetLambda0(); + const amrex::Real omega0 = 2.0 * MathConst::pi * phys_const.c / lambda0; + const bool linear_polarization = laser.LinearPolarization(); + + int * const ion_lev = soa_ion.GetIntData(PlasmaIdx::ion_lev).data(); + const amrex::Real * const x_prev = soa_ion.GetRealData(PlasmaIdx::x_prev).data(); + const amrex::Real * const y_prev = soa_ion.GetRealData(PlasmaIdx::y_prev).data(); + const amrex::Real * const uxp = soa_ion.GetRealData(PlasmaIdx::ux_half_step).data(); + const amrex::Real * const uyp = soa_ion.GetRealData(PlasmaIdx::uy_half_step).data(); + const amrex::Real * const psip =soa_ion.GetRealData(PlasmaIdx::psi_half_step).data(); + const auto * idcpup = soa_ion.GetIdCPUData().data(); + + // Make Ion Mask and load ADK prefactors + // Ion Mask is necessary to only resize electron particle tile once + amrex::Gpu::DeviceVector ion_mask(ptile_ion.numParticles(), 0); + uint8_t* AMREX_RESTRICT p_ion_mask = ion_mask.data(); + amrex::Gpu::DeviceScalar num_new_electrons(0); + uint32_t* AMREX_RESTRICT p_num_new_electrons = num_new_electrons.dataPtr(); + amrex::Real* AMREX_RESTRICT adk_prefactor = m_adk_prefactor.data(); + amrex::Real* AMREX_RESTRICT adk_exp_prefactor = m_adk_exp_prefactor.data(); + amrex::Real* AMREX_RESTRICT adk_power = m_adk_power.data(); + amrex::Real* AMREX_RESTRICT laser_adk_prefactor = m_laser_adk_prefactor.data(); + + long num_ions = ptile_ion.numParticles(); + + // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time + // and calculates ionization probability. If ionization occurs, it increments + // `p_num_new_electrons` to calculate the number of ionized electrons. + // It also constructs a mask with 1 boolean per macro-ion: 1 if ionized, 0 otherwise. + amrex::AnyCTO( + amrex::TypeList< + amrex::CompileTimeOptions<0, 1, 2, 3> + >{}, { + Hipace::m_depos_order_xy + }, + [&] (auto cto_func) { + amrex::ParallelForRNG(num_ions, cto_func); // enables the use of `amrex::Random` within the loop + + }, + [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine, + auto depos_order_xy) { + + // Avoid temp slice + const amrex::Real xp = x_prev[ip]; + const amrex::Real yp = y_prev[ip]; + + if (amrex::ConstParticleIDWrapper(idcpup[ip]) < 0 || + !laser_bounds.contains(xp, yp)) return; + + Complex A = 0; + Complex A_dx = 0; + Complex A_dzeta = 0; + + doLaserGatherShapeN(xp, yp, A, A_dx, A_dzeta, laser_arr, + dx_inv, dy_inv, dzeta_inv, x_pos_offset, y_pos_offset); + + // Convert from vector potential to electric field. Units are fixed later. + const Complex Et = I * A * omega0 + A_dzeta * phys_const.c; // transverse component + const Complex El = - A_dx * phys_const.c; // longitudinal component + + // Get amplitude of the electric field envelope and normalize to correct SI unit. + amrex::Real Ep = std::sqrt( amrex::abs(Et*Et) + amrex::abs(El*El) ); + Ep *= phys_const.m_e * phys_const.c / phys_const.q_e * E0; + + // Compute probability of ionization p + const amrex::Real gammap = (1.0_rt + uxp[ip] * uxp[ip] * clightsq + + uyp[ip] * uyp[ip] * clightsq + + psip[ip]* psip[ip] ) / ( 2.0_rt * psip[ip] ); + const int ion_lev_loc = ion_lev[ip]; + // gamma / (psi + 1) to complete dt for QSA + amrex::Real w_dtau_dc = gammap / psip[ip] * adk_prefactor[ion_lev_loc] * + std::pow(Ep, adk_power[ion_lev_loc]) * + std::exp( adk_exp_prefactor[ion_lev_loc]/Ep ); + + amrex::Real const w_dtau_ac = w_dtau_dc * + (linear_polarization ? std::sqrt(Ep * laser_adk_prefactor[ion_lev_loc]) : 1._rt); + + amrex::Real p = 1._rt - std::exp( - w_dtau_ac ); + + amrex::Real random_draw = amrex::Random(engine); + if (random_draw < p) + { + ion_lev[ip] += 1; + p_ion_mask[ip] = 1; + amrex::Gpu::Atomic::Add( p_num_new_electrons, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` + } + }); + amrex::Gpu::streamSynchronize(); + + if (num_new_electrons.dataValue() == 0) continue; + + if(Hipace::m_verbose >= 3) { + amrex::Print() << "Number of ionized Plasma Particles (laser): " + << num_new_electrons.dataValue() << "\n"; + } + + + // Resize electron particle tile + const auto old_size = ptile_elec.numParticles(); + const auto new_size = old_size + num_new_electrons.dataValue(); + ptile_elec.resize(new_size); + + // Load electron soa and aos after resize + auto arrdata_ion = ptile_ion.GetStructOfArrays().realarray(); + auto arrdata_elec = ptile_elec.GetStructOfArrays().realarray(); + auto int_arrdata_elec = ptile_elec.GetStructOfArrays().intarray(); + auto idcpu_elec = ptile_elec.GetStructOfArrays().GetIdCPUData().data(); + auto idcpu_ion = ptile_ion.GetStructOfArrays().GetIdCPUData().data(); + + const int init_ion_lev = m_product_pc->m_init_ion_lev; + + amrex::Gpu::DeviceScalar ip_elec(0); + uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); + + // This kernel adds the new ionized electrons to the Plasma Particle Container + amrex::ParallelFor(num_ions, + [=] AMREX_GPU_DEVICE (long ip) { + + if(p_ion_mask[ip] != 0) { + const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` + const long pidx = pid + old_size; + + // Copy ion data to new electron + amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2; // sets the ionized electron ID to 2 (valid/invalid) for the new electron + amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = + amrex::ParticleCPUWrapper{idcpu_ion[pidx]}; // current level + arrdata_elec[PlasmaIdx::x ][pidx] = arrdata_ion[PlasmaIdx::x ][ip]; + arrdata_elec[PlasmaIdx::y ][pidx] = arrdata_ion[PlasmaIdx::y ][ip]; + + arrdata_elec[PlasmaIdx::w ][pidx] = arrdata_ion[PlasmaIdx::w ][ip]; + arrdata_elec[PlasmaIdx::ux ][pidx] = 0._rt; + arrdata_elec[PlasmaIdx::uy ][pidx] = 0._rt; arrdata_elec[PlasmaIdx::psi ][pidx] = 1._rt; arrdata_elec[PlasmaIdx::x_prev ][pidx] = arrdata_ion[PlasmaIdx::x_prev][ip]; arrdata_elec[PlasmaIdx::y_prev ][pidx] = arrdata_ion[PlasmaIdx::y_prev][ip]; @@ -441,7 +659,7 @@ IonizationModule (const int lev, } }); - // synchronize before ion_mask and ip_elec go out of scope + // Synchronize before ion_mask and ip_elec go out of scope amrex::Gpu::streamSynchronize(); } } @@ -463,7 +681,7 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) // Loop over particle boxes for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti) { - // loading the data + // Loading the data const auto ptd = pti.GetParticleTile().getParticleTileData(); amrex::Long const num_particles = pti.numParticles(); @@ -484,13 +702,13 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) if (!ptd.id(ip).is_valid() || x*x + y*y > insitu_radius_sq) { return amrex::IdentityTuple(ReduceTuple{}, reduce_op); } - // particle's Lorentz factor + // Particle's Lorentz factor const amrex::Real gamma = (1.0_rt + ux*ux + uy*uy + psi*psi)/(2.0_rt*psi); - // the *c from uz cancels with the /c from the proper velocity conversion + // The *c from uz cancels with the /c from the proper velocity conversion const amrex::Real uz = (gamma - psi); - // weight with quasi-static weighting factor + // Weight with quasi-static weighting factor const amrex::Real w = ptd.rdata(PlasmaIdx::w)[ip] * gamma/psi; - // no quasi-static weighting factor to calculate quasi-static energy + // No quasi-static weighting factor to calculate quasi-static energy const amrex::Real energy = ptd.rdata(PlasmaIdx::w)[ip] * (gamma - 1._rt); return { // Tuple contains: w, // 0 sum(w) @@ -539,16 +757,16 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am HIPACE_PROFILE("PlasmaParticleContainer::InSituWriteToFile()"); #ifdef HIPACE_USE_OPENPMD - // create subdirectory + // Create subdirectory openPMD::auxiliary::create_directories(m_insitu_file_prefix); #endif - // zero pad the rank number; + // Zero pad the rank number; std::string::size_type n_zeros = 4; std::string rank_num = std::to_string(amrex::ParallelDescriptor::MyProc()); std::string pad_rank_num = std::string(n_zeros-std::min(rank_num.size(), n_zeros),'0')+rank_num; - // open file + // Open file std::ofstream ofs{m_insitu_file_prefix + "/reduced_" + m_name + "." + pad_rank_num + ".txt", std::ofstream::out | std::ofstream::app | std::ofstream::binary}; @@ -558,8 +776,8 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am geom.CellSizeArray().product() : 1; // dx * dy * dz in normalized units, 1 otherwise const int is_normalized_units = Hipace::m_normalized_units; - // specify the structure of the data later available in python - // avoid pointers to temporary objects as second argument, stack variables are ok + // Specify the structure of the data later available in python + // Avoid pointers to temporary objects as second argument, stack variables are ok const amrex::Vector all_data{ {"time" , &time}, {"step" , &step}, @@ -607,16 +825,16 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am }; if (ofs.tellp() == 0) { - // write JSON header containing a NumPy structured datatype + // Write JSON header containing a NumPy structured datatype insitu_utils::write_header(all_data, ofs); } - // write binary data according to datatype in header + // Write binary data according to datatype in header insitu_utils::write_data(all_data, ofs); - // close file + // Close file ofs.close(); - // assert no file errors + // Assert no file errors #ifdef HIPACE_USE_OPENPMD AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ofs, "Error while writing insitu plasma diagnostics"); #else @@ -624,7 +842,7 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am "Maybe the specified subdirectory does not exist"); #endif - // reset arrays for insitu data + // Reset arrays for insitu data for (auto& x : m_insitu_rdata) x = 0.; for (auto& x : m_insitu_idata) x = 0; for (auto& x : m_insitu_sum_rdata) x = 0.; diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index 9744f7fd47..189247bb5e 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -228,7 +228,7 @@ InitParticles (const amrex::RealVect& a_u_std, unsigned int uiy = amrex::min(ny-1,amrex::max(0,iy)); unsigned int uiz = amrex::min(nz-1,amrex::max(0,iz)); - // ordering of axes from fastest to slowest: + // Ordering of axes from fastest to slowest: // x // y // z (not used) @@ -379,15 +379,12 @@ InitParticles (const amrex::RealVect& a_u_std, void PlasmaParticleContainer:: -InitIonizationModule (const amrex::Geometry& geom, PlasmaParticleContainer* product_pc, - const amrex::Real background_density_SI) +InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_density_SI) { HIPACE_PROFILE("PlasmaParticleContainer::InitIonizationModule()"); using namespace amrex::literals; - if (!m_can_ionize) return; - const bool normalized_units = Hipace::m_normalized_units; if (normalized_units) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(background_density_SI!=0, @@ -395,14 +392,13 @@ InitIonizationModule (const amrex::Geometry& geom, PlasmaParticleContainer* prod "be specified via 'hipace.background_density_SI'"); } - m_product_pc = product_pc; amrex::ParmParse pp(m_name); std::string physical_element; getWithParser(pp, "element", physical_element); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ion_map_ids.count(physical_element) != 0, "There are no ionization energies available for this element. " "Please update src/utils/IonizationEnergiesTable.H using write_atomic_data_cpp.py"); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE((std::abs(product_pc->m_charge / m_charge +1) < 1e-3), + AMREX_ALWAYS_ASSERT_WITH_MESSAGE((std::abs(m_product_pc->m_charge / m_charge +1) < 1e-3), "Ion and Ionization product charges have to be opposite"); // Get atomic number and ionization energies from file const int ion_element_id = ion_map_ids[physical_element]; @@ -418,15 +414,14 @@ InitIonizationModule (const amrex::Geometry& geom, PlasmaParticleContainer* prod // without Gamma function const PhysConst phys_const = make_constants_SI(); const amrex::Real alpha = 0.0072973525693_rt; - const amrex::Real r_e = 2.8179403227e-15_rt; const amrex::Real a3 = alpha * alpha * alpha; const amrex::Real a4 = a3 * alpha; - const amrex::Real wa = a3 * phys_const.c / r_e; - const amrex::Real Ea = phys_const.m_e * phys_const.c * phys_const.c / phys_const.q_e * a4 / r_e; + const amrex::Real wa = a3 * phys_const.c / PhysConstSI::r_e; + const amrex::Real Ea = phys_const.m_e * phys_const.c * phys_const.c / phys_const.q_e * a4 / PhysConstSI::r_e; const amrex::Real UH = table_ionization_energies[0]; const amrex::Real l_eff = std::sqrt(UH/h_ionization_energies[0]) - 1._rt; - // plasma frequency in SI units to denormalize ionization + // Plasma frequency in SI units to denormalize ionization const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * PhysConstSI::q_e*PhysConstSI::q_e / (PhysConstSI::ep0 * PhysConstSI::m_e) ); @@ -435,21 +430,24 @@ InitIonizationModule (const amrex::Geometry& geom, PlasmaParticleContainer* prod m_adk_power.resize(ion_atomic_number); m_adk_prefactor.resize(ion_atomic_number); m_adk_exp_prefactor.resize(ion_atomic_number); + m_laser_adk_prefactor.resize(ion_atomic_number); amrex::Gpu::PinnedVector h_adk_power(ion_atomic_number); amrex::Gpu::PinnedVector h_adk_prefactor(ion_atomic_number); amrex::Gpu::PinnedVector h_adk_exp_prefactor(ion_atomic_number); + amrex::Gpu::PinnedVector h_laser_adk_prefactor(ion_atomic_number); for (int i=0; i