Interpolation routines for tracers with mapped_z. (AMReX-Codes#3714)

This code is from @asalmgren and @RevathiJambunathan . The proposed changes: - [ ] fix a bug or incorrect behavior in AMReX - [x] add new capabilities to AMReX - [ ] changes answers in the test suite to more than roundoff level - [ ] are likely to significantly affect the results of downstream AMReX users - [ ] include documentation in the code and/or rst files, if appropriate --------- Co-authored-by: Revathi Jambunathan <[email protected]> Co-authored-by: Ann Almgren <[email protected]>
RevathiJambunathan · Jan 23, 2024 · 3e2a3c2 · 3e2a3c2
1 parent 73b2155
commit 3e2a3c2
Showing 1 changed file with 262 additions and 0 deletions.
diff --git a/Src/Particle/AMReX_TracerParticle_mod_K.H b/Src/Particle/AMReX_TracerParticle_mod_K.H
@@ -13,6 +13,12 @@
 
 namespace amrex {
 
+//
+// **********************************************************************
+// Regular coordinates
+// **********************************************************************
+//
+
 /**
  \brief Linearly interpolates the mesh data to the particle position from cell-centered data.
 */
@@ -135,5 +141,261 @@ void linear_interpolate_to_particle (const P& p,
  } // d
 }
 
+//
+// **********************************************************************
+// Terrain-fitted coordinates
+// **********************************************************************
+//
+
+/**
+ \brief Linearly interpolates the mesh data to the particle position from cell-centered data
+ on a terrain-fitted grid.
+*/
+template <typename P>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void cic_interpolate_mapped_z (const P& p,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+ const amrex::Array4<amrex::Real const>& data_arr,
+ const amrex::Array4<amrex::Real const>& height_arr,
+ amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
+{
+ cic_interpolate_cc_mapped_z(p, plo, dxi, data_arr, height_arr, val, M);
+}
+
+/**
+ \brief Linearly interpolates the mesh data to the particle position from cell-centered data
+ on a terrain-fitted grid.
+*/
+template <typename P>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void cic_interpolate_cc_mapped_z (const P& p,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+ const amrex::Array4<amrex::Real const>& data_arr,
+ const amrex::Array4<amrex::Real const>& height_arr,
+ amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
+{
+ int icomp = 0;
+ int ncomp_per_array = M;
+ int num_arrays = 1;
+ IntVect is_nodal = amrex::IntVect::TheZeroVector();
+ linear_interpolate_to_particle_z(p, plo, dxi, &data_arr, height_arr,
+ val, &is_nodal, icomp, ncomp_per_array, num_arrays);
+}
+
+/**
+ \brief Linearly interpolates the mesh data to the particle position from node-centered data.
+ on a terrain-fitted grid.
+*/
+template <typename P>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void cic_interpolate_nd_mapped_z (const P& p,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+ const amrex::Array4<amrex::Real const>& data_arr,
+ const amrex::Array4<amrex::Real const>& height_arr,
+ amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
+{
+ int icomp = 0;
+ int ncomp_per_array = M;
+ int num_arrays = 1;
+ IntVect is_nodal = amrex::IntVect::TheUnitVector();
+ linear_interpolate_to_particle_z(p, plo, dxi, &data_arr, height_arr,
+ val, &is_nodal, icomp, ncomp_per_array, num_arrays);
+}
+
+/**
+ \brief Linearly interpolates the mesh data to the particle position from face-centered data
+ on a terrain-fitted grid.
+ The nth component of the data_arr array is nodal in the nth direction, and cell-centered in the others.
+*/
+template <typename P>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void mac_interpolate_mapped_z (const P& p,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+ amrex::GpuArray<amrex::Array4<amrex::Real const>,AMREX_SPACEDIM> const& data_arr,
+ const amrex::Array4<amrex::Real const>& height_arr,
+ amrex::ParticleReal * val)
+{
+ int icomp = 0;
+ int ncomp_per_array = 1;
+ int num_arrays = AMREX_SPACEDIM;
+ IntVect is_nodal[AMREX_SPACEDIM];
+ for (int d=0; d < AMREX_SPACEDIM; ++d) {
+ is_nodal[d] = amrex::IntVect::TheZeroVector();
+ is_nodal[d][d] = 1;
+ }
+ linear_interpolate_to_particle_z(p, plo, dxi, data_arr.data(), height_arr,
+ val, &is_nodal[0], icomp, ncomp_per_array, num_arrays);
+}
+
+/**
+ \brief Linearly interpolates the mesh data to the particle position from mesh data.
+ This general form can handle an arbitrary number of Array4s, each with different staggerings
+ on a terrain-fitted grid.
+*/
+template <typename P>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void linear_interpolate_to_particle_z (const P& p,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+ const Array4<amrex::Real const>* data_arr,
+ const amrex::Array4<amrex::Real const>& height_arr,
+ amrex::ParticleReal * val,
+ const IntVect* is_nodal,
+ int start_comp, int ncomp, int num_arrays)
+{
+#if (AMREX_SPACEDIM == 1)
+ amrex::ignore_unused(p, plo, dxi, data_arr, height_arr, val, is_nodal, start_comp, ncomp, num_arrays);
+ amrex::Abort(" Terrain fitted grid interpolation is not supported in 1D\n");
+#else
+ AMREX_ASSERT(val != nullptr);
+
+ int ctr = 0;
+
+ for (int d = 0; d < num_arrays; d++)
+ {
+ AMREX_D_TERM(amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast<Real>(!is_nodal[d][0])*Real(0.5);,
+ amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast<Real>(!is_nodal[d][1])*Real(0.5);,);
+
+ int const i0 = static_cast<int>(amrex::Math::floor(lx));
+ int k = 0; // temporary for 2D
+
+#if (AMREX_SPACEDIM == 2)
+ amrex::ignore_unused(ly);
+ int const j = p.idata(0);
+
+ amrex::Real hlo_xlo = amrex::Real(0.25) * ( height_arr(i0 , j , k)
+ + height_arr(i0 + (!(is_nodal[d][0])) , j , k)
+ + height_arr(i0 , j + (!is_nodal[d][1]) , k)
+ + height_arr(i0 + (!(is_nodal[d][0])) , j + (!is_nodal[d][1]) , k) );
+
+ amrex::Real hlo_xhi = 0.25 * ( height_arr(i0 + 1 , j , k )
+ + height_arr(i0 + 1 + (!(is_nodal[d][0])) , j , k )
+ + height_arr(i0 + 1 , j + (!is_nodal[d][1]), k )
+ + height_arr(i0 + 1 + (!(is_nodal[d][0])) , j + (!is_nodal[d][1]), k ) );
+
+
+ amrex::Real const xint = lx - static_cast<Real>(i0);
+ amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
+ amrex::Real height_at_px = sx[0] * hlo_xlo + sx[1] * hlo_xhi;
+
+ int const j0 = (amrex::Real(p.pos(1)) >= height_at_px) ? j : j-1;
+
+ int yctr = 0;
+ amrex::Real ht[4];
+ for (int ii=0; ii < 2; ++ii) {
+ for (int jj=0; jj < 2; ++jj) {
+ ht[yctr] = amrex::Real(0.25) * ( height_arr(i0 + ii , j0 + jj , k )
+ + height_arr(i0 + ii + (!(is_nodal[d][0])) , j0 + jj , k )
+ + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k )
+ + height_arr(i0 + ii + (!(is_nodal[d][0])) , j0 + jj + (!is_nodal[d][1]), k ) );
+ ++yctr;
+ }
+ }
+ amrex::Real hint_ilo = (p.pos(1) - ht[0]) / (ht[1] - ht[0]);
+ amrex::Real hint_ihi = (p.pos(1) - ht[2]) / (ht[3] - ht[2]);
+
+ amrex::Real sy[] = { amrex::Real(1.) - hint_ilo, amrex::Real(1.) - hint_ihi,
+ hint_ilo, hint_ihi};
+
+#elif (AMREX_SPACEDIM == 3)
+
+ int const j0 = static_cast<int>(amrex::Math::floor(ly));
+ k = p.idata(0);
+ amrex::Real const xint = lx - static_cast<Real>(i0);
+ amrex::Real const yint = ly - static_cast<Real>(j0);
+ amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
+ amrex::Real sy[] = { amrex::Real(1.) - yint, yint};
+
+ amrex::Real hlo[4];
+ int ilo = 0;
+ amrex::Real height_at_pxy = 0.;
+ for (int ii = 0; ii < 2; ++ii) {
+ for (int jj = 0; jj < 2; ++jj) {
+ hlo[ilo] = amrex::Real(0.125)
+ * ( height_arr(i0 + ii , j0 + jj , k )
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k )
+ + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k )
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k )
+ + height_arr(i0 + ii , j0 + jj , k + (!is_nodal[d][2]))
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k + (!is_nodal[d][2]))
+ + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k + (!is_nodal[d][2]))
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k + (!is_nodal[d][2]))
+ );
+ height_at_pxy += hlo[ilo] * sx[ii] * sy[jj];
+ ++ilo;
+ }
+ }
+
+ int const k0 = (amrex::Real(p.pos(2)) >= height_at_pxy ) ? k : k-1;
+
+ // to make z0-z7 concise
+ int zctr = 0;
+ amrex::Real ht[8];
+ for (int ii = 0; ii < 2; ++ii) {
+ for (int jj = 0; jj < 2; ++jj) {
+ for (int kk = 0; kk < 2; ++kk) {
+ ht[zctr] = amrex::Real(0.125) *
+ ( height_arr(i0 + ii , j0 + jj , k0 + kk )
+ + height_arr(i0 + ii , j0 + jj , k0 + kk + (!is_nodal[d][2]))
+ + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k0 + kk )
+ + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k0 + kk + (!is_nodal[d][2]))
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k0 + kk )
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k0 + kk + (!is_nodal[d][2]))
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k0 + kk )
+ + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k0 + kk + (!is_nodal[d][2]))
+ );
+ ++zctr;
+ }}}
+
+ amrex::Real hint_ilojlo = ( p.pos(2) - ht[0] ) / (ht[1] - ht[0]);
+ amrex::Real hint_ilojhi = ( p.pos(2) - ht[2] ) / (ht[3] - ht[2]);
+ amrex::Real hint_ihijlo = ( p.pos(2) - ht[4] ) / (ht[5] - ht[4]);
+ amrex::Real hint_ihijhi = ( p.pos(2) - ht[6] ) / (ht[7] - ht[6]);
+
+ amrex::Real sz[] = { amrex::Real(1.) - hint_ilojlo,
+ amrex::Real(1.) - hint_ihijlo,
+ amrex::Real(1.) - hint_ilojhi,
+ amrex::Real(1.) - hint_ihijhi,
+ hint_ilojlo,
+ hint_ihijlo,
+ hint_ilojhi,
+ hint_ihijhi};
+#endif
+ for (int comp = start_comp; comp < ncomp; ++comp) {
+ val[ctr] = amrex::ParticleReal(0.);
+#if (AMREX_SPACEDIM == 2)
+ // interpolate 4 points in y
+ int k0 = 0;
+ int sy_ctr = 0;
+ for (int jj = 0; jj <= 1; ++jj) {
+ for (int ii = 0; ii <=1; ++ii) {
+ val[ctr] += static_cast<ParticleReal>( (data_arr[d])(i0+ii, j0+jj, k0 ,comp)*sx[ii]*sy[sy_ctr] );
+ ++sy_ctr;
+ }
+ }
+#elif (AMREX_SPACEDIM == 3)
+ // not tested yet
+ int sz_ctr = 0;
+ for (int kk = 0; kk <= 1; ++kk) {
+ for (int jj = 0; jj <= 1; ++jj) {
+ for (int ii = 0; ii <= 1; ++ii) {
+ val[ctr] += static_cast<ParticleReal>(
+ (data_arr[d])(i0+ii, j0+jj, k0 + kk, comp)*sx[ii]*sy[jj]*sz[sz_ctr]);
+ ++sz_ctr;
+ }
+ }
+ }
+#endif
+ ctr++;
+ } // ncomp
+ } // d
+#endif
+}
+
+
 } // namespace amrex
 #endif // include guard