Skip to content

Commit

Permalink
fix to fast time stepping and to terrain diffusion of QKE (#1798)
Browse files Browse the repository at this point in the history 
* fix to fast time stepping and to terrain diffusion of QKE

* uncomment qfx2_z

* dt should be dtau

* fix github actions
  • Loading branch information
@asalmgren
asalmgren authored Sep 11, 2024
1 parent e7e76f3 commit 8813259
Show file tree
Hide file tree
Showing 6 changed files with 130 additions and 28 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/c-linter.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ jobs:
run: python3 Submodules/cpp-linter-action/run_on_changed_files.py ${{ github.event.pull_request.base.sha }} ${{ github.event.pull_request.head.sha }} -header-filter=Source/ -ignore-files="AMReX|GoogleTest" -run-linter

- name: Archive clang tidy report
uses: actions/upload-artifact@v1
uses: actions/upload-artifact@v4
with:
name: clang-tidy-report
path: clang-tidy-report.txt
50 changes: 49 additions & 1 deletion Source/Diffusion/DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
Array4< Real>& qfx1_x,
Array4< Real>& qfx1_y,
Array4< Real>& qfx1_z,
Array4< Real>& /*qfx2_z*/,
Array4< Real>& qfx2_z,
Array4< Real>& diss,
const Array4<const Real>& mu_turb,
const SolverChoice &solverChoice,
Expand Down Expand Up @@ -319,6 +319,18 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
} else {
zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
}
});
// Constant rho*alpha & Turb model
} else if (l_turb) {
Expand Down Expand Up @@ -427,6 +439,18 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
} else {
zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
}
});
// Constant alpha & no LES/PBL model
} else if(l_consA) {
Expand Down Expand Up @@ -531,6 +555,18 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
} else {
zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
}
});
// Constant rho*alpha & no LES/PBL model
} else {
Expand Down Expand Up @@ -635,6 +671,18 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
} else {
zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = zflux(i,j,k,qty_index);
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = zflux(i,j,k,qty_index);
}
});
}

Expand Down
20 changes: 20 additions & 0 deletions Source/Initialization/ERF_init1d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,6 +278,26 @@ ERF::initHSE (int lev)
base_state[lev].ParallelCopy(new_base_state,0,0,3,1,1);
}

// Make sure to fill the ghost cells of the base state
if (lev > 0) {
// Interp all three components: rho, p, pi
int icomp = 0; int bccomp = 0; int ncomp = 3;

PhysBCFunctNoOp null_bc;
Interpolater* mapper = &cell_cons_interp;

Real time = 0.;

Vector<MultiFab*> fmf = {&base_state[lev ], &base_state[lev ]};
Vector<MultiFab*> cmf = {&base_state[lev-1], &base_state[lev-1]};
Vector<Real> ftime = {time, time};
Vector<Real> ctime = {time, time};
FillPatchTwoLevels(base_state[lev], time,
cmf, ctime, fmf, ftime,
icomp, icomp, ncomp, geom[lev-1], geom[lev],
null_bc, 0, null_bc, 0, refRatio(lev-1),
mapper, domain_bcs_type, bccomp);
}
}

void
Expand Down
21 changes: 21 additions & 0 deletions Source/Initialization/ERF_init_from_input_sounding.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,27 @@ ERF::init_from_input_sounding (int lev)
input_sounding_data);

} //mfi

// Make sure to fill the ghost cells of the base state
if (lev > 0) {
// Interp all three components: rho, p, pi
int icomp = 0; int bccomp = 0; int ncomp = 3;

PhysBCFunctNoOp null_bc;
Interpolater* mapper = &cell_cons_interp;

Real time = 0.;

Vector<MultiFab*> fmf = {&base_state[lev ], &base_state[lev ]};
Vector<MultiFab*> cmf = {&base_state[lev-1], &base_state[lev-1]};
Vector<Real> ftime = {time, time};
Vector<Real> ctime = {time, time};
FillPatchTwoLevels(base_state[lev], time,
cmf, ctime, fmf, ftime,
icomp, icomp, ncomp, geom[lev-1], geom[lev],
null_bc, 0, null_bc, 0, refRatio(lev-1),
mapper, domain_bcs_type, bccomp);
}
}

