Skip to content

Commit 0dfb8f1

Browse files
asalmgrenasalmgren
committed
Try this version
1 parent 3c48ceb commit 0dfb8f1

File tree

1 file changed

+138
-134
lines changed

1 file changed

+138
-134
lines changed

‎Src/EB/AMReX_EB_StateRedistribute.cpp‎

Lines changed: 138 additions & 134 deletions
Original file line numberDiff line numberDiff line change
@@ -117,70 +117,76 @@ MLStateRedistribute ( Box const& bx, int ncomp,
117117

118118
for (int n = 0; n < ncomp; n++)
119119
{
120-
121-
// Define Qhat (from Berger and Guliani)
122-
// Here we initialize Qhat to equal U_in on all cells in bxg3 so that
123-
// in the event we need to use Qhat 3 cells out from the bx limits
124-
// in a modified slope computation, we have a value of Qhat to use.
125-
// But we only modify Qhat inside bxg2
126-
amrex::ParallelFor(bxg3,
127-
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
128-
{
129-
if (vfrac(i,j,k) > 0.0 && bxg2.contains(IntVect(AMREX_D_DECL(i,j,k)))
130-
&& domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k))))
120+
// Define Qhat (from Berger and Guliani)
121+
// Here we initialize Qhat to equal U_in on all cells in bxg3 so that
122+
// in the event we need to use Qhat 3 cells out from the bx limits
123+
// in a modified slope computation, we have a value of Qhat to use.
124+
// But we only modify Qhat inside bxg2
125+
amrex::ParallelFor(bxg3,
126+
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
131127
{
132-
if (itracker(i,j,k,0) == 0)
133-
{
134-
qt(i,j,k,0 ) = alpha(i,j,k,0) * U_in(i,j,k,n) * vfrac(i,j,k) / nbhd_vol(i,j,k);
135-
Qhat(i,j,k,n) = qt(i,j,k,0);
136-
}
137-
else
128+
if (vfrac(i,j,k) > 0.0 && bxg2.contains(IntVect(AMREX_D_DECL(i,j,k)))
129+
&& domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k))))
138130
{
139-
Qhat(i,j,k,n) = 0.;
140-
141-
// This loops over (i,j,k) and the neighbors of (i,j,k)
142-
for (int i_nbor = 0; i_nbor <= itracker(i,j,k,0); i_nbor++)
131+
if (itracker(i,j,k,0) == 0)
143132
{
144-
int r = i; int s = j; int t = k;
145-
Real fac = alpha(i,j,k,0);
146-
if (i_nbor > 0) {
147-
r = i+imap[itracker(i,j,k,i_nbor)];
148-
s = j+jmap[itracker(i,j,k,i_nbor)];
149-
t = k+kmap[itracker(i,j,k,i_nbor)];
150-
fac = alpha(i,j,k,1) / nrs(r,s,t);
151-
}
133+
qt(i,j,k,0 ) = alpha(i,j,k,0) * U_in(i,j,k,n) * vfrac(i,j,k) / nbhd_vol(i,j,k);
134+
Qhat(i,j,k,n) = qt(i,j,k,0);
135+
}
136+
else
137+
{
138+
Qhat(i,j,k,n) = 0.;
152139

153-
if (domain_per_grown.contains(IntVect(AMREX_D_DECL(r,s,t))))
140+
// This loops over (i,j,k) and the neighbors of (i,j,k)
141+
for (int i_nbor = 0; i_nbor <= itracker(i,j,k,0); i_nbor++)
154142
{
155-
qt(i,j,k,i_nbor) = fac * U_in(r,s,t,n) * vfrac(r,s,t) / nbhd_vol(i,j,k);
156-
Qhat(i,j,k,n) += qt(i,j,k,i_nbor);
143+
int r = i; int s = j; int t = k;
144+
Real fac = alpha(i,j,k,0);
145+
if (i_nbor > 0) {
146+
r = i+imap[itracker(i,j,k,i_nbor)];
147+
s = j+jmap[itracker(i,j,k,i_nbor)];
148+
t = k+kmap[itracker(i,j,k,i_nbor)];
149+
fac = alpha(i,j,k,1) / nrs(r,s,t);
150+
}
151+
152+
if (domain_per_grown.contains(IntVect(AMREX_D_DECL(r,s,t))))
153+
{
154+
qt(i,j,k,i_nbor) = fac * U_in(r,s,t,n) * vfrac(r,s,t) / nbhd_vol(i,j,k);
155+
Qhat(i,j,k,n) += qt(i,j,k,i_nbor);
156+
}
157157
}
158158
}
159+
} else {
160+
Qhat(i,j,k,n) = U_in(i,j,k,n);
159161
}
160-
} else {
161-
Qhat(i,j,k,n) = U_in(i,j,k,n);
162-
}
163-
});
162+
});
164163

