diff --git a/Src/Particle/AMReX_TracerParticle_mod_K.H b/Src/Particle/AMReX_TracerParticle_mod_K.H index 9df0864d33c..472057a408e 100644 --- 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 +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void cic_interpolate_mapped_z (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const amrex::Array4& data_arr, + const amrex::Array4& 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 +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void cic_interpolate_cc_mapped_z (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const amrex::Array4& data_arr, + const amrex::Array4& 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 +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void cic_interpolate_nd_mapped_z (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const amrex::Array4& data_arr, + const amrex::Array4& 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 +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void mac_interpolate_mapped_z (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + amrex::GpuArray,AMREX_SPACEDIM> const& data_arr, + const amrex::Array4& 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 +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void linear_interpolate_to_particle_z (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const Array4* data_arr, + const amrex::Array4& 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(!is_nodal[d][0])*Real(0.5);, + amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast(!is_nodal[d][1])*Real(0.5);,); + + int const i0 = static_cast(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(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(amrex::Math::floor(ly)); + k = p.idata(0); + amrex::Real const xint = lx - static_cast(i0); + amrex::Real const yint = ly - static_cast(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( (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( + (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