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

Output field coarsening #580

Merged
merged 8 commits into from
Sep 7, 2021
Merged

Output field coarsening #580

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
4203ddd
remove FtoS
AlexanderSinn Jul 28, 2021
363025f
remove FtoS
AlexanderSinn Jul 28, 2021
f8bb55f
add coarsening
AlexanderSinn Jul 30, 2021
bc702a0
lf
AlexanderSinn Jul 30, 2021
c58def4
fix doc and max
AlexanderSinn Jul 30, 2021
4937912
fix doc
AlexanderSinn Jul 30, 2021
5f10ebd
do field subsampling
AlexanderSinn Aug 9, 2021
50cd89c
add test
Aug 10, 2021
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
4 changes: 4 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 @@ -319,6 +319,10 @@ Diagnostic parameters
x-axis, respectively. In case of an even number of grid points, the value will be averaged
between the two inner grid points.

* ``diagnostic.coarsening`` (3 `int`) optional (default `1 1 1`)
Coarsening ratio of field output in x, y and z direction respectively. For x and y, the average
value is taken. For z, the extra slices are skipped in the output.

* ``diagnostic.field_data`` (`string`) optional (default `all`)
Names of the fields written to file, separated by a space. The field names need to be `all`,
`none` or a subset of `ExmBy EypBx Ez Bx By Bz jx jy jz jx_beam jy_beam jz_beam rho Psi`.
Expand Down
4 changes: 2 additions & 2 deletions src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1286,8 +1286,8 @@ Hipace::ResizeFDiagFAB (const int it)
void
Hipace::FillDiagnostics (const int lev, int i_slice)
{
m_fields.Copy(lev, i_slice, FieldCopyType::StoF, 0, 0, Comps[WhichSlice::This]["N"],
m_diags.getF(lev), m_diags.sliceDir(), Geom(lev));
m_fields.Copy(lev, i_slice, 0, 0, Comps[WhichSlice::This]["N"],
m_diags.getF(lev), m_diags.sliceDir(), m_diags.getCoarsening(lev), Geom(lev));
}

void
Expand Down
4 changes: 4 additions & 0 deletions src/diagnostics/Diagnostic.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,9 @@ public:
/** return slice direction of the diagnostics */
int sliceDir () {return m_slice_dir;}

/** return coarsening ratio of diagnostics */
amrex::IntVect getCoarsening (int lev) { return m_diag_coarsen[lev]; }

/** \brief return box which possibly was trimmed in case of slice IO
*
* \param[in] box_3d box to be possibly trimmed to a slice box
Expand All @@ -65,6 +68,7 @@ private:
amrex::Vector<amrex::FArrayBox> m_F;
DiagType m_diag_type; /**< Type of diagnostics (xyz xz yz) */
int m_slice_dir; /**< Slicing direction */
amrex::Vector<amrex::IntVect> m_diag_coarsen; /**< xyz coarsening ratio (positive) */
amrex::Vector<std::string> m_comps_output; /**< Component names to Write to output file */
amrex::Vector<std::string> m_output_beam_names; /**< Component names to Write to output file */
int m_nfields; /**< Number of physical fields to write */
Expand Down
19 changes: 18 additions & 1 deletion src/diagnostics/Diagnostic.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

Diagnostic::Diagnostic (int nlev)
: m_F(nlev),
m_diag_coarsen(nlev),
m_geom_io(nlev)
{
amrex::ParmParse ppd("diagnostic");
Expand All @@ -22,6 +23,18 @@ Diagnostic::Diagnostic (int nlev)
amrex::Abort("Unknown diagnostics type: must be xyz, xz or yz.");
}

for(int ilev = 0; ilev<nlev; ++ilev) {
amrex::Array<int,3> diag_coarsen_arr{1,1,1};
// set all levels the same for now
ppd.query("coarsening", diag_coarsen_arr);
if(m_slice_dir == 0 || m_slice_dir == 1) {
diag_coarsen_arr[m_slice_dir] = 1;
}
m_diag_coarsen[ilev] = amrex::IntVect(diag_coarsen_arr);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_diag_coarsen[ilev].min() >= 1,
"Coarsening ratio must be >= 1");
}