165-
//
166-
// ****************************************************************************************
167-
//
164+
//
165+
// ****************************************************************************************
166+
//
168167

169-
amrex::ParallelFor(bxg1,
170-
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
171-
{
172-
if (vfrac(i,j,k) > 0.0)
168+
amrex::ParallelFor(bxg1,
169+
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
173170
{
174-
if (itracker(i,j,k,0) == 0)
171+
if (vfrac(i,j,k) > 0.0)
175172
{
176-
if (domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k))))
173+
// NO NEIGBHORS EXCEPT SELF
174+
if (itracker(i,j,k,0) == 0)
177175
{
178-
U_out(i,j,k,n) = alpha(i,j,k,0)*qt(i,j,k,0);
179-
}
180-
}
181-
else // Multiple neighbors
182-
{
183-
// This loops over (i,j,k) and the neighbors of (i,j,k)
176+
if (domain_per_grown.contains(IntVect(AMREX_D_DECL(i,j,k))) )
177+
{
178+
Real fac = alpha(i,j,k,0) * nrs(i,j,k);
179+
180+
//
181+
// Update U_out(i,j,k) with qt(i,j,k,0)
182+
//
183+
U_out(i,j,k,n) += fac * qt(i,j,k,0) / nrs(i,j,k);
184+
185+
} // contains
186+
187+
} else {
188+
189+
// MULTIPLE NEIGBHORS
184190
for (int i_nbor = 0; i_nbor <= itracker(i,j,k,0); i_nbor++)
185191
{
186192
int r = i; int s = j; int t = k;
@@ -213,76 +219,76 @@ MLStateRedistribute ( Box const& bx, int ncomp,
213219
for(int iii(-1); iii<=1; iii++) {
214220
if (flag(i,j,k).isConnected(iii,jjj,kkk))
215221
{
216-
int rr = i+iii; int ss = j+jjj; int tt = k+kkk;
222+
int rr = i+iii; int ss = j+jjj; int tt = k+kkk;
217223

218-
x_max = amrex::max(x_max, cent_hat(rr,ss,tt,0)+static_cast<Real>(iii));
219-
x_min = amrex::min(x_min, cent_hat(rr,ss,tt,0)+static_cast<Real>(iii));
220-
y_max = amrex::max(y_max, cent_hat(rr,ss,tt,1)+static_cast<Real>(jjj));
221-
y_min = amrex::min(y_min, cent_hat(rr,ss,tt,1)+static_cast<Real>(jjj));
224+
x_max = amrex::max(x_max, cent_hat(rr,ss,tt,0)+static_cast<Real>(iii));
225+
x_min = amrex::min(x_min, cent_hat(rr,ss,tt,0)+static_cast<Real>(iii));
226+
y_max = amrex::max(y_max, cent_hat(rr,ss,tt,1)+static_cast<Real>(jjj));
227+
y_min = amrex::min(y_min, cent_hat(rr,ss,tt,1)+static_cast<Real>(jjj));
222228
#if (AMREX_SPACEDIM == 3)
223-
z_max = amrex::max(z_max, cent_hat(rr,ss,tt,2)+static_cast<Real>(kkk));
224-
z_min = amrex::min(z_min, cent_hat(rr,ss,tt,2)+static_cast<Real>(kkk));
229+
z_max = amrex::max(z_max, cent_hat(rr,ss,tt,2)+static_cast<Real>(kkk));
230+
z_min = amrex::min(z_min, cent_hat(rr,ss,tt,2)+static_cast<Real>(kkk));
225231
#endif
226232
}
227-
AMREX_D_TERM(},},})
233+
AMREX_D_TERM(},},})
228234

