diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index d3bd9555b5..28cd3d91d5 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -1165,6 +1165,11 @@ Field diagnostics will be deposited into individual fields accessible as ``w``, ``ux_`` or ``ux^2_`` (similarly for ``uy`` and ``uz``) in ``diagnostic.field_data``. +* ``hipace.temperature_depos_order`` (`int`) optional (default `2`) + When ``hipace.deposit_temp_individual`` is turned on, + this option specifies the shape order of the deposited fields. + Currently, 0,1,2,3 are implemented. + In-situ diagnostics ^^^^^^^^^^^^^^^^^^^ diff --git a/src/Hipace.H b/src/Hipace.H index 240e357775..0a6c499148 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -242,6 +242,8 @@ public: inline static bool m_deposit_rho_individual = false; /** Whether to deposit temperature for every individual plasma for diagnostics */ inline static bool m_deposit_temp_individual = false; + /** Order of the temperature deposition shape factor */ + inline static int m_temperature_depos_order = 2; /** Whether to interpolate the neutralizing background to MR levels 1 and 2 instead of depositing it */ inline static bool m_interpolate_neutralizing_background = false; /** Whether to use tiling for particle operations */ diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 5d551e7dab..c7ffaa63e3 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -130,6 +130,7 @@ Hipace::ReadParameters () queryWithParser(pph, "deposit_rho_individual", m_deposit_rho_individual); m_deposit_temp_individual = m_diags.needsTempIndividual(); queryWithParser(pph, "deposit_temp_individual", m_deposit_temp_individual); + queryWithParser(pph, "temperature_depos_order", m_temperature_depos_order); queryWithParser(pph, "interpolate_neutralizing_background", m_interpolate_neutralizing_background); bool do_mfi_sync = false; diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index ae9cc83977..78f69dbc51 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -63,6 +63,11 @@ Fields::AllocData ( // Need 1 extra guard cell transversally for transverse derivative int nguards_xy = (Hipace::m_depos_order_xy + 1) / 2 + 1; + // Check the temperature deposition order, if enabled + if (Hipace::m_deposit_temp_individual && + Hipace::m_temperature_depos_order > Hipace::m_depos_order_xy) { + nguards_xy = (Hipace::m_temperature_depos_order + 1) / 2 + 1; + } m_slices_nguards = amrex::IntVect{nguards_xy, nguards_xy, 0}; m_explicit = Hipace::m_explicit; diff --git a/src/particles/deposition/PlasmaDepositCurrent.cpp b/src/particles/deposition/PlasmaDepositCurrent.cpp index 35384e4515..e37d94675f 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.cpp +++ b/src/particles/deposition/PlasmaDepositCurrent.cpp @@ -131,7 +131,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, [=] AMREX_GPU_DEVICE (int ip, auto ptd, auto /*depos_order*/, auto /*can_ionize*/, - auto /*use_laserr*/) + auto /*use_laser*/) { // only deposit plasma currents on or below their according MR level return ptd.id(ip).is_valid() && (lev == 0 || ptd.cpu(ip) >= lev); diff --git a/src/particles/deposition/TemperatureDeposition.cpp b/src/particles/deposition/TemperatureDeposition.cpp index 282b406563..5f7ecbfb9d 100644 --- a/src/particles/deposition/TemperatureDeposition.cpp +++ b/src/particles/deposition/TemperatureDeposition.cpp @@ -54,22 +54,60 @@ DepositTemperature (PlasmaParticleContainer& plasma, // extract laser properties and boolean for the presence of the laser and for ionization const PhysConst pc = get_phys_const(); const int aabs = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1; - const bool can_ionize = plasma.m_can_ionize; - const bool use_laser = Hipace::m_use_laser; const amrex::Real laser_norm = (plasma.m_charge/pc.q_e) * (pc.m_e/plasma.m_mass) * (plasma.m_charge/pc.q_e) * (pc.m_e/plasma.m_mass); // Loop over particles - SharedMemoryDeposition<3, 3, true>( - int(pti.numParticles()), + amrex::AnyCTO( + // use compile-time options + amrex::TypeList< + amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order + amrex::CompileTimeOptions, // can_ionize + amrex::CompileTimeOptions // use_laser + >{}, { + Hipace::m_temperature_depos_order, + plasma.m_can_ionize, + Hipace::m_use_laser + }, + // call deposition function + // The three functions passed as arguments to this lambda + // are defined below as the next arguments. + [&](auto is_valid, auto get_start_cell, auto deposit){ + constexpr auto ctos = deposit.GetOptions(); + constexpr int depos_order = ctos[0]; + constexpr int use_laser = ctos[2]; + constexpr int stencil_size = depos_order + 1; + + if constexpr (use_laser) { + SharedMemoryDeposition( + int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), + isl_fab.box(), pti.GetParticleTile().getParticleTileData(), + amrex::GpuArray{aabs}, + amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); + } else { + SharedMemoryDeposition( + int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), + isl_fab.box(), pti.GetParticleTile().getParticleTileData(), + amrex::GpuArray{}, + amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); + } + }, // is_valid // return whether the particle is valid and should deposit - [=] AMREX_GPU_DEVICE (int ip, auto ptd) + [=] AMREX_GPU_DEVICE (int ip, auto ptd, + auto /*depos_order*/, + auto /*can_ionize*/, + auto /*use_laser*/) { // only deposit on or below their according MR level return ptd.id(ip).is_valid() && (lev == 0 || ptd.cpu(ip) >= lev); }, - [=] AMREX_GPU_DEVICE (int ip, auto ptd) -> amrex::IntVectND<2> + // get_start_cell + // return the lowest cell index that the particle deposits into + [=] AMREX_GPU_DEVICE (int ip, auto ptd, + auto depos_order, + auto /*can_ionize*/, + auto /*use_laser*/) -> amrex::IntVectND<2> { const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); @@ -78,17 +116,21 @@ DepositTemperature (PlasmaParticleContainer& plasma, const amrex::Real ymid = (yp - y_pos_offset) * dy_inv; auto [shape_x, i] = - compute_single_shape_factor(xmid, 0); + compute_single_shape_factor(xmid, 0); auto [shape_y, j] = - compute_single_shape_factor(ymid, 0); + compute_single_shape_factor(ymid, 0); return {i-1, j-1}; }, + // do_deposit // deposit of weight, momentum (ux, uy, uz) and their squares (uxsq, uysq, uzsq) [=] AMREX_GPU_DEVICE (int ip, auto ptd, - Array3 arr, - auto cache_idx, auto depos_idx) noexcept + Array3 arr, + auto cache_idx, auto depos_idx, + auto depos_order, + auto can_ionize, + auto use_laser) noexcept { const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); @@ -101,7 +143,7 @@ DepositTemperature (PlasmaParticleContainer& plasma, ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; } doLaserGatherShapeN<2>(xp, yp, Aabssqp, arr, cache_idx[0], - dx_inv, dy_inv, x_pos_offset, y_pos_offset); + dx_inv, dy_inv, x_pos_offset, y_pos_offset); Aabssqp *= laser_norm_ion; } @@ -116,25 +158,24 @@ DepositTemperature (PlasmaParticleContainer& plasma, const amrex::Real xmid = (xp - x_pos_offset) * dx_inv; const amrex::Real ymid = (yp - y_pos_offset) * dy_inv; - // --- Compute shape factors - // x direction - auto [shape_x, i] = compute_single_shape_factor(xmid, 0); - // y direction - auto [shape_y, j] = compute_single_shape_factor(ymid, 0); - - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[0]), wp); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[1]), wp*uxp); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[2]), wp*uyp); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[3]), wp*uzp); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[4]), wp*uxp*uxp); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[5]), wp*uyp*uyp); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[6]), wp*uzp*uzp); - }, - isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{aabs}, - amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq} - ); + for (int iy=0; iy <= depos_order; ++iy) { + for (int ix=0; ix <= depos_order; ++ix) { + // --- Compute shape factors + // x direction + auto [shape_x, i] = compute_single_shape_factor(xmid, ix); + // y direction + auto [shape_y, j] = compute_single_shape_factor(ymid, iy); + + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[0]), shape_x*shape_y*wp); + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[1]), shape_x*shape_y*wp*uxp); + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[2]), shape_x*shape_y*wp*uyp); + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[3]), shape_x*shape_y*wp*uzp); + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[4]), shape_x*shape_y*wp*uxp*uxp); + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[5]), shape_x*shape_y*wp*uyp*uyp); + amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[6]), shape_x*shape_y*wp*uzp*uzp); + } + } + }); Array3 field_arr = isl_fab.array(); // Normalize the components of momentum (ux, uy, uz) and their squares (uxsq, uysq, uzsq) diff --git a/tests/laser_ionization.1Rank.sh b/tests/laser_ionization.1Rank.sh index eab7146909..83499c80ce 100755 --- a/tests/laser_ionization.1Rank.sh +++ b/tests/laser_ionization.1Rank.sh @@ -30,6 +30,7 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_laser_ionization \ my_constants.a0 = 0.00885126 \ plasmas.do_push = 0 \ hipace.file_prefix=$TEST_NAME/linear \ + hipace.temperature_depos_order = 0 \ plasmas.insitu_file_prefix = $TEST_NAME/insitu_linear mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_laser_ionization \ @@ -37,6 +38,7 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_laser_ionization \ lasers.polarization = circular \ plasmas.do_push = 0 \ hipace.file_prefix=$TEST_NAME/circular \ + hipace.temperature_depos_order = 0 \ plasmas.insitu_file_prefix = $TEST_NAME/insitu_circular # Compare the result with theory