ppd.queryarr("field_data", m_comps_output);
const amrex::Vector<std::string> all_field_comps
{"ExmBy", "EypBx", "Ez", "Bx", "By", "Bz", "jx", "jx_beam", "jy", "jy_beam", "jz",
Expand Down Expand Up @@ -82,6 +95,8 @@ Diagnostic::AllocData (int lev, const amrex::Box& bx, int nfields, amrex::Geomet
// trim the 3D box to slice box for slice IO
amrex::Box F_bx = TrimIOBox(bx);

F_bx.coarsen(m_diag_coarsen[lev]);

m_F.push_back(amrex::FArrayBox(F_bx, m_nfields, amrex::The_Pinned_Arena()));

m_geom_io[lev] = geom;
Expand All @@ -94,14 +109,16 @@ Diagnostic::AllocData (int lev, const amrex::Box& bx, int nfields, amrex::Geomet
domain.setBig(m_slice_dir, icenter);
m_geom_io[lev] = amrex::Geometry(domain, &prob_domain, geom.Coord());
}
m_geom_io[lev].coarsen(m_diag_coarsen[lev]);
}

void
Diagnostic::ResizeFDiagFAB (const amrex::Box box, const int lev)
{
amrex::Box io_box = TrimIOBox(box);
io_box.coarsen(m_diag_coarsen[lev]);
m_F[lev].resize(io_box, m_nfields);
}
}

amrex::Box
Diagnostic::TrimIOBox (const amrex::Box box_3d)
Expand Down
10 changes: 4 additions & 6 deletions src/fields/Fields.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,6 @@

class Hipace;

/** \brief Direction of field copies: from 3D F to 2D slice S or the other way round */
enum struct FieldCopyType { FtoS, StoF };

/** \brief describes which slice with respect to the currently calculated is used */
struct WhichSlice {
enum slice { Next=0, This, Previous1, Previous2, RhoIons, N };
Expand Down Expand Up @@ -104,16 +101,17 @@ public:
*
* \param[in] lev MR level
* \param[in] i_slice z slice in which to write the data
* \param[in] copy_type whether to copy from the xy slice to the main array or opposite
* \param[in] slice_comp first component of the xy slice to copy from/to
* \param[in] full_comp first component of the full array to copy from/to
* \param[in] ncomp number of components to copy
* \param[in,out] fab full FArrayBox
* \param[in] slice_dir slicing direction. 0=x, 1=y, -1=no slicing (full 3D)
* \param[in] diag_coarsen coarsening ratio of diagnostics
* \param[in] geom main geometry
*/
void Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int full_comp,
int ncomp, amrex::FArrayBox& fab, int slice_dir, amrex::Geometry geom);
void Copy (int lev, int i_slice, int slice_comp, int full_comp,
int ncomp, amrex::FArrayBox& fab, int slice_dir,
amrex::IntVect diag_coarsen, amrex::Geometry geom);