229-
// If we need to grow the stencil, we let it be -nx:nx in the x-direction,
230-
// for example. Note that nx,ny,nz are either 1 or 2
231-
if ( (x_max-x_min) < slope_stencil_min_width ) { nx = 2; }
232-
if ( (y_max-y_min) < slope_stencil_min_width ) { ny = 2; }
235+
// If we need to grow the stencil, we let it be -nx:nx in the x-direction,
236+
// for example. Note that nx,ny,nz are either 1 or 2
237+
if ( (x_max-x_min) < slope_stencil_min_width ) { nx = 2; }
238+
if ( (y_max-y_min) < slope_stencil_min_width ) { ny = 2; }
233239
#if (AMREX_SPACEDIM == 3)
234-
if ( (z_max-z_min) < slope_stencil_min_width ) { nz = 2; }
240+
if ( (z_max-z_min) < slope_stencil_min_width ) { nz = 2; }
235241
#endif
236-
bool extdir_ilo = (d_bcrec_ptr[n].lo(0) == amrex::BCType::ext_dir ||
237-
d_bcrec_ptr[n].lo(0) == amrex::BCType::hoextrap);
238-
bool extdir_ihi = (d_bcrec_ptr[n].hi(0) == amrex::BCType::ext_dir ||
239-
d_bcrec_ptr[n].hi(0) == amrex::BCType::hoextrap);
240-
bool extdir_jlo = (d_bcrec_ptr[n].lo(1) == amrex::BCType::ext_dir ||
241-
d_bcrec_ptr[n].lo(1) == amrex::BCType::hoextrap);
242-
bool extdir_jhi = (d_bcrec_ptr[n].hi(1) == amrex::BCType::ext_dir ||
243-
d_bcrec_ptr[n].hi(1) == amrex::BCType::hoextrap);
242+
bool extdir_ilo = (d_bcrec_ptr[n].lo(0) == amrex::BCType::ext_dir ||
243+
d_bcrec_ptr[n].lo(0) == amrex::BCType::hoextrap);
244+
bool extdir_ihi = (d_bcrec_ptr[n].hi(0) == amrex::BCType::ext_dir ||
245+
d_bcrec_ptr[n].hi(0) == amrex::BCType::hoextrap);
246+
bool extdir_jlo = (d_bcrec_ptr[n].lo(1) == amrex::BCType::ext_dir ||
247+
d_bcrec_ptr[n].lo(1) == amrex::BCType::hoextrap);
248+
bool extdir_jhi = (d_bcrec_ptr[n].hi(1) == amrex::BCType::ext_dir ||
249+
d_bcrec_ptr[n].hi(1) == amrex::BCType::hoextrap);
244250
#if (AMREX_SPACEDIM == 3)
245-
bool extdir_klo = (d_bcrec_ptr[n].lo(2) == amrex::BCType::ext_dir ||
246-
d_bcrec_ptr[n].lo(2) == amrex::BCType::hoextrap);
247-
bool extdir_khi = (d_bcrec_ptr[n].hi(2) == amrex::BCType::ext_dir ||
248-
d_bcrec_ptr[n].hi(2) == amrex::BCType::hoextrap);
251+
bool extdir_klo = (d_bcrec_ptr[n].lo(2) == amrex::BCType::ext_dir ||
252+
d_bcrec_ptr[n].lo(2) == amrex::BCType::hoextrap);
253+
bool extdir_khi = (d_bcrec_ptr[n].hi(2) == amrex::BCType::ext_dir ||
254+
d_bcrec_ptr[n].hi(2) == amrex::BCType::hoextrap);
249255
#endif
250256

251-
// Compute slopes of Qhat (which is the sum of the qt's) then use
252-
// that for each qt separately
253-
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> slopes_eb;
254-
if (nx*ny*nz == 1) {
255-
// Compute slope using 3x3x3 stencil
256-
slopes_eb = amrex_calc_slopes_extdir_eb(
257-
i,j,k,n,Qhat,cent_hat,vfrac,
258-
AMREX_D_DECL(fcx,fcy,fcz),flag,
259-
AMREX_D_DECL(extdir_ilo, extdir_jlo, extdir_klo),
260-
AMREX_D_DECL(extdir_ihi, extdir_jhi, extdir_khi),
261-
AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
262-
AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
263-
max_order);
264-
} else {
265-
// Compute slope using grown stencil (no larger than 5x5x5)
266-
slopes_eb = amrex_calc_slopes_extdir_eb_grown(
267-
i,j,k,n,AMREX_D_DECL(nx,ny,nz),
268-
Qhat,cent_hat,vfrac,
269-
AMREX_D_DECL(fcx,fcy,fcz),flag,
270-
AMREX_D_DECL(extdir_ilo, extdir_jlo, extdir_klo),
271-
AMREX_D_DECL(extdir_ihi, extdir_jhi, extdir_khi),
272-
AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
273-
AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
274-
max_order);
275-
}
276-
277-
// We do the limiting separately because this limiter limits the slope based on the values
278-
// extrapolated to the cell centroid (cent_hat) locations - unlike the limiter in amrex
279-
// which bases the limiting on values extrapolated to the face centroids.
280-
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> lim_slope =
281-
amrex_calc_centroid_limiter(i,j,k,n,Qhat,flag,slopes_eb,cent_hat);
282-
283-
AMREX_D_TERM(lim_slope[0] *= slopes_eb[0];,
284-
lim_slope[1] *= slopes_eb[1];,
285-
lim_slope[2] *= slopes_eb[2];);
257+
// Compute slopes of Qhat (which is the sum of the qt's) then use
258+
// that for each qt separately
259+
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> slopes_eb;
260+
if (nx*ny*nz == 1) {
261+
// Compute slope using 3x3x3 stencil
262+
slopes_eb = amrex_calc_slopes_extdir_eb(
263+
i,j,k,n,Qhat,cent_hat,vfrac,
264+
AMREX_D_DECL(fcx,fcy,fcz),flag,
265+
AMREX_D_DECL(extdir_ilo, extdir_jlo, extdir_klo),
266+
AMREX_D_DECL(extdir_ihi, extdir_jhi, extdir_khi),
267+
AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
268+
AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
269+
max_order);
270+
} else {
271+
// Compute slope using grown stencil (no larger than 5x5x5)
272+
slopes_eb = amrex_calc_slopes_extdir_eb_grown(
273+
i,j,k,n,AMREX_D_DECL(nx,ny,nz),
274+
Qhat,cent_hat,vfrac,
275+
AMREX_D_DECL(fcx,fcy,fcz),flag,
276+
AMREX_D_DECL(extdir_ilo, extdir_jlo, extdir_klo),
277+
AMREX_D_DECL(extdir_ihi, extdir_jhi, extdir_khi),
278+
AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
279+
AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
280+
max_order);
281+
}
282+
283+
// We do the limiting separately because this limiter limits the slope based on the values
284+
// extrapolated to the cell centroid (cent_hat) locations - unlike the limiter in amrex
285+
// which bases the limiting on values extrapolated to the face centroids.
286+
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> lim_slope =
287+
amrex_calc_centroid_limiter(i,j,k,n,Qhat,flag,slopes_eb,cent_hat);
288+
289+
AMREX_D_TERM(lim_slope[0] *= slopes_eb[0];,
290+
lim_slope[1] *= slopes_eb[1];,
291+
lim_slope[2] *= slopes_eb[2];);
286292

