Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean-up in-situ diagnostic #1052

Merged
merged 3 commits into from
Jan 10, 2024
Merged

Clean-up in-situ diagnostic #1052

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d21b6bf
clean up insitu diagnostic
AlexanderSinn Dec 18, 2023
a3ad446
use amrex::constexpr_for
AlexanderSinn Dec 19, 2023
dc55b9a
add comments
AlexanderSinn Jan 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions docs/source/run/parameters.rst
View file
Open in desktop
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -529,6 +529,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
Open in desktop
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 @@ -723,12 +724,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 @@ -1286,14 +1282,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 @@ -1302,40 +1292,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
Open in desktop
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
Open in desktop
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
Loading