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

Interpolate jx_beam and jy_beam to level 1 #936

Merged
merged 1 commit into from
May 7, 2023
Merged

Interpolate jx_beam and jy_beam to level 1 #936

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
9 changes: 9 additions & 0 deletions src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -524,6 +524,15 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step)
if (islice < m_3D_geom[lev].Domain().smallEnd(Direction::z) ||
islice > m_3D_geom[lev].Domain().bigEnd(Direction::z)) {
continue;
} else if (islice == m_3D_geom[lev].Domain().bigEnd(Direction::z)) {
// first slice of level 1 (islice goes backwards)
// iterpolate jx_beam and jy_beam from level 0 to level 1
m_fields.LevelUp(m_3D_geom, lev, WhichSlice::Previous1, "jx_beam");
m_fields.LevelUp(m_3D_geom, lev, WhichSlice::Previous1, "jy_beam");
m_fields.LevelUp(m_3D_geom, lev, WhichSlice::This, "jx_beam");
m_fields.LevelUp(m_3D_geom, lev, WhichSlice::This, "jy_beam");
m_fields.duplicate(lev, WhichSlice::This, {"jx" , "jy" },
WhichSlice::This, {"jx_beam", "jy_beam"});
}
}

Expand Down
12 changes: 11 additions & 1 deletion src/fields/Fields.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,14 +183,24 @@ public:
*
* \param[in] geom Geometry
* \param[in] lev current level
* \param[in] component which can be Psi or rho
* \param[in] component which can be Psi or rho etc.
* \param[in] outer_edge start writing interpolated values at domain + outer_edge
* \param[in] inner_edge stop writing interpolated values at domain + inner_edge
*/
void InterpolateFromLev0toLev1 (amrex::Vector<amrex::Geometry> const& geom, const int lev,
std::string component, const amrex::IntVect outer_edge,
const amrex::IntVect inner_edge);

/** \brief Interpolate the full field from the coarse grid (lev-1) to the fine grid (lev).
*
* \param[in] geom Geometry
* \param[in] lev current level
* \param[in] component which can be jx_beam or jy_beam etc.
* \param[in] which_slice slice of the field to interpolate
*/
void LevelUp (amrex::Vector<amrex::Geometry> const& geom, const int lev,
const int which_slice, const std::string component);

/** \brief Compute ExmBy and EypBx on the slice container from J by solving a Poisson equation
* ExmBy and EypBx are solved in the same function because both rely on Psi.
*
Expand Down
34 changes: 34 additions & 0 deletions src/fields/Fields.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -695,6 +695,40 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector<amrex::Geometry> const& geom, c
}


void
Fields::LevelUp (amrex::Vector<amrex::Geometry> const& geom, const int lev,
const int which_slice, const std::string component)
{
if (lev == 0) return; // only interpolate field to lev 1
HIPACE_PROFILE("Fields::LevelUp()");
constexpr int interp_order = 2;
Copy link
Member

Choose a reason for hiding this comment

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

A interp_order stencil was already applied from beam particles to the grid. this is an additional smoothing step that is unnecessary, in particular because we are interpolating to a fine grid to resolve finer structures. A lower interpolation order here should reduce the difference with jx_beam obtained by depositing directly to the fine grid.

Suggested change
constexpr int interp_order = 2;
constexpr int interp_order = 1;

Copy link
Member Author

Choose a reason for hiding this comment

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

Currently all interpolation to level one is done with order two. Should the other functions be changed as well or just this one.

Copy link
Member

Choose a reason for hiding this comment

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

Probably all of them, but we should review this case by case (might sometimes depend on whether it's a source or a field). But you're right, it's not only about this one. Could you open a separate PR setting all others to e.g. 1 (or have a variable controlling them, if that's quick, but we'll probably scan it just in the next few days/weeks before merging the PR and set it to the best value)? Let me know if that works for you, and if so we can merge this PR as is for consistency.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, it is easier to change and test everything in one go after merging this PR with order two.

Copy link
Member

Choose a reason for hiding this comment

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

Let me take it back: we don't need to scan it. I believe that hard-coding it to 1 everywhere is what makes most sense. Then we could just compare old (2) vs. new (1) and see if it's a significant change. But I see no reason how 0 (interpolating from wrong position) or 3 (huge smoothing, as seen from level 1) could ever make sense. So no need to expose it as a variable (unless it's easier in the code).

Copy link
Member

Choose a reason for hiding this comment

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

ok

Copy link
Member Author

Choose a reason for hiding this comment

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

Let me take it back: we don't need to scan it. I believe that hard-coding it to 1 everywhere is what makes most sense. Then we could just compare old (2) vs. new (1) and see if it's a significant change. But I see no reason how 0 (interpolating from wrong position) or 3 (huge smoothing, as seen from level 1) could ever make sense. So no need to expose it as a variable (unless it's easier in the code).

I agree


auto field_coarse_interp = interpolated_field_xy<interp_order, amrex::MultiFab>{
getField(lev-1, which_slice, component), geom[lev-1]};
amrex::MultiFab field_fine = getField(lev, which_slice, component);

for (amrex::MFIter mfi( field_fine, DfltMfi); mfi.isValid(); ++mfi)
{
auto arr_field_coarse_interp = field_coarse_interp.array(mfi);
const Array2<amrex::Real> arr_field_fine = field_fine.array(mfi);

const amrex::Real dx = geom[lev].CellSize(0);
const amrex::Real dy = geom[lev].CellSize(1);
const amrex::Real offset0 = GetPosOffset(0, geom[lev], geom[lev].Domain());
const amrex::Real offset1 = GetPosOffset(1, geom[lev], geom[lev].Domain());

amrex::ParallelFor(field_fine[mfi].box(),
[=] AMREX_GPU_DEVICE (int i, int j , int) noexcept
{
// interpolate the full field
const amrex::Real x = i * dx + offset0;
const amrex::Real y = j * dy + offset1;
arr_field_fine(i,j) = arr_field_coarse_interp(x,y);
});
}
}


void
Fields::SolvePoissonExmByAndEypBx (amrex::Vector<amrex::Geometry> const& geom, const int lev)
{
Expand Down