Skip to content

Commit

Permalink
Clean-up in-situ diagnostic (#1052)
Browse files Browse the repository at this point in the history 
* clean up insitu diagnostic

* use amrex::constexpr_for

* add comments
  • Loading branch information
@AlexanderSinn
AlexanderSinn authored Jan 10, 2024
1 parent c40f0fb commit c41482d
Show file tree
Hide file tree
Showing 9 changed files with 214 additions and 299 deletions.
5 changes: 5 additions & 0 deletions docs/source/run/parameters.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -932,6 +932,11 @@ Thereby, "[]" stands for averaging over all particles in the current slice,
Averages and totals over all slices are also provided for convenience under the
respective ``average`` and ``total`` subcategories.

For the field in-situ diagnostics, the following quantities are calculated per slice and stored:
``[Ex^2], [Ey^2], [Ez^2], [Bx^2], [By^2], [Bz^2], [ExmBy^2], [EypBx^2], [Ez*jz_beam]``.
Thereby, "[]" stands for averaging over all cells in the current slice.
These quantities can be used to calculate the energy stored in the fields.

Additionally, some metadata is also available:
``time, step, n_slices, charge, mass, z_lo, z_hi, normalized_density_factor``.
``time`` and ``step`` refers to the physical time of the simulation and step number of the
Expand Down
2 changes: 2 additions & 0 deletions src/fields/Fields.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -521,6 +521,8 @@ private:
bool m_explicit = false;
/** If any plasma species has a neutralizing background */
bool m_any_neutral_background = false;
/** Number of real field properties for in-situ per-slice reduced diagnostics. */
static constexpr int m_insitu_nrp = 9;
/** How often the insitu field diagnostics should be computed and written
* Default is 0, meaning no output */
int m_insitu_period {0};
Expand Down
66 changes: 22 additions & 44 deletions src/fields/Fields.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include "utils/Constants.H"
#include "utils/GPUUtil.H"
#include "utils/InsituUtil.H"
#include "utils/TemplateUtil.H"
#include "particles/particles_utils/ShapeFactors.H"
#ifdef HIPACE_USE_OPENPMD
# include <openPMD/auxiliary/Filesystem.hpp>
Expand Down Expand Up @@ -205,8 +206,8 @@ Fields::AllocData (
"Must choose a different field insitu file prefix compared to the full diagnostics");
#endif
// Allocate memory for in-situ diagnostics
m_insitu_rdata.resize(geom.Domain().length(2)*9, 0.);
m_insitu_sum_rdata.resize(9, 0.);
m_insitu_rdata.resize(geom.Domain().length(2)*m_insitu_nrp, 0.);
m_insitu_sum_rdata.resize(m_insitu_nrp, 0.);
}
}

