From 88132596437a4bdd7581bc3bf3d47ce53c1f3b6e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 11 Sep 2024 12:17:00 -0700 Subject: [PATCH] fix to fast time stepping and to terrain diffusion of QKE (#1798) * fix to fast time stepping and to terrain diffusion of QKE * uncomment qfx2_z * dt should be dtau * fix github actions --- .github/workflows/c-linter.yml | 2 +- Source/Diffusion/DiffusionSrcForState_T.cpp | 50 ++++++++++++++++++- Source/Initialization/ERF_init1d.cpp | 20 ++++++++ .../ERF_init_from_input_sounding.cpp | 21 ++++++++ Source/TimeIntegration/ERF_fast_rhs_N.cpp | 17 +++---- Source/TimeIntegration/ERF_fast_rhs_T.cpp | 48 ++++++++++++------ 6 files changed, 130 insertions(+), 28 deletions(-) diff --git a/.github/workflows/c-linter.yml b/.github/workflows/c-linter.yml index b178182f2..2a7d447b7 100644 --- a/.github/workflows/c-linter.yml +++ b/.github/workflows/c-linter.yml @@ -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 diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 9bb0e1eec..1c0eae48a 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -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& mu_turb, const SolverChoice &solverChoice, @@ -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) { @@ -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) { @@ -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 { @@ -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); + } }); } diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 71d5ecbdf..45a629c44 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -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 fmf = {&base_state[lev ], &base_state[lev ]}; + Vector cmf = {&base_state[lev-1], &base_state[lev-1]}; + Vector ftime = {time, time}; + Vector 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 diff --git a/Source/Initialization/ERF_init_from_input_sounding.cpp b/Source/Initialization/ERF_init_from_input_sounding.cpp index 9c9901ad9..0daf2b496 100644 --- a/Source/Initialization/ERF_init_from_input_sounding.cpp +++ b/Source/Initialization/ERF_init_from_input_sounding.cpp @@ -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 fmf = {&base_state[lev ], &base_state[lev ]}; + Vector cmf = {&base_state[lev-1], &base_state[lev-1]}; + Vector ftime = {time, time}; + Vector 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); + } } /** diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index 03a09e928..36d9e2ba4 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -415,9 +415,6 @@ 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 @@ -425,12 +422,12 @@ void erf_fast_rhs_N (int step, int nrk, 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); @@ -453,8 +450,8 @@ 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); } } @@ -462,7 +459,7 @@ 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) { - 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) { diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 5f9b90328..7cf306284 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -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) { @@ -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); @@ -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); @@ -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); @@ -548,19 +562,21 @@ 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) { @@ -568,7 +584,7 @@ void erf_fast_rhs_T (int step, int nrk, } } } - 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) {