Skip to content

Commit 4bc33c7

Browse files
asalmgrenasalmgren
authored
Modify how we interpolate velocity from faces to particle position; f… (#4598)
…or vertical velocities we always want to use the cell of the particle location ## Summary The previous interpolation scheme tested on where the particle location was relative to the plane between the face velocities; for the vertical velocities we don't want to do this because we always want to use the velocities at the k location of the particle position (otherwise there is a discrepancy between the particle being identified as in cell k but being below the threshhold for using the vertical velocities associated with cell k) ## Additional background Note this only affects cases where the velocity components are on faces and where the mesh is "distorted", i.e. with terrain-following coordinates ## Checklist The proposed changes: - [X] fix a bug or incorrect behavior in AMReX - [ ] add new capabilities to AMReX - [ ] changes answers in the test suite to more than roundoff level - [ ] are likely to significantly affect the results of downstream AMReX users - [ ] include documentation in the code and/or rst files, if appropriate
1 parent 68a31f9 commit 4bc33c7

File tree

1 file changed

+78
-33
lines changed

1 file changed

+78
-33
lines changed

‎Src/Particle/AMReX_TracerParticle_mod_K.H‎

Lines changed: 78 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -278,25 +278,46 @@ void linear_interpolate_to_particle_z (const P& p,
278278
amrex::ignore_unused(ly);
279279
int const j = p.idata(0);
280280

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-
281+
int j0 = j;
293282

294283
amrex::Real const xint = lx - static_cast<Real>(i0);
295284
amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
296-
amrex::Real height_at_px = sx[0] * hlo_xlo + sx[1] * hlo_xhi;
297285

298-
int const j0 = (amrex::Real(p.pos(1)) >= height_at_px) ? j : j-1;
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 ) );
299307

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+
//
300321
int yctr = 0;
301322
amrex::Real ht[4];
302323
for (int ii=0; ii < 2; ++ii) {
@@ -324,28 +345,52 @@ void linear_interpolate_to_particle_z (const P& p,
324345
amrex::Real sx[] = { amrex::Real(1.) - xint, xint};
325346
amrex::Real sy[] = { amrex::Real(1.) - yint, yint};
326347

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;
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+
}
344380
}
345-
}
346381

347-
int const k0 = (amrex::Real(p.pos(2)) >= height_at_pxy ) ? k : k-1;
382+
if (amrex::Real(p.pos(2)) < height_at_pxy ) {k0 = k-1;}
383+
384+
} // if not vertical direction
348385

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+
//
349394
int zctr = 0;
350395
amrex::Real ht[8];
351396
for (int ii = 0; ii < 2; ++ii) {

0 commit comments

Comments
 (0)