287293
for (int r_nbor = 0; r_nbor <= itracker(i,j,k,0); r_nbor++)
288294
{
@@ -320,23 +326,23 @@ MLStateRedistribute ( Box const& bx, int ncomp,
320326

321327
if (as_crse) {
322328

323-
// Covered (by fine) to uncovered (by fine)
324-
if ( bx.contains(IntVect(AMREX_D_DECL(r,s,t))) &&
325-
flag_as_crse( r, s, t) == amrex_yafluxreg_crse_fine_boundary_cell &&
326-
flag_as_crse(ii,jj,kk) == amrex_yafluxreg_fine_cell )
327-
{
328-
drho_as_crse(r,s,t,n) -= fac*update/nrs(r,s,t) * fac_for_deltaR;
329-
}
330-
331-
// Uncovered (by fine) to covered (by fine)
332-
if ( bx.contains(IntVect(AMREX_D_DECL(ii,jj,kk))) &&
333-
flag_as_crse( r, s, t) == amrex_yafluxreg_fine_cell &&
334-
flag_as_crse(ii,jj,kk) == amrex_yafluxreg_crse_fine_boundary_cell )
335-
{
336-
drho_as_crse(ii,jj,kk,n) += fac * update / nrs(r,s,t) *
337-
(vfrac(r,s,t) / vfrac(ii,jj,kk)) * fac_for_deltaR;
338-
}
339-
} // as_crse
329+
// Covered (by fine) to uncovered (by fine)
330+
if ( bx.contains(IntVect(AMREX_D_DECL(r,s,t))) &&
331+
flag_as_crse( r, s, t) == amrex_yafluxreg_crse_fine_boundary_cell &&
332+
flag_as_crse(ii,jj,kk) == amrex_yafluxreg_fine_cell )
333+
{
334+
drho_as_crse(r,s,t,n) -= fac*update/nrs(r,s,t) * fac_for_deltaR;
335+
}
336+
337+
// Uncovered (by fine) to covered (by fine)
338+
if ( bx.contains(IntVect(AMREX_D_DECL(ii,jj,kk))) &&
339+
flag_as_crse( r, s, t) == amrex_yafluxreg_fine_cell &&
340+
flag_as_crse(ii,jj,kk) == amrex_yafluxreg_crse_fine_boundary_cell )
341+
{
342+
drho_as_crse(ii,jj,kk,n) += fac * update / nrs(r,s,t) *
343+
(vfrac(r,s,t) / vfrac(ii,jj,kk)) * fac_for_deltaR;
344+
}
345+
} // as_crse
340346

341347
if (as_fine) {
342348

@@ -352,16 +358,14 @@ MLStateRedistribute ( Box const& bx, int ncomp,
352358
dm_as_fine(r,s,t,n) += fac*update/nrs(r,s,t) * vfrac(r,s,t) * fac_for_deltaR;
353359
}
354360
} // as_fine
361+
355362
} // r_nbor
356363
} // bx contains
357364
} // i_nbor
358-
} // more than one nbor
359-
} // vfrac
360-
});
365+
} // ELSE
366+
} // vfrac
361367

362-
//
363-
// ****************************************************************************************
364-
//
368+
});
365369

366370
} // loop over n = 0, ncomp-1
367371

0 commit comments

Comments
 (0)