/**
Expand Down
17 changes: 7 additions & 10 deletions Source/TimeIntegration/ERF_fast_rhs_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -415,22 +415,19 @@ void erf_fast_rhs_N (int step, int nrk,
Box b2d = tbz; // Copy constructor
b2d.setRange(2,0);

// dt is the timestep for the RK stage, so dtau = facinv * dt
Real dt = dtau / facinv;

{
BL_PROFILE("fast_rhs_b2d_loop");
#ifdef AMREX_USE_GPU
auto const lo = lbound(bx);
auto const hi = ubound(bx);
ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
{
// w at bottom boundary of grid is 0 if at domain boundary, otherwise w = w_old + dt * slow_rhs
RHS_a(i,j,lo.z) = dt * slow_rhs_rho_w(i,j,lo.z);
// w at bottom boundary of grid is 0 if at domain boundary, otherwise w = w_old + dtau * slow_rhs
RHS_a(i,j,lo.z) = dtau * slow_rhs_rho_w(i,j,lo.z);

// w at top boundary of grid is 0 if at domain boundary, otherwise w = w_old + dt * slow_rhs
// w at top boundary of grid is 0 if at domain boundary, otherwise w = w_old + dtau * slow_rhs
// TODO TODO: Note that if we ever change this, we will need to include it in avg_zmom at the top
RHS_a(i,j,hi.z+1) = dt * slow_rhs_rho_w(i,j,hi.z+1);
RHS_a(i,j,hi.z+1) = dtau * slow_rhs_rho_w(i,j,hi.z+1);

// w = specified Dirichlet value at k = lo.z
soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
Expand All @@ -453,16 +450,16 @@ void erf_fast_rhs_N (int step, int nrk,
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
// w at bottom boundary of grid is 0 if at domain boundary, otherwise w_old + dt * slow_rhs
RHS_a (i,j,lo.z) = stage_zmom(i,j,lo.z) + dt * slow_rhs_rho_w(i,j,lo.z);
// w at bottom boundary of grid is 0 if at domain boundary, otherwise w_old + dtau * slow_rhs
RHS_a (i,j,lo.z) = stage_zmom(i,j,lo.z) + dtau * slow_rhs_rho_w(i,j,lo.z);
soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
}
}
// Note that if we ever change this, we will need to include it in avg_zmom at the top
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
RHS_a (i,j,hi.z+1) = dt * slow_rhs_rho_w(i,j,hi.z+1);
RHS_a (i,j,hi.z+1) = dtau * slow_rhs_rho_w(i,j,hi.z+1);
}
}
for (int k = lo.z+1; k <= hi.z+1; ++k) {
Expand Down
48 changes: 32 additions & 16 deletions Source/TimeIntegration/ERF_fast_rhs_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,6 +231,12 @@ void erf_fast_rhs_T (int step, int nrk,
// *********************************************************************
{
BL_PROFILE("fast_rhs_xymom_T");

const auto& tbx_lo = lbound(tbx);
const auto& tbx_hi = ubound(tbx);
const auto& tby_lo = lbound(tby);
const auto& tby_hi = ubound(tby);

ParallelFor(tbx, tby,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Expand All @@ -257,6 +263,11 @@ void erf_fast_rhs_T (int step, int nrk,

new_drho_u(i, j, k) = old_drho_u(i,j,k) + dtau * fast_rhs_rho_u
+ dtau * slow_rhs_rho_u(i,j,k);
if (k == tbx_lo.z && k != 0) {
new_drho_u(i,j,k-1) = new_drho_u(i,j,k);
} else if (k == tbx_hi.z) {
new_drho_u(i,j,k+1) = new_drho_u(i,j,k);
}

avg_xmom(i,j,k) += facinv*new_drho_u(i,j,k);

Expand Down Expand Up @@ -288,6 +299,12 @@ void erf_fast_rhs_T (int step, int nrk,
new_drho_v(i, j, k) = old_drho_v(i,j,k) + dtau * fast_rhs_rho_v
+ dtau * slow_rhs_rho_v(i,j,k);

if (k == tby_lo.z && k != 0) {
new_drho_v(i,j,k-1) = new_drho_v(i,j,k);
} else if (k == tby_hi.z) {
new_drho_v(i,j,k+1) = new_drho_v(i,j,k);
}

avg_ymom(i,j,k) += facinv*new_drho_v(i,j,k);

cur_ymom(i,j,k) = stage_ymom(i,j,k) + new_drho_v(i,j,k);
Expand Down Expand Up @@ -519,17 +536,14 @@ void erf_fast_rhs_T (int step, int nrk,
auto const domlo = lbound(domain);
auto const domhi = ubound(domain);

// dt is the timestep for the RK stage, so dtau = facinv * dt
Real dt = dtau / facinv;

{
BL_PROFILE("fast_rhs_b2d_loop_t");
#ifdef AMREX_USE_GPU
ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
{
// w_klo, w_khi given by specified Dirichlet values
RHS_a(i,j,lo.z ) = dt * slow_rhs_rho_w(i,j,lo.z);
RHS_a(i,j,hi.z+1) = dt * slow_rhs_rho_w(i,j,hi.z+1);
RHS_a(i,j,lo.z ) = dtau * slow_rhs_rho_w(i,j,lo.z);
RHS_a(i,j,hi.z+1) = dtau * slow_rhs_rho_w(i,j,hi.z+1);

// w = specified Dirichlet value at k = lo.z
soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
Expand All @@ -548,27 +562,29 @@ void erf_fast_rhs_T (int step, int nrk,
#else
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
RHS_a(i,j,lo.z) = dt * slow_rhs_rho_w(i,j,lo.z);
for (int i = lo.x; i <= hi.x; ++i)
{
RHS_a(i,j,lo.z) = dtau * slow_rhs_rho_w(i,j,lo.z);
soln_a(i,j,lo.z) = RHS_a(i,j,lo.z) * inv_coeffB_a(i,j,lo.z);
}
}
}

for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
RHS_a(i,j,hi.z+1) = dt * slow_rhs_rho_w(i,j,hi.z+1);
}
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i)
{
RHS_a(i,j,hi.z+1) = dtau * slow_rhs_rho_w(i,j,hi.z+1);
soln_a(i,j,hi.z+1) = RHS_a(i,j,hi.z+1) * inv_coeffB_a(i,j,hi.z+1);
}
}
for (int k = lo.z+1; k <= hi.z+1; ++k) {

for (int k = lo.z+1; k <= hi.z; ++k) {
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
soln_a(i,j,k) = (RHS_a(i,j,k)-coeffA_a(i,j,k)*soln_a(i,j,k-1)) * inv_coeffB_a(i,j,k);
}
}
}
for (int k = hi.z; k >= lo.z; --k) {
for (int k = hi.z; k > lo.z; --k) {
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
Expand Down

0 comments on commit 8813259

Please sign in to comment.