Block-Structured AMR Software Framework
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 
14 namespace amrex {
15 
16 //
17 // **********************************************************************
18 // Regular coordinates
19 // **********************************************************************
20 //
21 
25 template <typename P>
26 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
27 void 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 
36 template <typename P>
37 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
38 void 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 
54 template <typename P>
55 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
56 void 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 
73 template <typename P>
74 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
75 void 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 
98 template <typename P>
99 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
100 void 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 
165 template <typename P>
166 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
167 void 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 
181 template <typename P>
182 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
183 void 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 
202 template <typename P>
203 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
204 void 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 
224 template <typename P>
225 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
226 void 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 
250 template <typename P>
251 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
252 void 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 
421 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
422 void 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 
463 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
464 void 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 
504 template <typename P>
505 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
506 void 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 
524 template <typename P>
525 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
526 void 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
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
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 constexpr AMREX_FORCE_INLINE 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:682
amrex::IntVectND< AMREX_SPACEDIM >::TheZeroVector
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE 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:672
abs
#define abs(x)
Definition: complex-type.h:85
amrex::SundialsUserFun::f
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
amrex::openbc::M
static constexpr int M
Definition: AMReX_OpenBC.H:13
amrex::openbc::P
static constexpr int P
Definition: AMReX_OpenBC.H:14
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:111
amrex::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition: AMReX.cpp:221
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:33
amrex::GpuArray::data
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * data() const noexcept
Definition: AMReX_Array.H:55