Expand Down Expand Up @@ -725,12 +726,7 @@ Fields::SetBoundaryCondition (amrex::Vector<amrex::Geometry> const& geom, const
const amrex::Real x = (i * dx + poff_x) * scale;
const amrex::Real y = (j * dy + poff_y) * scale;
if (x*x + y*y > cutoff_sq) {
return MultipoleTuple{0._rt,
0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt,
0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt,
0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt,
0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt
};
return MakeZeroTuple(MultipoleTuple{});
}
amrex::Real s_v = arr_staging_area(i, j);
return GetMultipoleCoeffs(s_v, x, y);
Expand Down Expand Up @@ -1288,14 +1284,8 @@ Fields::InSituComputeDiags (int step, amrex::Real time, int islice, const amrex:

amrex::MultiFab& slicemf = getSlices(lev);

// Tuple contains:
// 0, 1, 2, 3, 4, 5, 6, 7 8
// <Ex^2>, <Ey^2>, <Ez^2>, <Bx^2>, <By^2>, <Bz^2>, <ExmBy^2>, <EypBx^2>, <Ez*jz_beam>
amrex::ReduceOps<amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum> reduce_op;
amrex::ReduceData<amrex::Real, amrex::Real, amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real, amrex::Real> reduce_data(reduce_op);
TypeMultiplier<amrex::ReduceOps, amrex::ReduceOpSum[m_insitu_nrp]> reduce_op;
TypeMultiplier<amrex::ReduceData, amrex::Real[m_insitu_nrp]> reduce_data(reduce_op);
using ReduceTuple = typename decltype(reduce_data)::Type;

for ( amrex::MFIter mfi(slicemf, DfltMfi); mfi.isValid(); ++mfi ) {
Expand All @@ -1304,40 +1294,28 @@ Fields::InSituComputeDiags (int step, amrex::Real time, int islice, const amrex:
mfi.tilebox(), reduce_data,
[=] AMREX_GPU_DEVICE (int i, int j, int) -> ReduceTuple
{
return {
pow<2>(arr(i,j,ExmBy) + arr(i,j,By) * clight),
pow<2>(arr(i,j,EypBx) - arr(i,j,Bx) * clight),
pow<2>(arr(i,j,Ez)),
pow<2>(arr(i,j,Bx)),
pow<2>(arr(i,j,By)),
pow<2>(arr(i,j,Bz)),
pow<2>(arr(i,j,ExmBy)),
pow<2>(arr(i,j,EypBx)),
arr(i,j,Ez)*arr(i,j,jz_beam)
return { // Tuple contains:
pow<2>(arr(i,j,ExmBy) + arr(i,j,By) * clight), // 0 [Ex^2]
pow<2>(arr(i,j,EypBx) - arr(i,j,Bx) * clight), // 1 [Ey^2]
pow<2>(arr(i,j,Ez)), // 2 [Ez^2]
pow<2>(arr(i,j,Bx)), // 3 [Bx^2]
pow<2>(arr(i,j,By)), // 4 [By^2]
pow<2>(arr(i,j,Bz)), // 5 [Bz^2]
pow<2>(arr(i,j,ExmBy)), // 6 [ExmBy^2]
pow<2>(arr(i,j,EypBx)), // 7 [EypBx^2]
arr(i,j,Ez)*arr(i,j,jz_beam) // 8 [Ez*jz_beam]
};
});
}

ReduceTuple a = reduce_data.value();

m_insitu_rdata[islice ] = amrex::get<0>(a)*dxdydz;
m_insitu_rdata[islice+1*nslices] = amrex::get<1>(a)*dxdydz;
m_insitu_rdata[islice+2*nslices] = amrex::get<2>(a)*dxdydz;
m_insitu_rdata[islice+3*nslices] = amrex::get<3>(a)*dxdydz;
m_insitu_rdata[islice+4*nslices] = amrex::get<4>(a)*dxdydz;
m_insitu_rdata[islice+5*nslices] = amrex::get<5>(a)*dxdydz;
m_insitu_rdata[islice+6*nslices] = amrex::get<6>(a)*dxdydz;
m_insitu_rdata[islice+7*nslices] = amrex::get<7>(a)*dxdydz;
m_insitu_rdata[islice+8*nslices] = amrex::get<8>(a)*dxdydz;
m_insitu_sum_rdata[0] += amrex::get<0>(a)*dxdydz;
m_insitu_sum_rdata[1] += amrex::get<1>(a)*dxdydz;
m_insitu_sum_rdata[2] += amrex::get<2>(a)*dxdydz;
m_insitu_sum_rdata[3] += amrex::get<3>(a)*dxdydz;
m_insitu_sum_rdata[4] += amrex::get<4>(a)*dxdydz;
m_insitu_sum_rdata[5] += amrex::get<5>(a)*dxdydz;
m_insitu_sum_rdata[6] += amrex::get<6>(a)*dxdydz;
m_insitu_sum_rdata[7] += amrex::get<7>(a)*dxdydz;
m_insitu_sum_rdata[8] += amrex::get<8>(a)*dxdydz;
amrex::constexpr_for<0, m_insitu_nrp>(
[&] (auto idx) {
m_insitu_rdata[islice + idx.value * nslices] = amrex::get<idx.value>(a)*dxdydz;
m_insitu_sum_rdata[idx.value] += amrex::get<idx.value>(a)*dxdydz;
}
);
}

void
Expand Down
49 changes: 5 additions & 44 deletions src/fields/OpenBoundary.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#ifndef OPEN_BOUNDARY_H_
#define OPEN_BOUNDARY_H_

#include "utils/TemplateUtil.H"

#include <AMReX_AmrCore.H>
#include <cmath>

Expand All @@ -27,50 +29,9 @@ amrex::Real pow (amrex::Real base) {
return 0._rt; //shut up compiler
}

using MultipoleTuple = amrex::GpuTuple<
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real>;

using MultipoleReduceOpList = amrex::TypeList<
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
amrex::ReduceOpSum>;

using MultipoleReduceTypeList = amrex::TypeList<
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real, amrex::Real, amrex::Real,
amrex::Real>;
using MultipoleTuple = TypeMultiplier<amrex::GpuTuple, amrex::Real[37]>;
using MultipoleReduceOpList = TypeMultiplier<amrex::TypeList, amrex::ReduceOpSum[37]>;
using MultipoleReduceTypeList = TypeMultiplier<amrex::TypeList, amrex::Real[37]>;

// To solve a poisson equation (d^2/dx^2 + d^2/dy^2)phi = source with open boundary conditions for
// phi(x,y), the source field at (x',y') is integrated together with the Green's function
Expand Down
23 changes: 4 additions & 19 deletions src/particles/beam/BeamParticleContainer.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -301,27 +301,12 @@ private:

int m_nslices; /**< number of z slices of the domain */
/** Number of real beam properties for in-situ per-slice reduced diagnostics. */
int m_insitu_nrp {18};
static constexpr int m_insitu_nrp = 18;
/** Number of int beam properties for in-situ per-slice reduced diagnostics. */
int m_insitu_nip {1};
/** Per-slice real beam properties:
* 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
* sum(w), [x], [x^2], [y], [y^2], [z], [z^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2],
*
* 13, 14, 15, 16, 17
* [x*ux], [y*uy], [z*uz], [ga], [ga^2]
* where [] means average over all particles within slice.
* Per-slice emittance: sqrt( abs( ([x^2]-[x]^2) * ([ux^2]-[ux]^2) - ([x*ux]-[x][ux])^2 ) ).
* Projected emittance: Same as above AFTER averaging all these quantities over slices.
* Energy spread: sqrt([ga^2]-[ga]^2), and same as above.
* Momentum spread: sqrt([uz^2]-[uz]^2), and same as above.
*/
static constexpr int m_insitu_nip = 1;
/** Per-slice real beam properties */
amrex::Vector<amrex::Real> m_insitu_rdata;
/** Per-slice int beam properties:
* 0
* Np
* Np: number of particles in this slice
*/
/** Per-slice int beam properties */
amrex::Vector<int> m_insitu_idata;
/** Sum of all per-slice real beam properties */
amrex::Vector<amrex::Real> m_insitu_sum_rdata;
Expand Down
Loading

0 comments on commit c41482d

Please sign in to comment.