Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_TracerParticle_mod_K.H
Go to the documentation of this file.
1#ifndef AMREX_TRACERPARTICLE_MOD_K_H
2#define AMREX_TRACERPARTICLE_MOD_K_H
3
4#include <AMReX_Config.H>
5#include <AMReX_FArrayBox.H>
6#include <AMReX_Box.H>
7#include <AMReX_Gpu.H>
8#include <AMReX_Geometry.H>
9#include <AMReX_REAL.H>
10#include <AMReX_IntVect.H>
11#include <AMReX_TracerParticles.H>
12#include <cmath>
13
14namespace amrex {
15
16//
17// **********************************************************************
18// Regular coordinates
19// **********************************************************************
20//
21
25template <typename P>
26AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
27void cic_interpolate (const P& p,
28 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
29 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
30 const amrex::Array4<amrex::Real const>& data_arr,
31 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
32{
33 cic_interpolate_cc(p, plo, dxi, data_arr, val, M);
34}
35
36template <typename P>
37AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
38void cic_interpolate_cc (const P& p,
39 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
40 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
41 const amrex::Array4<amrex::Real const>& data_arr,
42 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
43{
44 int start_comp = 0;
45 int ncomp_per_array = M;
46 int num_arrays = 1;
47 IntVect is_nodal = amrex::IntVect::TheZeroVector();
48 linear_interpolate_to_particle (p, plo, dxi, &data_arr, val, &is_nodal, start_comp, ncomp_per_array, num_arrays);
49}
50
54template <typename P>
55AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
56void cic_interpolate_nd (const P& p,
57 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
58 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
59 const amrex::Array4<amrex::Real const>& data_arr,
60 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
61{
62 int start_comp = 0;
63 int ncomp_per_array = M;
64 int num_arrays = 1;
65 IntVect is_nodal = amrex::IntVect::TheUnitVector();
66 linear_interpolate_to_particle (p, plo, dxi, &data_arr, val, &is_nodal, start_comp, ncomp_per_array, num_arrays);
67}
68
73template <typename P>
74AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
75void mac_interpolate (const P& p,
76 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
77 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
78 amrex::GpuArray<amrex::Array4<amrex::Real const>,AMREX_SPACEDIM> const& data_arr,
79 amrex::ParticleReal * val)
80{
81 int start_comp = 0;
82 int ncomp_per_array = 1;
83 int num_arrays = AMREX_SPACEDIM;
84 IntVect is_nodal[AMREX_SPACEDIM];
85 for (int d=0; d < AMREX_SPACEDIM; ++d) {
86 is_nodal[d] = amrex::IntVect::TheZeroVector();
87 is_nodal[d][d] = 1;
88 }
89 linear_interpolate_to_particle (p, plo, dxi, data_arr.data(), val, &is_nodal[0], start_comp, ncomp_per_array, num_arrays);
90}
91
92
93
98template <typename P>
99AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
100void linear_interpolate_to_particle (const P& p,
101 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
102 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
103 const Array4<amrex::Real const>* data_arr,
104 amrex::ParticleReal * val,
105 const IntVect* is_nodal,
106 int start_comp, int ncomp, int num_arrays)
107{
108 AMREX_ASSERT(val != nullptr);
109
110 int ctr = 0;
111
112 for (int d = 0; d < num_arrays; d++)
113 {
114 AMREX_D_TERM(amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast<Real>(!is_nodal[d][0])*Real(0.5);,
115 amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast<Real>(!is_nodal[d][1])*Real(0.5);,
116 amrex::Real lz = (Real(p.pos(2))-plo[2])*dxi[2] - static_cast<Real>(!is_nodal[d][2])*Real(0.5));
117
118 // (i0,j0,k0) is the lower corner of the box needed for interpolation
119 // i0 = (i-1) if particle is lower than center of cell i
120 // i0 = (i ) if particle is higher than center of cell i
121 AMREX_D_TERM(int const i0 = static_cast<int>(amrex::Math::floor(lx));,
122 int const j0 = static_cast<int>(amrex::Math::floor(ly));,
123 int const k0 = static_cast<int>(amrex::Math::floor(lz)));
124
125 AMREX_D_TERM(amrex::Real const xint = lx - static_cast<Real>(i0);,
126 amrex::Real const yint = ly - static_cast<Real>(j0);,
127 amrex::Real const zint = lz - static_cast<Real>(k0));
128
129 amrex::Real sx[] = {amrex::Real(1.0) - xint, xint};
130#if (AMREX_SPACEDIM > 1)
131 amrex::Real sy[] = {amrex::Real(1.0) - yint, yint};
132#endif
133#if (AMREX_SPACEDIM > 2)
134 amrex::Real sz[] = {amrex::Real(1.0) - zint, zint};
135#endif
136
137 for (int comp = start_comp; comp < start_comp + ncomp; ++comp) {
138 val[ctr] = ParticleReal(0.0);
139#if (AMREX_SPACEDIM > 2)
140 for (int kk = 0; kk <=1; ++kk) {
141#endif
142
143#if (AMREX_SPACEDIM > 1)
144 for (int jj = 0; jj <= 1; ++jj) {
145#endif
146 for (int ii = 0; ii <= 1; ++ii) {
147 val[ctr] += static_cast<ParticleReal>((data_arr[d])(IntVect(AMREX_D_DECL(i0+ii, j0+jj, k0+kk)), comp) *
148 AMREX_D_TERM(sx[ii],*sy[jj],*sz[kk]));
149 AMREX_D_TERM(},},});
150 ctr++;
151 } // ncomp
152 } // d
153}
154
155//
156// **********************************************************************
157// Terrain-fitted coordinates
158// **********************************************************************
159//
160
165template <typename P>
166AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
167void cic_interpolate_mapped_z (const P& p,
168 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
169 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
170 const amrex::Array4<amrex::Real const>& data_arr,
171 const amrex::Array4<amrex::Real const>& height_arr,
172 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
173{
174 cic_interpolate_cc_mapped_z(p, plo, dxi, data_arr, height_arr, val, M);
175}
176
181template <typename P>
182AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
183void cic_interpolate_cc_mapped_z (const P& p,
184 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
185 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
186 const amrex::Array4<amrex::Real const>& data_arr,
187 const amrex::Array4<amrex::Real const>& height_arr,
188 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
189{
190 int icomp = 0;
191 int ncomp_per_array = M;
192 int num_arrays = 1;
193 IntVect is_nodal = amrex::IntVect::TheZeroVector();
194 linear_interpolate_to_particle_z(p, plo, dxi, &data_arr, height_arr,
195 val, &is_nodal, icomp, ncomp_per_array, num_arrays);
196}
197
202template <typename P>
203AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
204void cic_interpolate_nd_mapped_z (const P& p,
205 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
206 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
207 const amrex::Array4<amrex::Real const>& data_arr,
208 const amrex::Array4<amrex::Real const>& height_arr,
209 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
210{
211 int icomp = 0;
212 int ncomp_per_array = M;
213 int num_arrays = 1;
214 IntVect is_nodal = amrex::IntVect::TheUnitVector();
215 linear_interpolate_to_particle_z(p, plo, dxi, &data_arr, height_arr,
216 val, &is_nodal, icomp, ncomp_per_array, num_arrays);
217}
218
224template <typename P>
225AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
226void mac_interpolate_mapped_z (const P& p,
227 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
228 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
229 amrex::GpuArray<amrex::Array4<amrex::Real const>,AMREX_SPACEDIM> const& data_arr,
230 const amrex::Array4<amrex::Real const>& height_arr,
231 amrex::ParticleReal * val)
232{
233 int icomp = 0;
234 int ncomp_per_array = 1;
235 int num_arrays = AMREX_SPACEDIM;
236 IntVect is_nodal[AMREX_SPACEDIM];
237 for (int d=0; d < AMREX_SPACEDIM; ++d) {
238 is_nodal[d] = amrex::IntVect::TheZeroVector();
239 is_nodal[d][d] = 1;
240 }
241 linear_interpolate_to_particle_z(p, plo, dxi, data_arr.data(), height_arr,
242 val, &is_nodal[0], icomp, ncomp_per_array, num_arrays);
243}
244
250template <typename P>
251AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
252void linear_interpolate_to_particle_z (const P& p,
253 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
254 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
255 const Array4<amrex::Real const>* data_arr,
256 const amrex::Array4<amrex::Real const>& height_arr,
257 amrex::ParticleReal * val,
258 const IntVect* is_nodal,
259 int start_comp, int ncomp, int num_arrays)
260{
261#if (AMREX_SPACEDIM == 1)
262 amrex::ignore_unused(p, plo, dxi, data_arr, height_arr, val, is_nodal, start_comp, ncomp, num_arrays);
263 amrex::Abort(" Terrain fitted grid interpolation is not supported in 1D\n");
264#else
265 AMREX_ASSERT(val != nullptr);
266
267 int ctr = 0;
268
269 for (int d = 0; d < num_arrays; d++)
270 {
271 AMREX_D_TERM(amrex::Real lx = (Real(p.pos(0))-plo[0])*dxi[0] - static_cast<Real>(!is_nodal[d][0])*Real(0.5);,
272 amrex::Real ly = (Real(p.pos(1))-plo[1])*dxi[1] - static_cast<Real>(!is_nodal[d][1])*Real(0.5);,);
273
274 int const i0 = static_cast<int>(amrex::Math::floor(lx));
275 int k = 0;
276
277#if (AMREX_SPACEDIM == 2)
278 amrex::ignore_unused(ly);
279 int const j = p.idata(0);
280
281 int j0 = j;
282
283 amrex::Real const xint = lx - static_cast<Real>(i0);
284 amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
285
286 // hlo holds the height of the face-based velocity at the relevant face center, e.g.
287 // for u it holds the height of the x-face centers of the two cells that will contribute
288 // for v it holds the height of the y-face centers of the two cells that will contribute
289 //
290 // height_at_px is then the value of hlo interpolated to the (i) of the particle location
291 //
292 // We only use height_at_px to determine which cells (in the j-direction) we should use for interpolation
293 // -- but ONLY for the lateral velocities because otherwise we risk using vertical velocity below the sfc
294 //
295 if (d != AMREX_SPACEDIM-1) {
296 amrex::Real hlo_xlo = amrex::Real(0.25)
297 * ( height_arr(i0 , j , k)
298 + height_arr(i0 + (!(is_nodal[d][0])) , j , k)
299 + height_arr(i0 , j + (!is_nodal[d][1]) , k)
300 + height_arr(i0 + (!(is_nodal[d][0])) , j + (!is_nodal[d][1]) , k) );
301
302 amrex::Real hlo_xhi = amrex::Real(0.25)
303 * ( height_arr(i0 + 1 , j , k )
304 + height_arr(i0 + 1 + (!(is_nodal[d][0])) , j , k )
305 + height_arr(i0 + 1 , j + (!is_nodal[d][1]), k )
306 + height_arr(i0 + 1 + (!(is_nodal[d][0])) , j + (!is_nodal[d][1]), k ) );
307
308 amrex::Real height_at_px = sx[0] * hlo_xlo + sx[1] * hlo_xhi;
309
310 if (amrex::Real(p.pos(1)) < height_at_px) {j0 = j-1;}
311
312 } // if not vertical direction
313
314 // ht holds the heights of the face-based velocity at the relevant face centers, e.g.
315 // for u it holds the height of the x-face centers of the four cells that will contribute
316 // for v it holds the height of the y-face centers of the four cells that will contribute
317 //
318 // ht differs from hlo in that it uses the correct j0 which may be different from k above
319 // AND it computes all 4 values needed for 2D interpolation (hlo held only the lower 2)
320 //
321 int yctr = 0;
322 amrex::Real ht[4];
323 for (int ii=0; ii < 2; ++ii) {
324 for (int jj=0; jj < 2; ++jj) {
325 ht[yctr] = amrex::Real(0.25)
326 * ( height_arr(i0 + ii , j0 + jj , k )
327 + height_arr(i0 + ii + (!(is_nodal[d][0])) , j0 + jj , k )
328 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k )
329 + height_arr(i0 + ii + (!(is_nodal[d][0])) , j0 + jj + (!is_nodal[d][1]), k ) );
330 ++yctr;
331 }
332 }
333 amrex::Real hint_ilo = (p.pos(1) - ht[0]) / (ht[1] - ht[0]);
334 amrex::Real hint_ihi = (p.pos(1) - ht[2]) / (ht[3] - ht[2]);
335
336 amrex::Real sy[] = { amrex::Real(1.) - hint_ilo, amrex::Real(1.) - hint_ihi,
337 hint_ilo, hint_ihi};
338
339#elif (AMREX_SPACEDIM == 3)
340
341 int const j0 = static_cast<int>(amrex::Math::floor(ly));
342 k = p.idata(0);
343 amrex::Real const xint = lx - static_cast<Real>(i0);
344 amrex::Real const yint = ly - static_cast<Real>(j0);
345 amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
346 amrex::Real sy[] = { amrex::Real(1.) - yint, yint};
347
348 int k0 = k;
349
350 // hlo holds the height of the face-based velocity at the relevant face center, e.g.
351 // for u it holds the height of the x-face centers of the four cells that will contribute
352 // for v it holds the height of the y-face centers of the four cells that will contribute
353 // for w it holds the height of the z-face centers of the four cells that will contribute
354 //
355 // height_at_pxy is then the value of hlo interpolated to the (i,j) of the particle location
356 //
357 // We use height_at_pxy to determine which cells (in the k-direction) we should use for interpolation
358 // -- but ONLY for the lateral velocities because otherwise we risk using vertical velocity below the sfc
359 //
360
361 if (d != AMREX_SPACEDIM-1) {
362 amrex::Real hlo[4];
363 int ilo = 0;
364 amrex::Real height_at_pxy = 0.;
365 for (int ii = 0; ii < 2; ++ii) {
366 for (int jj = 0; jj < 2; ++jj) {
367 hlo[ilo] = amrex::Real(0.125)
368 * ( height_arr(i0 + ii , j0 + jj , k )
369 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k )
370 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k )
371 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k )
372 + height_arr(i0 + ii , j0 + jj , k + (!is_nodal[d][2]))
373 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k + (!is_nodal[d][2]))
374 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k + (!is_nodal[d][2]))
375 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k + (!is_nodal[d][2]))
376 );
377 height_at_pxy += hlo[ilo] * sx[ii] * sy[jj];
378 ++ilo;
379 }
380 }
381
382 if (amrex::Real(p.pos(2)) < height_at_pxy ) {k0 = k-1;}
383
384 } // if not vertical direction
385
386 // ht holds the heights of the face-based velocity at the relevant face centers, e.g.
387 // for u it holds the height of the x-face centers of the eight cells that will contribute
388 // for v it holds the height of the y-face centers of the eight cells that will contribute
389 // for w it holds the height of the z-face centers of the eight cells that will contribute
390 //
391 // ht differs from hlo in that it uses the correct k0 which may be different from k above
392 // AND it computes all 8 values needed for 3D interpolation (hlo held only the lower 4)
393 //
394 int zctr = 0;
395 amrex::Real ht[8];
396 for (int ii = 0; ii < 2; ++ii) {
397 for (int jj = 0; jj < 2; ++jj) {
398 for (int kk = 0; kk < 2; ++kk) {
399 ht[zctr] = amrex::Real(0.125) *
400 ( height_arr(i0 + ii , j0 + jj , k0 + kk )
401 + height_arr(i0 + ii , j0 + jj , k0 + kk + (!is_nodal[d][2]))
402 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k0 + kk )
403 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k0 + kk + (!is_nodal[d][2]))
404 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k0 + kk )
405 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k0 + kk + (!is_nodal[d][2]))
406 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k0 + kk )
407 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k0 + kk + (!is_nodal[d][2]))
408 );
409 ++zctr;
410 }}}
411
412 amrex::Real hint_ilojlo = ( p.pos(2) - ht[0] ) / (ht[1] - ht[0]);
413 amrex::Real hint_ilojhi = ( p.pos(2) - ht[2] ) / (ht[3] - ht[2]);
414 amrex::Real hint_ihijlo = ( p.pos(2) - ht[4] ) / (ht[5] - ht[4]);
415 amrex::Real hint_ihijhi = ( p.pos(2) - ht[6] ) / (ht[7] - ht[6]);
416
417 amrex::Real sz[] = { amrex::Real(1.) - hint_ilojlo,
418 amrex::Real(1.) - hint_ihijlo,
419 amrex::Real(1.) - hint_ilojhi,
420 amrex::Real(1.) - hint_ihijhi,
421 hint_ilojlo,
422 hint_ihijlo,
423 hint_ilojhi,
424 hint_ihijhi};
425#endif
426 for (int comp = start_comp; comp < start_comp + ncomp; ++comp) {
427 val[ctr] = amrex::ParticleReal(0.);
428#if (AMREX_SPACEDIM == 2)
429 int k0 = 0;
430 int sy_ctr = 0;
431 for (int jj = 0; jj <= 1; ++jj) {
432 for (int ii = 0; ii <=1; ++ii) {
433 val[ctr] += static_cast<ParticleReal>( (data_arr[d])(i0+ii, j0+jj, k0 ,comp)*sx[ii]*sy[sy_ctr] );
434 ++sy_ctr;
435 }
436 }
437#elif (AMREX_SPACEDIM == 3)
438 int sz_ctr = 0;
439 for (int kk = 0; kk <= 1; ++kk) {
440 for (int jj = 0; jj <= 1; ++jj) {
441 for (int ii = 0; ii <= 1; ++ii) {
442 val[ctr] += static_cast<ParticleReal>(
443 (data_arr[d])(i0+ii, j0+jj, k0 + kk, comp)*sx[ii]*sy[jj]*sz[sz_ctr]);
444 ++sz_ctr;
445 }
446 }
447 }
448#endif
449 ctr++;
450 } // ncomp
451 } // d
452#endif
453}
454
455
456//
457// **********************************************************************
458// General mapped coordinates
459// **********************************************************************
460//
461
466AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
467void particle_interp_decomp (Array2D<Real,1,7,1,7>& a, Array1D<int,1,7>& ip, int neq)
468{
469 amrex::Real amult,apm,apn;
470 for (int n=1; n<neq+1; n++)
471 {
472 ip(n) = n;
473 }
474
475 for (int i=1; i<neq; i++)
476 {
477 int ip1=i+1;
478 int k=i;
479 int ipi = ip(i);
480 apm = std::abs(a(ipi,i));
481
482 for (int j=ip1; j<neq+1; j++)
483 {
484 int ipj = ip(j);
485 apn=std::abs(a(ipj,i));
486 if(apm < apn){
487 apm = apn;
488 k=j;
489 }
490 }
491
492 int j=ip(k);
493 ip(k)=ip(i);
494 ip(i)=j;
495 for (int l=ip1; l<neq+1; l++)
496 {
497 int n=ip(l);
498 amult=a(n,i)/a(j,i);
499 a(n,i)=amult;
500 for (int kk = ip1; kk < neq+1; kk++)
501 {
502 a(n,kk) -= amult*a(j,kk);
503 }
504 }
505 }
506}
507
508AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
509void particle_interp_solve (Array2D<Real,1,7,1,7>& a, Array1D<Real,1,7>& b, Array1D<int,1,7>& ip, int neq)
510{
511 amrex::Real scr;
512 for (int l=2; l<neq+1; l++)
513 {
514 int n = ip(l);
515
516 for (int k=1; k<l; k++){
517 b(n) -= a(n,k) * b(ip(k));
518 }
519 }
520 b(ip(neq))=b(ip(neq))/a(ip(neq),neq);
521
522 for (int l=1; l<neq; l++)
523 {
524 int j=neq-l;
525 int jp1=j+1;
526 int n=ip(j);
527 for (int k=jp1; k<neq+1; k++){
528 b(n) -= a(n,k ) * b(ip(k));
529 }
530 b(n)=b(n)/a(n,j);
531 }
532 for (int n=1; n<neq+1; n++)
533 {
534 while(ip(n) != n) {
535 int j=ip(n);
536 scr=b(j);
537 ip(n)=ip(j);
538 b(j)=b(ip(j));
539 b(ip(j))=scr;
540 ip(j)=j;
541 }
542 }
543}
544
545
549template <typename P>
550AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
551void cic_interpolate_nd_mapped (const P& p,
552 const amrex::Array4<amrex::Real const>& data_arr,
553 const amrex::Array4<amrex::Real const>& loc_arr,
554 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
555{
556 int icomp = 0;
557 int ncomp_per_array = M;
558 int num_arrays = 1;
559 IntVect is_nodal = amrex::IntVect::TheUnitVector();
560 linear_interpolate_to_particle_mapped(p, &data_arr, loc_arr,
561 val, &is_nodal, icomp, ncomp_per_array, num_arrays);
562}
563
569template <typename P>
570AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
571void linear_interpolate_to_particle_mapped (const P& p,
572 const Array4<amrex::Real const>* data_arr,
573 const amrex::Array4<amrex::Real const>& loc_arr,
574 amrex::ParticleReal * val,
575 const IntVect* is_nodal,
576 int start_comp, int ncomp, int num_arrays)
577{
578#if (AMREX_SPACEDIM == 1)
579 amrex::ignore_unused(p, data_arr, loc_arr, val, is_nodal, start_comp, ncomp, num_arrays);
580 amrex::Abort(" General mapped interpolation is not supported in 1D\n");
581#else
582
583 AMREX_ASSERT(val != nullptr);
584 AMREX_ASSERT(num_arrays == 1);
585 for (int d = 0; d < num_arrays; d++)
586 {
587 // Assertion to ensure that only nodal data is supported
588 AMREX_ALWAYS_ASSERT_WITH_MESSAGE((is_nodal[d][0] == 1),"For general mapped coordinates, interpolation is supported only for node-centered data");
589 AMREX_ALWAYS_ASSERT_WITH_MESSAGE((is_nodal[d][1] == 1),"For general mapped coordinates, interpolation is supported only for node-centered data");
590#if (AMREX_SPACEDIM == 3)
591 AMREX_ALWAYS_ASSERT_WITH_MESSAGE((is_nodal[d][2] == 1),"For general mapped coordinates, interpolation is supported only for node-centered data");
592#endif
593 }
594 int n_of_data = 0;
595
596 int i = p.idata(0);
597 int j = p.idata(1);
598 for (int comp = start_comp; comp < start_comp + ncomp; ++comp) {
599
600#if (AMREX_SPACEDIM == 2)
601 // Value of data at surrounding nodes
602 amrex::Real f1 = (data_arr[n_of_data])(i ,j ,0,comp);
603 amrex::Real f2 = (data_arr[n_of_data])(i ,j+1,0,comp) - f1;
604 amrex::Real f3 = (data_arr[n_of_data])(i+1,j ,0,comp) - f1;
605 amrex::Real f4 = (data_arr[n_of_data])(i+1,j+1,0,comp) - f1;
606
607 // Node locations
608 amrex::Real x1 = loc_arr(i ,j ,0,0);
609 amrex::Real x2 = loc_arr(i ,j+1,0,0) - x1;
610 amrex::Real x3 = loc_arr(i+1,j ,0,0) - x1;
611 amrex::Real x4 = loc_arr(i+1,j+1,0,0) - x1;
612
613 amrex::Real y1 = loc_arr(i ,j ,0,1);
614 amrex::Real y2 = loc_arr(i ,j+1,0,1) - y1;
615 amrex::Real y3 = loc_arr(i+1,j ,0,1) - y1;
616 amrex::Real y4 = loc_arr(i+1,j+1,0,1) - y1;
617
618 // Particle location
619 amrex::Real x = p.pos(0) - x1;
620 amrex::Real y = p.pos(1) - y1;
621
622 amrex::Real det = x2*x4*y3*(y2 - y4) - x3*x4*y2*(y3 - y4) - x2*x3*(y2 - y3)*y4;
623
624 amrex::Real b = ( f4*(x2*y2*y3 - x3*y2*y3) + (f2*(x3 - x4)*y3)*y4 + f3*(-(x2*y2*y4) + x4*y2*y4) ) / det;
625 amrex::Real c = ( -f2*x3*x4*y3 + f4*x2*x3*(-y2 + y3) + f3*x2*x4*(y2 - y4) + f2*x3*x4*y4 ) / det;
626 amrex::Real d = ( f2*x4*y3 + f4*(x3*y2 - x2*y3) - f2*x3*y4 + f3*(-(x4*y2) + x2*y4) ) / det;
627
628 amrex::Real f = b*x + c*y + d*x*y;
629 val[comp-start_comp] = f1 + f;
630
631#elif (AMREX_SPACEDIM == 3)
632
633 int k = p.idata(2);
634
635 //
636 // Ordering of particles (lo,lo,lo) (hi,lo,lo) (lo,hi,lo) (hi,hi,lo) (lo,lo,hi) (hi,lo,hi) (lo,hi,hi) (hi,hi,hi)
637 // All locations are calculated as offsets from (lo,lo,lo), which is only need 7 values are needed
638 //
639 Array1D<Real,0,7> x_vals;
640 Array1D<Real,0,7> y_vals;
641 Array1D<Real,0,7> z_vals;
642 Array1D<Real,1,7> b;
643
644 amrex::Real x1, y1, z1, f1;
645
646 x1 = loc_arr(i ,j ,k ,0);
647 x_vals(1) = loc_arr(i+1,j ,k ,0) - x1;
648 x_vals(2) = loc_arr(i ,j+1,k ,0) - x1;
649 x_vals(3) = loc_arr(i+1,j+1,k ,0) - x1;
650 x_vals(4) = loc_arr(i ,j ,k+1,0) - x1;
651 x_vals(5) = loc_arr(i+1,j ,k+1,0) - x1;
652 x_vals(6) = loc_arr(i ,j+1,k+1,0) - x1;
653 x_vals(7) = loc_arr(i+1,j+1,k+1,0) - x1;
654
655 y1 = loc_arr(i ,j ,k ,1);
656 y_vals(1) = loc_arr(i+1,j ,k ,1) - y1;
657 y_vals(2) = loc_arr(i ,j+1,k ,1) - y1;
658 y_vals(3) = loc_arr(i+1,j+1,k ,1) - y1;
659 y_vals(4) = loc_arr(i ,j ,k+1,1) - y1;
660 y_vals(5) = loc_arr(i+1,j ,k+1,1) - y1;
661 y_vals(6) = loc_arr(i ,j+1,k+1,1) - y1;
662 y_vals(7) = loc_arr(i+1,j+1,k+1,1) - y1;
663
664 z1 = loc_arr(i ,j ,k ,2);
665 z_vals(1) = loc_arr(i+1,j ,k ,2) - z1;
666 z_vals(2) = loc_arr(i ,j+1,k ,2) - z1;
667 z_vals(3) = loc_arr(i+1,j+1,k ,2) - z1;
668 z_vals(4) = loc_arr(i ,j ,k+1,2) - z1;
669 z_vals(5) = loc_arr(i+1,j ,k+1,2) - z1;
670 z_vals(6) = loc_arr(i ,j+1,k+1,2) - z1;
671 z_vals(7) = loc_arr(i+1,j+1,k+1,2) - z1;
672
673 f1 = (data_arr[n_of_data])(i ,j ,k ,comp);
674 b(1) = (data_arr[n_of_data])(i+1,j ,k ,comp) - f1;
675 b(2) = (data_arr[n_of_data])(i ,j+1,k ,comp) - f1;
676 b(3) = (data_arr[n_of_data])(i+1,j+1,k ,comp) - f1;
677 b(4) = (data_arr[n_of_data])(i ,j ,k+1,comp) - f1;
678 b(5) = (data_arr[n_of_data])(i+1,j ,k+1,comp) - f1;
679 b(6) = (data_arr[n_of_data])(i ,j+1,k+1,comp) - f1;
680 b(7) = (data_arr[n_of_data])(i+1,j+1,k+1,comp) - f1;
681
682 Array2D<Real,1,7,1,7> a;
683 Array1D<int,1,7> ip;
684
685 int neq = 7;
686
687 for (int n=1; n<neq+1; n++) {
688 a(n,1) = x_vals(n);
689 a(n,2) = y_vals(n);
690 a(n,3) = z_vals(n);
691 a(n,4) = x_vals(n)*y_vals(n);
692 a(n,5) = x_vals(n)*z_vals(n);
693 a(n,6) = z_vals(n)*y_vals(n);
694 a(n,7) = x_vals(n)*y_vals(n)*z_vals(n);
695 }
696
697 particle_interp_decomp(a,ip,neq);
698 particle_interp_solve(a,b,ip,neq);
699
700 amrex::Real px, py, pz;
701 px = p.pos(0) - x1;
702 py = p.pos(1) - y1;
703 pz = p.pos(2) - z1;
704
705 amrex::Real f = b(1)*px + b(2)*py + b(3)*pz + b(4)*px*py + b(5)*px*pz + b(6)*py*pz + b(7)*px*py*pz;
706 val[comp-start_comp] = static_cast<ParticleReal>(f + f1);
707#endif
708
709 } // comp
710#endif
711}
712
713} // namespace amrex
714#endif // include guard
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
AMREX_ASSERT
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
AMReX_Box.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMReX_FArrayBox.H
AMReX_Geometry.H
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
AMReX_Gpu.H
AMReX_IntVect.H
AMReX_REAL.H
AMREX_D_TERM
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:172
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
Definition AMReX_SPACE.H:171
AMReX_TracerParticles.H
amrex::IntVectND< 3 >::TheUnitVector
__host__ static __device__ constexpr IntVectND< dim > TheUnitVector() noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:687
amrex::IntVectND< 3 >::TheZeroVector
__host__ static __device__ constexpr IntVectND< dim > TheZeroVector() noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:677
amrex
Definition AMReX_Amr.cpp:49
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:138
amrex::cic_interpolate
__host__ __device__ void cic_interpolate(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const amrex::Array4< amrex::Real const > &data_arr, amrex::ParticleReal *val, int M=3)
Linearly interpolates the mesh data to the particle position from cell-centered data.
Definition AMReX_TracerParticle_mod_K.H:27
amrex::linear_interpolate_to_particle_z
__host__ __device__ void linear_interpolate_to_particle_z(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > 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)
Linearly interpolates the mesh data to the particle position from mesh data. This general form can ha...
Definition AMReX_TracerParticle_mod_K.H:252
amrex::mac_interpolate_mapped_z
__host__ __device__ void mac_interpolate_mapped_z(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, amrex::GpuArray< amrex::Array4< amrex::Real const >, 3 > const &data_arr, const amrex::Array4< amrex::Real const > &height_arr, amrex::ParticleReal *val)
Linearly interpolates the mesh data to the particle position from face-centered data on a terrain-fit...
Definition AMReX_TracerParticle_mod_K.H:226
amrex::linear_interpolate_to_particle
__host__ __device__ void linear_interpolate_to_particle(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const Array4< amrex::Real const > *data_arr, amrex::ParticleReal *val, const IntVect *is_nodal, int start_comp, int ncomp, int num_arrays)
Linearly interpolates the mesh data to the particle position from mesh data. This general form can ha...
Definition AMReX_TracerParticle_mod_K.H:100
amrex::cic_interpolate_nd_mapped_z
__host__ __device__ void cic_interpolate_nd_mapped_z(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const amrex::Array4< amrex::Real const > &data_arr, const amrex::Array4< amrex::Real const > &height_arr, amrex::ParticleReal *val, int M=3)
Linearly interpolates the mesh data to the particle position from node-centered data....
Definition AMReX_TracerParticle_mod_K.H:204
amrex::linear_interpolate_to_particle_mapped
__host__ __device__ void linear_interpolate_to_particle_mapped(const P &p, const Array4< amrex::Real const > *data_arr, const amrex::Array4< amrex::Real const > &loc_arr, amrex::ParticleReal *val, const IntVect *is_nodal, int start_comp, int ncomp, int num_arrays)
Linearly interpolates the mesh data to the particle position on a general mapped grid....
Definition AMReX_TracerParticle_mod_K.H:571
amrex::cic_interpolate_nd_mapped
__host__ __device__ void cic_interpolate_nd_mapped(const P &p, const amrex::Array4< amrex::Real const > &data_arr, const amrex::Array4< amrex::Real const > &loc_arr, amrex::ParticleReal *val, int M=3)
Linearly interpolates the mesh data to the particle position from node-centered data on a general map...
Definition AMReX_TracerParticle_mod_K.H:551
amrex::cic_interpolate_cc_mapped_z
__host__ __device__ void cic_interpolate_cc_mapped_z(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const amrex::Array4< amrex::Real const > &data_arr, const amrex::Array4< amrex::Real const > &height_arr, amrex::ParticleReal *val, int M=3)
Linearly interpolates the mesh data to the particle position from cell-centered data on a terrain-fit...
Definition AMReX_TracerParticle_mod_K.H:183
amrex::cic_interpolate_cc
__host__ __device__ void cic_interpolate_cc(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const amrex::Array4< amrex::Real const > &data_arr, amrex::ParticleReal *val, int M=3)
Definition AMReX_TracerParticle_mod_K.H:38
amrex::IntVect
IntVectND< 3 > IntVect
Definition AMReX_BaseFwd.H:30
amrex::mac_interpolate
__host__ __device__ void mac_interpolate(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, amrex::GpuArray< amrex::Array4< amrex::Real const >, 3 > const &data_arr, amrex::ParticleReal *val)
Linearly interpolates the mesh data to the particle position from face-centered data....
Definition AMReX_TracerParticle_mod_K.H:75
amrex::particle_interp_decomp
__host__ __device__ void particle_interp_decomp(Array2D< Real, 1, 7, 1, 7 > &a, Array1D< int, 1, 7 > &ip, int neq)
Helper functions for 3D interpolation from values on nodes/vertices of arbitrary hexahedra to particl...
Definition AMReX_TracerParticle_mod_K.H:467
amrex::particle_interp_solve
__host__ __device__ void particle_interp_solve(Array2D< Real, 1, 7, 1, 7 > &a, Array1D< Real, 1, 7 > &b, Array1D< int, 1, 7 > &ip, int neq)
Definition AMReX_TracerParticle_mod_K.H:509
amrex::cic_interpolate_mapped_z
__host__ __device__ void cic_interpolate_mapped_z(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const amrex::Array4< amrex::Real const > &data_arr, const amrex::Array4< amrex::Real const > &height_arr, amrex::ParticleReal *val, int M=3)
Linearly interpolates the mesh data to the particle position from cell-centered data on a terrain-fit...
Definition AMReX_TracerParticle_mod_K.H:167
amrex::cic_interpolate_nd
__host__ __device__ void cic_interpolate_nd(const P &p, amrex::GpuArray< amrex::Real, 3 > const &plo, amrex::GpuArray< amrex::Real, 3 > const &dxi, const amrex::Array4< amrex::Real const > &data_arr, amrex::ParticleReal *val, int M=3)
Linearly interpolates the mesh data to the particle position from node-centered data.
Definition AMReX_TracerParticle_mod_K.H:56
amrex::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:230
amrex::Array1D
Definition AMReX_Array.H:193
amrex::Array2D
Definition AMReX_Array.H:341
amrex::Array4
Definition AMReX_Array4.H:61
amrex::GpuArray
Definition AMReX_Array.H:34
amrex::GpuArray::data
__host__ __device__ const T * data() const noexcept
Definition AMReX_Array.H:66