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.

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

Implement linear interpolation for both cell-centered and nodal data types #3638

Merged
merged 11 commits into from
Jan 31, 2023

Conversation

prkkumar
Copy link
Member

@prkkumar prkkumar commented Jan 25, 2023

This PR implements linear interpolation for both cell-centered and nodal arrays. Before this, we were doing nearest grid point interpolation for cell-centered directions, which was resulting in "stair-casing" in the aux field. This effect was particularly intensified with larger mesh refinement ratio (e.g. 4).

@prkkumar prkkumar marked this pull request as ready for review January 25, 2023 04:24
@prkkumar prkkumar requested review from EZoni and RemiLehe January 25, 2023 04:24
@prkkumar prkkumar changed the title use two points and correct weights for higher order interpolation Implement higher order interpolation for both cell-centered and nodal data types Jan 25, 2023
@prkkumar prkkumar changed the title Implement higher order interpolation for both cell-centered and nodal data types Implement linear interpolation for both cell-centered and nodal data types Jan 25, 2023
@EZoni
Copy link
Member

EZoni commented Jan 25, 2023

Closing and re-opening to trigger CI.

@EZoni EZoni closed this Jan 25, 2023
@EZoni EZoni reopened this Jan 25, 2023
@EZoni EZoni self-assigned this Jan 25, 2023
Comment on lines 60 to 70
wj = (sj == 0) ? 1.0_rt : (rj - amrex::Math::abs(j - (jc + jj) * rj))
wj = (sj == 0) ? (rj - amrex::Math::abs(j + 0.5 - (jc + jj + 0.5) * rj))
/ static_cast<amrex::Real>(rj)
: (rj - amrex::Math::abs(j - (jc + jj) * rj))
/ static_cast<amrex::Real>(rj);
wk = (sk == 0) ? 1.0_rt : (rk - amrex::Math::abs(k - (kc + kk) * rk))
wk = (sk == 0) ? (rk - amrex::Math::abs(k + 0.5 - (kc + kk + 0.5) * rk))
/ static_cast<amrex::Real>(rk)
: (rk - amrex::Math::abs(k - (kc + kk) * rk))
/ static_cast<amrex::Real>(rk);
wl = (sl == 0) ? 1.0_rt : (rl - amrex::Math::abs(l - (lc + ll) * rl))
wl = (sl == 0) ? (rl - amrex::Math::abs(l + 0.5 - (lc + ll + 0.5) * rl))
/ static_cast<amrex::Real>(rl)
: (rl - amrex::Math::abs(l - (lc + ll) * rl))
Copy link
Member

@EZoni EZoni Jan 25, 2023

Choose a reason for hiding this comment

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

What do you think about aligning the code in the following way? I feel like it'd make it quite more readable (it doesn't matter too much if we go over 100 characters in these few lines):

wj = (sj == 0) ? (rj - amrex::Math::abs(j + 0.5_rt - (jc + jj + 0.5_rt) * rj)) / static_cast<amrex::Real>(rj)
               : (rj - amrex::Math::abs(j          - (jc + jj         ) * rj)) / static_cast<amrex::Real>(rj);
wk = (sk == 0) ? (rk - amrex::Math::abs(k + 0.5_rt - (kc + kk + 0.5_rt) * rk)) / static_cast<amrex::Real>(rk)
               : (rk - amrex::Math::abs(k          - (kc + kk         ) * rk)) / static_cast<amrex::Real>(rk);
wl = (sl == 0) ? (rl - amrex::Math::abs(l + 0.5_rt - (lc + ll + 0.5_rt) * rl)) / static_cast<amrex::Real>(rl)
               : (rl - amrex::Math::abs(l          - (lc + ll         ) * rl)) / static_cast<amrex::Real>(rl);

Please remember to add the _rt suffix to the real constants, too.

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, that looks better. Thanks, @EZoni !

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
@EZoni EZoni self-requested a review January 25, 2023 21:49
Copy link
Member

@EZoni EZoni left a comment

Choose a reason for hiding this comment

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

Thank you for this PR, @prkkumar! I reset the CI benchmarks, obviously almost all of those using MR needed to be reset. I also commented one last suggestion on the interpolation function, let me know what you think.

Comment on lines 60 to 65
wj = (sj == 0) ? (rj - amrex::Math::abs(j + 0.5_rt - (jc + jj + 0.5_rt) * rj)) / static_cast<amrex::Real>(rj)
: (rj - amrex::Math::abs(j - (jc + jj ) * rj)) / static_cast<amrex::Real>(rj);
wk = (sk == 0) ? (rk - amrex::Math::abs(k + 0.5_rt - (kc + kk + 0.5_rt) * rk)) / static_cast<amrex::Real>(rk)
: (rk - amrex::Math::abs(k - (kc + kk ) * rk)) / static_cast<amrex::Real>(rk);
wl = (sl == 0) ? (rl - amrex::Math::abs(l + 0.5_rt - (lc + ll + 0.5_rt) * rl)) / static_cast<amrex::Real>(rl)
: (rl - amrex::Math::abs(l - (lc + ll ) * rl)) / static_cast<amrex::Real>(rl);
Copy link
Member

Choose a reason for hiding this comment

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

At this point, I wonder if it would be slightly better to remove the if conditions on sj, sk and sl from these nested loops, and have them only once instead, before the loops. sj, sk and sl do not vary depending on jj, kk and ll, so in principle there would be no need to check them within these nested loops.

I think we could do something like the following:

wj = (rj - amrex::Math::abs(j + hj - (jc + jj + hj) * rj)) / static_cast<amrex::Real>(rj);
wk = (rk - amrex::Math::abs(k + hk - (kc + kk + hk) * rk)) / static_cast<amrex::Real>(rk);
wl = (rl - amrex::Math::abs(l + hl - (lc + ll + hl) * rl)) / static_cast<amrex::Real>(rl);

where hj, hk and hl are defined, before the nested loops, as

amrex::Real hj = (sj == 0) ? 0.5_rt : 0._rt;
amrex::Real hk = (sk == 0) ? 0.5_rt : 0._rt;
amrex::Real hl = (sl == 0) ? 0.5_rt : 0._rt;

What do you think?

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, I like it this way. Thank you! I will make the change and push it. Also, thank you for resetting the CI benchmarks!

@EZoni EZoni self-requested a review January 26, 2023 18:39
@EZoni EZoni merged commit 3b46d91 into ECP-WarpX:development Jan 31, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants