Block-Structured AMR Software Framework
 
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 amrex::Real hlo_xlo = amrex::Real(0.25)
282 * ( height_arr(i0 , j , k)
283 + height_arr(i0 + (!(is_nodal[d][0])) , j , k)
284 + height_arr(i0 , j + (!is_nodal[d][1]) , k)
285 + height_arr(i0 + (!(is_nodal[d][0])) , j + (!is_nodal[d][1]) , k) );
286
287 amrex::Real hlo_xhi = amrex::Real(0.25)
288 * ( height_arr(i0 + 1 , j , k )
289 + height_arr(i0 + 1 + (!(is_nodal[d][0])) , j , k )
290 + height_arr(i0 + 1 , j + (!is_nodal[d][1]), k )
291 + height_arr(i0 + 1 + (!(is_nodal[d][0])) , j + (!is_nodal[d][1]), k ) );
292
293
294 amrex::Real const xint = lx - static_cast<Real>(i0);
295 amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
296 amrex::Real height_at_px = sx[0] * hlo_xlo + sx[1] * hlo_xhi;
297
298 int const j0 = (amrex::Real(p.pos(1)) >= height_at_px) ? j : j-1;
299
300 int yctr = 0;
301 amrex::Real ht[4];
302 for (int ii=0; ii < 2; ++ii) {
303 for (int jj=0; jj < 2; ++jj) {
304 ht[yctr] = amrex::Real(0.25)
305 * ( height_arr(i0 + ii , j0 + jj , k )
306 + height_arr(i0 + ii + (!(is_nodal[d][0])) , j0 + jj , k )
307 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k )
308 + height_arr(i0 + ii + (!(is_nodal[d][0])) , j0 + jj + (!is_nodal[d][1]), k ) );
309 ++yctr;
310 }
311 }
312 amrex::Real hint_ilo = (p.pos(1) - ht[0]) / (ht[1] - ht[0]);
313 amrex::Real hint_ihi = (p.pos(1) - ht[2]) / (ht[3] - ht[2]);
314
315 amrex::Real sy[] = { amrex::Real(1.) - hint_ilo, amrex::Real(1.) - hint_ihi,
316 hint_ilo, hint_ihi};
317
318#elif (AMREX_SPACEDIM == 3)
319
320 int const j0 = static_cast<int>(amrex::Math::floor(ly));
321 k = p.idata(0);
322 amrex::Real const xint = lx - static_cast<Real>(i0);
323 amrex::Real const yint = ly - static_cast<Real>(j0);
324 amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
325 amrex::Real sy[] = { amrex::Real(1.) - yint, yint};
326
327 amrex::Real hlo[4];
328 int ilo = 0;
329 amrex::Real height_at_pxy = 0.;
330 for (int ii = 0; ii < 2; ++ii) {
331 for (int jj = 0; jj < 2; ++jj) {
332 hlo[ilo] = amrex::Real(0.125)
333 * ( height_arr(i0 + ii , j0 + jj , k )
334 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k )
335 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k )
336 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k )
337 + height_arr(i0 + ii , j0 + jj , k + (!is_nodal[d][2]))
338 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k + (!is_nodal[d][2]))
339 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k + (!is_nodal[d][2]))
340 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k + (!is_nodal[d][2]))
341 );
342 height_at_pxy += hlo[ilo] * sx[ii] * sy[jj];
343 ++ilo;
344 }
345 }
346
347 int const k0 = (amrex::Real(p.pos(2)) >= height_at_pxy ) ? k : k-1;
348
349 int zctr = 0;
350 amrex::Real ht[8];
351 for (int ii = 0; ii < 2; ++ii) {
352 for (int jj = 0; jj < 2; ++jj) {
353 for (int kk = 0; kk < 2; ++kk) {
354 ht[zctr] = amrex::Real(0.125) *
355 ( height_arr(i0 + ii , j0 + jj , k0 + kk )
356 + height_arr(i0 + ii , j0 + jj , k0 + kk + (!is_nodal[d][2]))
357 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k0 + kk )
358 + height_arr(i0 + ii , j0 + jj + (!is_nodal[d][1]), k0 + kk + (!is_nodal[d][2]))
359 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k0 + kk )
360 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj , k0 + kk + (!is_nodal[d][2]))
361 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k0 + kk )
362 + height_arr(i0 + ii + (!is_nodal[d][0]), j0 + jj + (!is_nodal[d][1]), k0 + kk + (!is_nodal[d][2]))
363 );
364 ++zctr;
365 }}}
366
367 amrex::Real hint_ilojlo = ( p.pos(2) - ht[0] ) / (ht[1] - ht[0]);
368 amrex::Real hint_ilojhi = ( p.pos(2) - ht[2] ) / (ht[3] - ht[2]);
369 amrex::Real hint_ihijlo = ( p.pos(2) - ht[4] ) / (ht[5] - ht[4]);
370 amrex::Real hint_ihijhi = ( p.pos(2) - ht[6] ) / (ht[7] - ht[6]);
371
372 amrex::Real sz[] = { amrex::Real(1.) - hint_ilojlo,
373 amrex::Real(1.) - hint_ihijlo,
374 amrex::Real(1.) - hint_ilojhi,
375 amrex::Real(1.) - hint_ihijhi,
376 hint_ilojlo,
377 hint_ihijlo,
378 hint_ilojhi,
379 hint_ihijhi};
380#endif
381 for (int comp = start_comp; comp < start_comp + ncomp; ++comp) {
382 val[ctr] = amrex::ParticleReal(0.);
383#if (AMREX_SPACEDIM == 2)
384 int k0 = 0;
385 int sy_ctr = 0;
386 for (int jj = 0; jj <= 1; ++jj) {
387 for (int ii = 0; ii <=1; ++ii) {
388 val[ctr] += static_cast<ParticleReal>( (data_arr[d])(i0+ii, j0+jj, k0 ,comp)*sx[ii]*sy[sy_ctr] );
389 ++sy_ctr;
390 }
391 }
392#elif (AMREX_SPACEDIM == 3)
393 int sz_ctr = 0;
394 for (int kk = 0; kk <= 1; ++kk) {
395 for (int jj = 0; jj <= 1; ++jj) {
396 for (int ii = 0; ii <= 1; ++ii) {
397 val[ctr] += static_cast<ParticleReal>(
398 (data_arr[d])(i0+ii, j0+jj, k0 + kk, comp)*sx[ii]*sy[jj]*sz[sz_ctr]);
399 ++sz_ctr;
400 }
401 }
402 }
403#endif
404 ctr++;
405 } // ncomp
406 } // d
407#endif
408}
409
410
411//
412// **********************************************************************
413// General mapped coordinates
414// **********************************************************************
415//
416
421AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
422void particle_interp_decomp (Array2D<Real,1,7,1,7>& a, Array1D<int,1,7>& ip, int neq)
423{
424 amrex::Real amult,apm,apn;
425 for (int n=1; n<neq+1; n++)
426 {
427 ip(n) = n;
428 }
429
430 for (int i=1; i<neq; i++)
431 {
432 int ip1=i+1;
433 int k=i;
434 int ipi = ip(i);
435 apm = std::abs(a(ipi,i));
436
437 for (int j=ip1; j<neq+1; j++)
438 {
439 int ipj = ip(j);
440 apn=std::abs(a(ipj,i));
441 if(apm < apn){
442 apm = apn;
443 k=j;
444 }
445 }
446
447 int j=ip(k);
448 ip(k)=ip(i);
449 ip(i)=j;
450 for (int l=ip1; l<neq+1; l++)
451 {
452 int n=ip(l);
453 amult=a(n,i)/a(j,i);
454 a(n,i)=amult;
455 for (int kk = ip1; kk < neq+1; kk++)
456 {
457 a(n,kk) -= amult*a(j,kk);
458 }
459 }
460 }
461}
462
463AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
464void particle_interp_solve (Array2D<Real,1,7,1,7>& a, Array1D<Real,1,7>& b, Array1D<int,1,7>& ip, int neq)
465{
466 amrex::Real scr;
467 for (int l=2; l<neq+1; l++)
468 {
469 int n = ip(l);
470
471 for (int k=1; k<l; k++){
472 b(n) -= a(n,k) * b(ip(k));
473 }
474 }
475 b(ip(neq))=b(ip(neq))/a(ip(neq),neq);
476
477 for (int l=1; l<neq; l++)
478 {
479 int j=neq-l;
480 int jp1=j+1;
481 int n=ip(j);
482 for (int k=jp1; k<neq+1; k++){
483 b(n) -= a(n,k ) * b(ip(k));
484 }
485 b(n)=b(n)/a(n,j);
486 }
487 for (int n=1; n<neq+1; n++)
488 {
489 while(ip(n) != n) {
490 int j=ip(n);
491 scr=b(j);
492 ip(n)=ip(j);
493 b(j)=b(ip(j));
494 b(ip(j))=scr;
495 ip(j)=j;
496 }
497 }
498}
499
500
504template <typename P>
505AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
506void cic_interpolate_nd_mapped (const P& p,
507 const amrex::Array4<amrex::Real const>& data_arr,
508 const amrex::Array4<amrex::Real const>& loc_arr,
509 amrex::ParticleReal * val, int M = AMREX_SPACEDIM)
510{
511 int icomp = 0;
512 int ncomp_per_array = M;
513 int num_arrays = 1;
514 IntVect is_nodal = amrex::IntVect::TheUnitVector();
515 linear_interpolate_to_particle_mapped(p, &data_arr, loc_arr,
516 val, &is_nodal, icomp, ncomp_per_array, num_arrays);
517}
518
524template <typename P>
525AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
526void linear_interpolate_to_particle_mapped (const P& p,
527 const Array4<amrex::Real const>* data_arr,
528 const amrex::Array4<amrex::Real const>& loc_arr,
529 amrex::ParticleReal * val,
530 const IntVect* is_nodal,
531 int start_comp, int ncomp, int num_arrays)
532{
533#if (AMREX_SPACEDIM == 1)
534 amrex::ignore_unused(p, data_arr, loc_arr, val, is_nodal, start_comp, ncomp, num_arrays);
535 amrex::Abort(" General mapped interpolation is not supported in 1D\n");
536#else
537
538 AMREX_ASSERT(val != nullptr);
539 AMREX_ASSERT(num_arrays == 1);
540 for (int d = 0; d < num_arrays; d++)
541 {
542 // Assertion to ensure that only nodal data is supported
543 AMREX_ALWAYS_ASSERT_WITH_MESSAGE((is_nodal[d][0] == 1),"For general mapped coordinates, interpolation is supported only for node-centered data");
544 AMREX_ALWAYS_ASSERT_WITH_MESSAGE((is_nodal[d][1] == 1),"For general mapped coordinates, interpolation is supported only for node-centered data");
545#if (AMREX_SPACEDIM == 3)
546 AMREX_ALWAYS_ASSERT_WITH_MESSAGE((is_nodal[d][2] == 1),"For general mapped coordinates, interpolation is supported only for node-centered data");
547#endif
548 }
549 int n_of_data = 0;
550
551 int i = p.idata(0);
552 int j = p.idata(1);
553 for (int comp = start_comp; comp < start_comp + ncomp; ++comp) {
554
555#if (AMREX_SPACEDIM == 2)
556 // Value of data at surrounding nodes
557 amrex::Real f1 = (data_arr[n_of_data])(i ,j ,0,comp);
558 amrex::Real f2 = (data_arr[n_of_data])(i ,j+1,0,comp) - f1;
559 amrex::Real f3 = (data_arr[n_of_data])(i+1,j ,0,comp) - f1;
560 amrex::Real f4 = (data_arr[n_of_data])(i+1,j+1,0,comp) - f1;
561
562 // Node locations
563 amrex::Real x1 = loc_arr(i ,j ,0,0);
564 amrex::Real x2 = loc_arr(i ,j+1,0,0) - x1;
565 amrex::Real x3 = loc_arr(i+1,j ,0,0) - x1;
566 amrex::Real x4 = loc_arr(i+1,j+1,0,0) - x1;
567
568 amrex::Real y1 = loc_arr(i ,j ,0,1);
569 amrex::Real y2 = loc_arr(i ,j+1,0,1) - y1;
570 amrex::Real y3 = loc_arr(i+1,j ,0,1) - y1;
571 amrex::Real y4 = loc_arr(i+1,j+1,0,1) - y1;
572
573 // Particle location
574 amrex::Real x = p.pos(0) - x1;
575 amrex::Real y = p.pos(1) - y1;
576
577 amrex::Real det = x2*x4*y3*(y2 - y4) - x3*x4*y2*(y3 - y4) - x2*x3*(y2 - y3)*y4;
578
579 amrex::Real b = ( f4*(x2*y2*y3 - x3*y2*y3) + (f2*(x3 - x4)*y3)*y4 + f3*(-(x2*y2*y4) + x4*y2*y4) ) / det;
580 amrex::Real c = ( -f2*x3*x4*y3 + f4*x2*x3*(-y2 + y3) + f3*x2*x4*(y2 - y4) + f2*x3*x4*y4 ) / det;
581 amrex::Real d = ( f2*x4*y3 + f4*(x3*y2 - x2*y3) - f2*x3*y4 + f3*(-(x4*y2) + x2*y4) ) / det;
582
583 amrex::Real f = b*x + c*y + d*x*y;
584 val[comp-start_comp] = f1 + f;
585
586#elif (AMREX_SPACEDIM == 3)
587
588 int k = p.idata(2);
589
590 //
591 // 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)
592 // All locations are calculated as offsets from (lo,lo,lo), which is only need 7 values are needed
593 //
594 Array1D<Real,0,7> x_vals;
595 Array1D<Real,0,7> y_vals;
596 Array1D<Real,0,7> z_vals;
597 Array1D<Real,1,7> b;
598
599 amrex::Real x1, y1, z1, f1;
600
601 x1 = loc_arr(i ,j ,k ,0);
602 x_vals(1) = loc_arr(i+1,j ,k ,0) - x1;
603 x_vals(2) = loc_arr(i ,j+1,k ,0) - x1;
604 x_vals(3) = loc_arr(i+1,j+1,k ,0) - x1;
605 x_vals(4) = loc_arr(i ,j ,k+1,0) - x1;
606 x_vals(5) = loc_arr(i+1,j ,k+1,0) - x1;
607 x_vals(6) = loc_arr(i ,j+1,k+1,0) - x1;
608 x_vals(7) = loc_arr(i+1,j+1,k+1,0) - x1;
609
610 y1 = loc_arr(i ,j ,k ,1);
611 y_vals(1) = loc_arr(i+1,j ,k ,1) - y1;
612 y_vals(2) = loc_arr(i ,j+1,k ,1) - y1;
613 y_vals(3) = loc_arr(i+1,j+1,k ,1) - y1;
614 y_vals(4) = loc_arr(i ,j ,k+1,1) - y1;
615 y_vals(5) = loc_arr(i+1,j ,k+1,1) - y1;
616 y_vals(6) = loc_arr(i ,j+1,k+1,1) - y1;
617 y_vals(7) = loc_arr(i+1,j+1,k+1,1) - y1;
618
619 z1 = loc_arr(i ,j ,k ,2);
620 z_vals(1) = loc_arr(i+1,j ,k ,2) - z1;
621 z_vals(2) = loc_arr(i ,j+1,k ,2) - z1;
622 z_vals(3) = loc_arr(i+1,j+1,k ,2) - z1;
623 z_vals(4) = loc_arr(i ,j ,k+1,2) - z1;
624 z_vals(5) = loc_arr(i+1,j ,k+1,2) - z1;
625 z_vals(6) = loc_arr(i ,j+1,k+1,2) - z1;
626 z_vals(7) = loc_arr(i+1,j+1,k+1,2) - z1;
627
628 f1 = (data_arr[n_of_data])(i ,j ,k ,comp);
629 b(1) = (data_arr[n_of_data])(i+1,j ,k ,comp) - f1;
630 b(2) = (data_arr[n_of_data])(i ,j+1,k ,comp) - f1;
631 b(3) = (data_arr[n_of_data])(i+1,j+1,k ,comp) - f1;
632 b(4) = (data_arr[n_of_data])(i ,j ,k+1,comp) - f1;
633 b(5) = (data_arr[n_of_data])(i+1,j ,k+1,comp) - f1;
634 b(6) = (data_arr[n_of_data])(i ,j+1,k+1,comp) - f1;
635 b(7) = (data_arr[n_of_data])(i+1,j+1,k+1,comp) - f1;
636
637 Array2D<Real,1,7,1,7> a;
638 Array1D<int,1,7> ip;
639
640 int neq = 7;
641
642 for (int n=1; n<neq+1; n++) {
643 a(n,1) = x_vals(n);
644 a(n,2) = y_vals(n);
645 a(n,3) = z_vals(n);
646 a(n,4) = x_vals(n)*y_vals(n);
647 a(n,5) = x_vals(n)*z_vals(n);
648 a(n,6) = z_vals(n)*y_vals(n);
649 a(n,7) = x_vals(n)*y_vals(n)*z_vals(n);
650 }
651
652 particle_interp_decomp(a,ip,neq);
653 particle_interp_solve(a,b,ip,neq);
654
655 amrex::Real px, py, pz;
656 px = p.pos(0) - x1;
657 py = p.pos(1) - y1;
658 pz = p.pos(2) - z1;
659
660 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;
661 val[comp-start_comp] = static_cast<ParticleReal>(f + f1);
662#endif
663
664 } // comp
665#endif
666}
667
668} // namespace amrex
669#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:129
cic_interpolate_mapped_z
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)
Linearly interpolates the mesh data to the particle position from cell-centered data on a terrain-fit...
Definition AMReX_TracerParticle_mod_K.H:167
linear_interpolate_to_particle_z
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)
Linearly interpolates the mesh data to the particle position from mesh data. This general form can ha...
Definition AMReX_TracerParticle_mod_K.H:252
cic_interpolate_nd_mapped_z
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)
Linearly interpolates the mesh data to the particle position from node-centered data....
Definition AMReX_TracerParticle_mod_K.H:204
mac_interpolate_mapped_z
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)
Linearly interpolates the mesh data to the particle position from face-centered data on a terrain-fit...
Definition AMReX_TracerParticle_mod_K.H:226
cic_interpolate_nd_mapped
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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=AMREX_SPACEDIM)
Linearly interpolates the mesh data to the particle position from node-centered data on a general map...
Definition AMReX_TracerParticle_mod_K.H:506
particle_interp_decomp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:422
particle_interp_solve
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:464
cic_interpolate_cc_mapped_z
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)
Linearly interpolates the mesh data to the particle position from cell-centered data on a terrain-fit...
Definition AMReX_TracerParticle_mod_K.H:183
linear_interpolate_to_particle_mapped
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:526
AMReX_TracerParticles.H
amrex::IntVectND< AMREX_SPACEDIM >::TheUnitVector
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE 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:680
amrex::IntVectND< AMREX_SPACEDIM >::TheZeroVector
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE 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:670
amrex
Definition AMReX_Amr.cpp:49
amrex::cic_interpolate_nd
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cic_interpolate_nd(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, amrex::ParticleReal *val, int M=AMREX_SPACEDIM)
Linearly interpolates the mesh data to the particle position from node-centered data.
Definition AMReX_TracerParticle_mod_K.H:56
amrex::linear_interpolate_to_particle
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void linear_interpolate_to_particle(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, 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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cic_interpolate(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, amrex::ParticleReal *val, int M=AMREX_SPACEDIM)
Linearly interpolates the mesh data to the particle position from cell-centered data.
Definition AMReX_TracerParticle_mod_K.H:27
amrex::mac_interpolate
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mac_interpolate(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, amrex::ParticleReal *val)
Linearly interpolates the mesh data to the particle position from face-centered data....
Definition AMReX_TracerParticle_mod_K.H:75
amrex::IntVect
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:127
amrex::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:230
amrex::cic_interpolate_cc
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cic_interpolate_cc(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, amrex::ParticleReal *val, int M=AMREX_SPACEDIM)
Definition AMReX_TracerParticle_mod_K.H:38
amrex::Array4
Definition AMReX_Array4.H:61
amrex::GpuArray
Definition AMReX_Array.H:34
amrex::GpuArray::data
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * data() const noexcept
Definition AMReX_Array.H:56