/** \brief Shift slices by 1 element: slices (1,2) are then stored in (2,3).
*
Expand Down
81 changes: 59 additions & 22 deletions src/fields/Fields.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -233,12 +233,13 @@ Fields::LongitudinalDerivative (const amrex::MultiFab& src1, const amrex::MultiF


void
Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int full_comp,
int ncomp, amrex::FArrayBox& fab, int slice_dir, amrex::Geometry geom)
Fields::Copy (int lev, int i_slice, int slice_comp, int full_comp,
int ncomp, amrex::FArrayBox& fab, int slice_dir,
amrex::IntVect diag_coarsen, amrex::Geometry geom)
{
using namespace amrex::literals;
HIPACE_PROFILE("Fields::Copy()");
auto& slice_mf = m_slices[lev][WhichSlice::This]; // copy from/to the current slice
auto& slice_mf = m_slices[lev][WhichSlice::This]; // copy from the current slice
amrex::Array4<amrex::Real> slice_array; // There is only one Box.
for (amrex::MFIter mfi(slice_mf); mfi.isValid(); ++mfi) {
auto& slice_fab = slice_mf[mfi];
Expand All @@ -250,42 +251,78 @@ Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int
}

amrex::Box const& vbx = fab.box();
if (vbx.smallEnd(Direction::z) <= i_slice and
vbx.bigEnd (Direction::z) >= i_slice)
if (vbx.smallEnd(Direction::z) <= i_slice / diag_coarsen[2] and
vbx.bigEnd (Direction::z) >= i_slice / diag_coarsen[2] and
i_slice % diag_coarsen[2] == 0)
{
amrex::Box copy_box = vbx;
copy_box.setSmall(Direction::z, i_slice);
copy_box.setBig (Direction::z, i_slice);
copy_box.setSmall(Direction::z, i_slice / diag_coarsen[2]);
copy_box.setBig (Direction::z, i_slice / diag_coarsen[2]);

amrex::Array4<amrex::Real> const& full_array = fab.array();

const amrex::IntVect ncells_global = geom.Domain().length();
const bool nx_even = ncells_global[0] % 2 == 0;
const bool ny_even = ncells_global[1] % 2 == 0;

if (copy_type == FieldCopyType::FtoS) {
const int coarse_x = diag_coarsen[0];
const int coarse_y = diag_coarsen[1];
const int coarse_z = diag_coarsen[2];

const int ncells_x = ncells_global[0];
const int ncells_y = ncells_global[1];

if (slice_dir ==-1 /* 3D data */){

amrex::ParallelFor(copy_box, ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
slice_array(i,j,k,n+slice_comp) = full_array(i,j,k,n+full_comp);
amrex::Real field_value = 0;
int n_values = 0;
for (int j_c = 0; j_c < coarse_y && (j*coarse_y+j_c) < ncells_y; ++j_c) {
for (int i_c = 0; i_c < coarse_x && (i*coarse_x+i_c) < ncells_x; ++i_c) {
field_value += slice_array( i*coarse_x+i_c,
j*coarse_y+j_c,
k*coarse_z,
n+slice_comp);
++n_values;
}
}
full_array(i,j,k,n+full_comp) = field_value / std::max(n_values,1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand correctly, this is large-order smoothing interpolation (where the stencil length is coarse_* for direction *), where the same coefficient is applies to all point of the stencil. This causes significant smoothing of the data, which could be unwanted. Besides, we probably don't have to re-write interpolation functions, we already have some (originally written for particles, but also used for field interpolation in #572).

In this case, I would recommend to do sub-sampling rather than large-stencil interpolation, where we would ideally just take 1 point in each direction. Of course, this is impractical because we want to conserve the centering of data. For this I would recommend to do as shown on the image below: interpolate between the two nearest grid points (in each direction), which is the minimum to fix the centering (this is what we do for slice I/O). Do you think you could implement this instead? Of course this is valid for x and y, and for 3D slice I/O.

IMG_20210809_104848

});
} else {

} else if (slice_dir == 0 /* yz slice */){

amrex::ParallelFor(copy_box, ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
amrex::Real field_value = 0;
int n_values = 0;
for (int j_c = 0; j_c < coarse_y && (j*coarse_y+j_c) < ncells_y; ++j_c) {
field_value += nx_even ?
0.5_rt * (slice_array(i-1,j*coarse_y+j_c,k*coarse_z,n+slice_comp) +
slice_array(i ,j*coarse_y+j_c,k*coarse_z,n+slice_comp))
: slice_array(i ,j*coarse_y+j_c,k*coarse_z,n+slice_comp);
++n_values;
}
full_array(i,j,k,n+full_comp) = field_value / std::max(n_values,1);
});

} else /* slice_dir == 1, xz slice */{

amrex::ParallelFor(copy_box, ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
if (slice_dir ==-1 /* 3D data */){
full_array(i,j,k,n+full_comp) = slice_array(i,j,k,n+slice_comp);
} else if (slice_dir == 0 /* yz slice */){
full_array(i,j,k,n+full_comp) =
nx_even ? 0.5_rt * (slice_array(i-1,j,k,n+slice_comp) +
slice_array(i,j,k,n+slice_comp))
: slice_array(i,j,k,n+slice_comp);
} else /* slice_dir == 1, xz slice */{
full_array(i,j,k,n+full_comp) =
ny_even ? 0.5_rt * ( slice_array(i,j-1,k,n+slice_comp) +
slice_array(i,j,k,n+slice_comp))
: slice_array(i,j,k,n+slice_comp);
amrex::Real field_value = 0;
int n_values = 0;
for (int i_c = 0; i_c < coarse_x && (i*coarse_x+i_c) < ncells_x; ++i_c) {
field_value += ny_even ?
0.5_rt * (slice_array(i*coarse_x+i_c,j-1,k*coarse_z,n+slice_comp) +
slice_array(i*coarse_x+i_c,j ,k*coarse_z,n+slice_comp))
: slice_array(i*coarse_x+i_c,j ,k*coarse_z,n+slice_comp);
++n_values;
}
full_array(i,j,k,n+full_comp) = field_value / std::max(n_values,1);
});
}
}
Expand Down