From 194cd609a2f2b76c5fe012bcfc9ca75e8e3206d1 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 3 Nov 2023 14:29:06 -0700 Subject: [PATCH 1/2] Refluxing2 (#1284) * update particle functionality * modify for refluxing plus some other cleanup such as terrain_type and coupling_type becoming enum structs * clean up and add breaks in switch statement --- Docs/sphinx_doc/Inputs.rst | 2 +- Docs/sphinx_doc/RegressionTests.rst | 2 +- Docs/sphinx_doc/theory/WetEquations.rst | 2 +- Exec/DevTests/MovingTerrain/inputs | 4 +- Exec/RegTests/ScalarAdvDiff/prob.cpp | 5 +- Source/Advection/Advection.H | 11 +- Source/Advection/AdvectionSrcForMom.cpp | 2 +- Source/Advection/AdvectionSrcForMom_N.H | 66 ++++-- Source/Advection/AdvectionSrcForMom_T.H | 54 +++-- Source/Advection/AdvectionSrcForScalars.H | 101 ++++++++ Source/Advection/AdvectionSrcForState.cpp | 223 +++++++++++------- Source/Advection/AdvectionSrcForState_N.H | 124 ++-------- Source/Advection/AdvectionSrcForState_T.H | 11 +- Source/Advection/Make.package | 1 + .../BoundaryConditions_zvel.cpp | 4 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 4 +- Source/BoundaryConditions/ERF_PhysBCFunct.H | 6 +- Source/DataStructs/DataStruct.H | 46 +++- Source/ERF.H | 12 +- Source/ERF.cpp | 50 ++-- Source/ERF_make_new_level.cpp | 6 +- Source/TimeIntegration/ERF_Advance.cpp | 16 +- Source/TimeIntegration/ERF_TimeStep.cpp | 4 +- Source/TimeIntegration/ERF_advance_dycore.cpp | 2 +- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 39 +-- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 26 +- Source/TimeIntegration/TI_fast_rhs_fun.H | 4 +- Source/TimeIntegration/TI_headers.H | 2 +- Source/TimeIntegration/TI_no_substep_fun.H | 2 +- Source/TimeIntegration/TI_slow_rhs_fun.H | 10 +- Source/Utils/Interpolation_UPW.H | 215 +---------------- Source/Utils/Interpolation_WENO.H | 100 +------- Source/Utils/Interpolation_WENO_Z.H | 149 +----------- Source/Utils/TerrainMetrics.cpp | 4 - .../DensityCurrent_detJ2_MT.i | 4 +- .../MovingTerrain_nosub/MovingTerrain_nosub.i | 4 +- .../MovingTerrain_sub/MovingTerrain_sub.i | 4 +- 37 files changed, 498 insertions(+), 823 deletions(-) create mode 100644 Source/Advection/AdvectionSrcForScalars.H diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 672ef6a4c..faaf43adf 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -860,7 +860,7 @@ List of Parameters | **erf.use_terrain** | use terrain-fitted | true / false | false | | | coordinates? | | | +-----------------------------+--------------------+--------------------+------------+ -| **erf.terrain_type** | static or moving? | 0 / 1 | 0 | +| **erf.terrain_type** | static or moving? | Static / Moving | Static | +-----------------------------+--------------------+--------------------+------------+ | **erf.terrain_smoothing** | specify terrain | 0, | 0 | | | following | 1, | | diff --git a/Docs/sphinx_doc/RegressionTests.rst b/Docs/sphinx_doc/RegressionTests.rst index d5c858677..56266be2d 100644 --- a/Docs/sphinx_doc/RegressionTests.rst +++ b/Docs/sphinx_doc/RegressionTests.rst @@ -45,7 +45,7 @@ The following problems are currently tested in the CI: | DensityCurrent_detJ2_MT | 256 4 64 | Symmetry | Periodic | SlipWall | None | use_terrain = true | | | | Outflow | | SlipWall | | uses zlevels | | | | Outflow | | SlipWall | | detJ = 2 everywhere | -| | | | | | | terrain_type = 1 | +| | | | | | | terrain_type = Moving | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | EkmanSpiral | 4 4 400 | Periodic | Periodic | NoSlipWall | Geo | +Coriolis | | | | | | SlipWall | | +gravity | diff --git a/Docs/sphinx_doc/theory/WetEquations.rst b/Docs/sphinx_doc/theory/WetEquations.rst index fc6c29d77..38288e484 100644 --- a/Docs/sphinx_doc/theory/WetEquations.rst +++ b/Docs/sphinx_doc/theory/WetEquations.rst @@ -38,7 +38,7 @@ The governing equations for this model are \frac{\partial \rho_d}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u}) \frac{\partial (\rho_d \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{u}) - - \frac{1}{1 + q_v + q_c} \nabla p^\prime - \nabla \cdot \tau + \mathbf{F} + \delta_{i,3}\mathbf{B} + \frac{1}{1 + q_v + q_c} ( \nabla p^\prime + \delta_{i,3}\mathbf{B} ) - \nabla \cdot \tau + \mathbf{F} \frac{\partial (\rho_d \theta_d)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \theta_d) + \nabla \cdot ( \rho_d \alpha_{T}\ \nabla \theta_d) + \frac{\theta_d L_v}{T_d C_{pd}} f_{cond} diff --git a/Exec/DevTests/MovingTerrain/inputs b/Exec/DevTests/MovingTerrain/inputs index 6c1e42970..8b0783fb4 100644 --- a/Exec/DevTests/MovingTerrain/inputs +++ b/Exec/DevTests/MovingTerrain/inputs @@ -19,8 +19,8 @@ zlo.type = "SlipWall" zhi.type = "SlipWall" # TERRRAIN GRID TYPE -erf.use_terrain = 1 # enable terrain stencils -erf.terrain_type = 1 # moving terrain +erf.use_terrain = true # enable terrain stencils +erf.terrain_type = Moving # moving terrain erf.terrain_smoothing = 2 # Sullivan 2004 approach # FOR NO SUBSTEPPING diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp b/Exec/RegTests/ScalarAdvDiff/prob.cpp index 643980677..a7d7dd560 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp @@ -118,8 +118,9 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0) / r0)); } else if (parms.prob_type == 13) { const Real r0_z = parms.rad_0 * (prob_hi[2] - prob_lo[2]); - const Real r2d_xz = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z)); //ellipse for mapfac shear validation - state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0_z) / r0_z)); + //ellipse for mapfac shear validation + const Real r2d_xz_ell = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z)); + state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz_ell, r0_z) / r0_z)); } else if (parms.prob_type == 14) { state(i, j, k, RhoScalar_comp) = std::cos(PI*x); } else { diff --git a/Source/Advection/Advection.H b/Source/Advection/Advection.H index 9b3121a57..1007f1f04 100644 --- a/Source/Advection/Advection.H +++ b/Source/Advection/Advection.H @@ -28,10 +28,12 @@ void AdvectionSrcForRhoAndTheta (const amrex::Box& bx, const amrex::Box& valid_b const amrex::Array4& mf_u, const amrex::Array4& mf_v, const AdvType horiz_adv_type, const AdvType vert_adv_type, - const int use_terrain); + const bool use_terrain); /** Compute advection tendency for all scalars other than density and potential temperature */ -void AdvectionSrcForScalars (const amrex::Box& bx, +void AdvectionSrcForScalars (int level, int finest_level, + const amrex::MFIter& mfi, + const amrex::Box& bx, const int icomp, const int ncomp, const amrex::Array4& avg_xmom, const amrex::Array4& avg_ymom, @@ -39,10 +41,11 @@ void AdvectionSrcForScalars (const amrex::Box& bx, const amrex::Array4& cell_prim, const amrex::Array4& src, const amrex::Array4& detJ, + const amrex::Real dt, const amrex::GpuArray& cellSize, const amrex::Array4& mf_m, const AdvType horiz_adv_type, const AdvType vert_adv_type, - const int use_terrain); + const bool use_terrain, const bool is_two_way_coupling); /** Compute advection tendencies for all components of momentum */ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, @@ -58,7 +61,7 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& mf_u, const amrex::Array4& mf_v, const AdvType horiz_adv_type, const AdvType vert_adv_type, - const int use_terrain, const int domhi_z); + const bool use_terrain, const int domhi_z); AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Advection/AdvectionSrcForMom.cpp b/Source/Advection/AdvectionSrcForMom.cpp index 2ca74102a..0af4657bd 100644 --- a/Source/Advection/AdvectionSrcForMom.cpp +++ b/Source/Advection/AdvectionSrcForMom.cpp @@ -51,7 +51,7 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, const Array4& mf_v, const AdvType horiz_adv_type, const AdvType vert_adv_type, - const int use_terrain, + const bool use_terrain, const int domhi_z) { BL_PROFILE_VAR("AdvectionSrcForMom", AdvectionSrcForMom); diff --git a/Source/Advection/AdvectionSrcForMom_N.H b/Source/Advection/AdvectionSrcForMom_N.H index a46807956..d20b504b6 100644 --- a/Source/Advection/AdvectionSrcForMom_N.H +++ b/Source/Advection/AdvectionSrcForMom_N.H @@ -42,22 +42,27 @@ AdvectionSrcForXMom_N (int i, int j, int k, amrex::Real interp_hi(0.), interp_lo(0.); rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0)); - rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0)); - interp_u_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo); + interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi); xflux_hi = rho_u_avg_hi * interp_hi; + + rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0)); + interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo); xflux_lo = rho_u_avg_lo * interp_lo; rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_v_inv(i,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0)); - rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0)); - interp_u_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo); + interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi); yflux_hi = rho_v_avg_hi * interp_hi; + + rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0)); + interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo); yflux_lo = rho_v_avg_lo * interp_lo; rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i-1, j, k+1)); - rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i-1, j, k )); - interp_u_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi); - interp_u_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo); + interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi); zflux_hi = rho_w_avg_hi * interp_hi; + + rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i-1, j, k )); + interp_u_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo); zflux_lo = rho_w_avg_lo * interp_lo; amrex::Real mfsq = 1 / (mf_u_inv(i,j,0) * mf_u_inv(i,j,0)); @@ -110,22 +115,27 @@ AdvectionSrcForYMom_N (int i, int j, int k, amrex::Real interp_hi(0.), interp_lo(0.); rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i+1, j-1, k) * mf_u_inv(i+1,j-1,0)); - rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0)); - interp_v_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo); + interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi); xflux_hi = rho_u_avg_hi * interp_hi; + + rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0)); + interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo); xflux_lo = rho_u_avg_lo * interp_lo; rho_v_avg_hi = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0)); - rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0)); - interp_v_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo); + interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi); yflux_hi = rho_v_avg_hi * interp_hi; + + rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0)); + interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo); yflux_lo = rho_v_avg_lo * interp_lo; rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i, j-1, k+1)); - rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i, j-1, k )); - interp_v_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi); - interp_v_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo); + interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi); zflux_hi = rho_w_avg_hi * interp_hi; + + rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i, j-1, k )); + interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,rho_w_avg_lo); zflux_lo = rho_w_avg_lo * interp_lo; amrex::Real mfsq = 1 / (mf_v_inv(i,j,0) * mf_v_inv(i,j,0)); @@ -186,15 +196,19 @@ AdvectionSrcForZMom_N (int i, int j, int k, amrex::Real interp_hi(0.), interp_lo(0.); rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0); - rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0); - interp_w_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo); + interp_w_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi); xflux_hi = rho_u_avg_hi * interp_hi; + + rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0); + interp_w_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo); xflux_lo = rho_u_avg_lo * interp_lo; rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0); - rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0); - interp_w_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo); + interp_w_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi); yflux_hi = rho_v_avg_hi * interp_hi; + + rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0); + interp_w_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo); yflux_lo = rho_v_avg_lo * interp_lo; // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1)); @@ -208,15 +222,15 @@ AdvectionSrcForZMom_N (int i, int j, int k, rho_w_avg_hi = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k+1)); if (k == domhi_z) { - interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd); + interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd); } else if (k == domhi_z-1 || k == 1) { if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) { - interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th); + interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th); } else { - interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,vert_adv_type); + interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,vert_adv_type); } } else { - interp_w_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi); + interp_w_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi); } zflux_hi = rho_w_avg_hi * interp_hi; } @@ -231,15 +245,15 @@ AdvectionSrcForZMom_N (int i, int j, int k, } else { rho_w_avg_lo = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k-1)); if (k == 1) { - interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd); + interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd); } else if (k == 2 || k == domhi_z) { if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) { - interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th); + interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th); } else { - interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type); + interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type); } } else { - interp_w_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo); + interp_w_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo); } zflux_lo = rho_w_avg_lo * interp_lo; } diff --git a/Source/Advection/AdvectionSrcForMom_T.H b/Source/Advection/AdvectionSrcForMom_T.H index a8b7ed79b..e61d22fce 100644 --- a/Source/Advection/AdvectionSrcForMom_T.H +++ b/Source/Advection/AdvectionSrcForMom_T.H @@ -47,9 +47,10 @@ AdvectionSrcForXMom_T (int i, int j, int k, // **************************************************************************************** rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i, j, k) * mf_u_inv(i ,j ,0)); - rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j ,0) + rho_u(i, j, k) * mf_u_inv(i ,j ,0)); + interp_u_h.InterpolateInX(i,j,k+1,0,interp_hi,rho_u_avg_hi); - interp_u_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo); + rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j ,0) + rho_u(i, j, k) * mf_u_inv(i ,j ,0)); + interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo); amrex::Real centFluxXXNext = rho_u_avg_hi * Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd) * interp_hi; amrex::Real centFluxXXPrev = rho_u_avg_lo * Compute_h_zeta_AtCellCenter(i-1,j ,k ,cellSizeInv,z_nd) * interp_lo; @@ -58,9 +59,10 @@ AdvectionSrcForXMom_T (int i, int j, int k, // Y-fluxes (at edges in k-direction) // **************************************************************************************** rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0)); - rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i ,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0)); + interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi); - interp_u_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo); + rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i ,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0)); + interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo); amrex::Real edgeFluxXYNext = rho_v_avg_hi * Compute_h_zeta_AtEdgeCenterK(i ,j+1,k ,cellSizeInv,z_nd) * interp_hi; amrex::Real edgeFluxXYPrev = rho_v_avg_lo * Compute_h_zeta_AtEdgeCenterK(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; @@ -71,8 +73,8 @@ AdvectionSrcForXMom_T (int i, int j, int k, Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i-1, j, k+1)); Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i-1, j, k )); - interp_u_v.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi); - interp_u_v.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo); + interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi); + interp_u_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo); amrex::Real edgeFluxXZNext = Omega_avg_hi * interp_hi; amrex::Real edgeFluxXZPrev = Omega_avg_lo * interp_lo; @@ -133,9 +135,10 @@ AdvectionSrcForYMom_T (int i, int j, int k, // x-fluxes (at edges in k-direction) // **************************************************************************************** rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i+1, j-1, k) * mf_u_inv(i+1,j-1,0)); - rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0)); + interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi); - interp_v_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo); + rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0)); + interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo); amrex::Real edgeFluxYXNext = rho_u_avg_hi * Compute_h_zeta_AtEdgeCenterK(i+1,j ,k ,cellSizeInv,z_nd) * interp_hi; amrex::Real edgeFluxYXPrev = rho_u_avg_lo * Compute_h_zeta_AtEdgeCenterK(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; @@ -144,9 +147,10 @@ AdvectionSrcForYMom_T (int i, int j, int k, // y-fluxes (at cell centers) // **************************************************************************************** rho_v_avg_hi = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0)); - rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0)); + interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi); - interp_v_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo); + rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0)); + interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo); amrex::Real centFluxYYNext = rho_v_avg_hi * Compute_h_zeta_AtCellCenter(i ,j ,k ,cellSizeInv,z_nd) * interp_hi; amrex::Real centFluxYYPrev = rho_v_avg_lo * Compute_h_zeta_AtCellCenter(i ,j-1,k ,cellSizeInv,z_nd) * interp_lo; @@ -157,8 +161,8 @@ AdvectionSrcForYMom_T (int i, int j, int k, Omega_avg_hi = 0.5 * (Omega(i, j, k+1) + Omega(i, j-1, k+1)); Omega_avg_lo = 0.5 * (Omega(i, j, k ) + Omega(i, j-1, k )); - interp_v_v.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi); - interp_v_v.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo); + interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi); + interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,Omega_avg_lo); amrex::Real edgeFluxYZNext = Omega_avg_hi * interp_hi; amrex::Real edgeFluxYZPrev = Omega_avg_lo * interp_lo; @@ -227,9 +231,10 @@ AdvectionSrcForZMom_T (int i, int j, int k, // x-fluxes (at edges in j-direction) // **************************************************************************************** rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0); - rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0); + interp_omega_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi); - interp_omega_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo); + rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0); + interp_omega_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo); amrex::Real edgeFluxZXNext = rho_u_avg_hi * Compute_h_zeta_AtEdgeCenterJ(i+1,j ,k ,cellSizeInv,z_nd) * interp_hi; amrex::Real edgeFluxZXPrev = rho_u_avg_lo * Compute_h_zeta_AtEdgeCenterJ(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; @@ -238,9 +243,10 @@ AdvectionSrcForZMom_T (int i, int j, int k, // y-fluxes (at edges in i-direction) // **************************************************************************************** rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0); - rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0); + interp_omega_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi); - interp_omega_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo); + rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0); + interp_omega_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo); amrex::Real edgeFluxZYNext = rho_v_avg_hi * Compute_h_zeta_AtEdgeCenterI(i ,j+1,k ,cellSizeInv,z_nd) * interp_hi; amrex::Real edgeFluxZYPrev = rho_v_avg_lo * Compute_h_zeta_AtEdgeCenterI(i ,j ,k ,cellSizeInv,z_nd) * interp_lo; @@ -260,15 +266,15 @@ AdvectionSrcForZMom_T (int i, int j, int k, centFluxZZNext *= w(i,j,k); } else { if (k == domhi_z) { - interp_omega_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi,AdvType::Centered_2nd); + interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,AdvType::Centered_2nd); } else if (k == domhi_z-1 || k == 1) { if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) { - interp_omega_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi,AdvType::Centered_4th); + interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,AdvType::Centered_4th); } else { - interp_omega_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi,vert_adv_type); + interp_omega_wall.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi,vert_adv_type); } } else { - interp_omega_v.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi); + interp_omega_v.InterpolateInZ(i,j,k+1,0,interp_hi,Omega_avg_hi); } centFluxZZNext *= interp_hi; } @@ -286,15 +292,15 @@ AdvectionSrcForZMom_T (int i, int j, int k, centFluxZZPrev *= w(i,j,k); } else { if (k == 1) { - interp_omega_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo,AdvType::Centered_2nd); + interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,AdvType::Centered_2nd); } else if (k == 2 || k == domhi_z) { if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) { - interp_omega_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo,AdvType::Centered_4th); + interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,AdvType::Centered_4th); } else { - interp_omega_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo,vert_adv_type); + interp_omega_wall.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo,vert_adv_type); } } else { - interp_omega_v.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo); + interp_omega_v.InterpolateInZ(i,j,k,0,interp_lo,Omega_avg_lo); } centFluxZZPrev *= interp_lo; } diff --git a/Source/Advection/AdvectionSrcForScalars.H b/Source/Advection/AdvectionSrcForScalars.H new file mode 100644 index 000000000..44e5865a3 --- /dev/null +++ b/Source/Advection/AdvectionSrcForScalars.H @@ -0,0 +1,101 @@ +#include +#include + +/** + * Wrapper function for computing the advective tendency w/ spatial order > 2. + */ +template +void +AdvectionSrcForScalarsWrapper_N(const amrex::Box& bx, + const int& ncomp, const int& icomp, + const amrex::GpuArray, AMREX_SPACEDIM> flx_arr, + const amrex::Array4& cell_prim, + const amrex::Array4& avg_xmom, + const amrex::Array4& avg_ymom, + const amrex::Array4& avg_zmom) +{ + // Instantiate structs for vert/horiz interp + InterpType_H interp_prim_h(cell_prim); + InterpType_V interp_prim_v(cell_prim); + + const amrex::Box xbx = amrex::surroundingNodes(bx,0); + const amrex::Box ybx = amrex::surroundingNodes(bx,1); + const amrex::Box zbx = amrex::surroundingNodes(bx,2); + + amrex::ParallelFor(xbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int cons_index = icomp + n; + const int prim_index = cons_index - 1; + + amrex::Real interpx(0.); + + interp_prim_h.InterpolateInX(i,j,k,prim_index,interpx,avg_xmom(i ,j ,k )); + (flx_arr[0])(i,j,k,n) = avg_xmom(i,j,k) * interpx; + }); + amrex::ParallelFor(ybx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int cons_index = icomp + n; + const int prim_index = cons_index - 1; + + amrex::Real interpy(0.); + interp_prim_h.InterpolateInY(i,j,k,prim_index,interpy,avg_ymom(i ,j ,k )); + + (flx_arr[1])(i,j,k,n) = avg_ymom(i,j,k) * interpy; + }); + amrex::ParallelFor(zbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int cons_index = icomp + n; + const int prim_index = cons_index - 1; + + amrex::Real interpz(0.); + + interp_prim_v.InterpolateInZ(i,j,k,prim_index,interpz,avg_zmom(i ,j ,k )); + + (flx_arr[2])(i,j,k,n) = avg_zmom(i,j,k) * interpz; + }); +} + +/** + * Wrapper function for templating the vertical advective tendency w/ spatial order > 2. + */ +template +void +AdvectionSrcForScalarsVert_N(const amrex::Box& bx, + const int& ncomp, const int& icomp, + const amrex::GpuArray, AMREX_SPACEDIM> flx_arr, + const amrex::Array4& cell_prim, + const amrex::Array4& avg_xmom, + const amrex::Array4& avg_ymom, + const amrex::Array4& avg_zmom, + const AdvType vert_adv_type) +{ + switch(vert_adv_type) { + case AdvType::Centered_2nd: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Upwind_3rd: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Centered_4th: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Upwind_5th: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Centered_6th: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, + flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + default: + AMREX_ASSERT_WITH_MESSAGE(false, "Unknown vertical advection scheme!"); + } +} diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index ae260a59c..4be573815 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -1,8 +1,10 @@ #include #include #include +#include #include #include +#include using namespace amrex; @@ -51,7 +53,7 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, const Array4& mf_v, const AdvType horiz_adv_type, const AdvType vert_adv_type, - const int use_terrain) + const bool use_terrain) { BL_PROFILE_VAR("AdvectionSrcForRhoAndTheta", AdvectionSrcForRhoAndTheta); auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; @@ -100,37 +102,43 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, }); // Template higher order methods } else { - if (horiz_adv_type == AdvType::Centered_2nd) { + switch(horiz_adv_type) { + case AdvType::Centered_2nd: AdvectionSrcForRhoThetaVert_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Upwind_3rd) { + break; + case AdvType::Upwind_3rd: AdvectionSrcForRhoThetaVert_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Centered_4th) { + break; + case AdvType::Centered_4th: AdvectionSrcForRhoThetaVert_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Upwind_5th) { + break; + case AdvType::Upwind_5th: AdvectionSrcForRhoThetaVert_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Centered_6th) { + break; + case AdvType::Centered_6th: AdvectionSrcForRhoThetaVert_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else { + break; + default: AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } } @@ -188,37 +196,43 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, }); // Template higher order methods (horizontal first) } else { - if (horiz_adv_type == AdvType::Centered_2nd) { + switch(horiz_adv_type) { + case AdvType::Centered_2nd: AdvectionSrcForRhoThetaVert_T(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Upwind_3rd) { + break; + case AdvType::Upwind_3rd: AdvectionSrcForRhoThetaVert_T(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Centered_4th) { + break; + case AdvType::Centered_4th: AdvectionSrcForRhoThetaVert_T(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Upwind_5th) { + break; + case AdvType::Upwind_5th: AdvectionSrcForRhoThetaVert_T(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else if (horiz_adv_type == AdvType::Centered_6th) { + break; + case AdvType::Centered_6th: AdvectionSrcForRhoThetaVert_T(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, Omega, avg_xmom, avg_ymom, avg_zmom, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, vert_adv_type); - } else { + break; + default: AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } } @@ -248,98 +262,141 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, */ void -AdvectionSrcForScalars (const Box& bx, const int icomp, const int ncomp, +AdvectionSrcForScalars (int level, int finest_level, + const MFIter& mfi, + const Box& bx, const int icomp, const int ncomp, const Array4& avg_xmom, const Array4& avg_ymom, const Array4& avg_zmom, const Array4& cell_prim, const Array4& advectionSrc, const Array4& detJ, + const Real dt, const GpuArray& cellSizeInv, const Array4& mf_m, const AdvType horiz_adv_type, const AdvType vert_adv_type, - const int use_terrain) + const bool use_terrain, + const bool two_way_coupling) { BL_PROFILE_VAR("AdvectionSrcForScalars", AdvectionSrcForScalars); - auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; + auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; + + amrex::GpuArray flux; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + const Box& efbx = amrex::surroundingNodes(bx, dir); + flux[dir].resize(efbx, ncomp, amrex::The_Async_Arena()); + flux[dir].setVal(0.); + } + + const GpuArray, AMREX_SPACEDIM> + flx_arr{{AMREX_D_DECL(flux[0].array(), flux[1].array(), flux[2].array())}}; + + const Box xbx = surroundingNodes(bx,0); + const Box ybx = surroundingNodes(bx,1); + const Box zbx = surroundingNodes(bx,2); // Inline with 2nd order for efficiency + // NOTE: we don't need to weight avg_xmom, avg_ymom, avg_zmom with terrain metrics + // because that was done when they were constructed in AdvectionSrcForRhoAndTheta if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) { - amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + amrex::ParallelFor(xbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; - - // NOTE: we don't need to weight avg_xmom, avg_ymom, avg_zmom with terrain metrics - // because that was done when they were constructed in AdvectionSrcForRhoAndTheta - const int cons_index = icomp + n; const int prim_index = cons_index - 1; - - Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); - - advectionSrc(i,j,k,cons_index) = - 0.5 * invdetJ * ( - ( avg_xmom(i+1,j,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i+1,j,k,prim_index)) - - avg_xmom(i ,j,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i-1,j,k,prim_index)) ) * dxInv * mfsq + - ( avg_ymom(i,j+1,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i,j+1,k,prim_index)) - - avg_ymom(i,j ,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i,j-1,k,prim_index)) ) * dyInv * mfsq + - ( avg_zmom(i,j,k+1) * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k+1,prim_index)) - - avg_zmom(i,j,k ) * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k-1,prim_index)) ) * dzInv); + (flx_arr[0])(i,j,k,n) = 0.5 * avg_xmom(i,j,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i-1,j,k,prim_index)); + }); + amrex::ParallelFor(ybx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int cons_index = icomp + n; + const int prim_index = cons_index - 1; + (flx_arr[1])(i,j,k,n) = 0.5 * avg_ymom(i,j,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i,j-1,k,prim_index)); + }); + amrex::ParallelFor(zbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int cons_index = icomp + n; + const int prim_index = cons_index - 1; + (flx_arr[2])(i,j,k,n) = 0.5 * avg_zmom(i,j,k) * (cell_prim(i,j,k,prim_index) + cell_prim(i,j,k-1,prim_index)); }); + // Template higher order methods (horizontal first) } else { - if (horiz_adv_type == AdvType::Centered_2nd) { - AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_adv_type); - } else if (horiz_adv_type == AdvType::Upwind_3rd) { - AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_adv_type); - } else if (horiz_adv_type == AdvType::Centered_4th) { - AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_adv_type); - } else if (horiz_adv_type == AdvType::Upwind_5th) { - AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_adv_type); - } else if (horiz_adv_type == AdvType::Centered_6th) { - AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_adv_type); - } else if (horiz_adv_type == AdvType::Weno_3) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (horiz_adv_type == AdvType::Weno_5) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (horiz_adv_type == AdvType::Weno_3Z) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (horiz_adv_type == AdvType::Weno_3MZQ) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (horiz_adv_type == AdvType::Weno_5Z) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else { + switch(horiz_adv_type) { + case AdvType::Centered_2nd: + AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + break; + case AdvType::Upwind_3rd: + AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + break; + case AdvType::Centered_4th: + AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + break; + case AdvType::Upwind_5th: + AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + break; + case AdvType::Centered_6th: + AdvectionSrcForScalarsVert_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom, vert_adv_type); + break; + case AdvType::Weno_3: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Weno_5: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Weno_3Z: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Weno_3MZQ: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + case AdvType::Weno_5Z: + AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, flx_arr, cell_prim, + avg_xmom, avg_ymom, avg_zmom); + break; + default: AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } } + + amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; + + Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); + + const int cons_index = icomp + n; + advectionSrc(i,j,k,cons_index) = - invdetJ * ( + ( (flx_arr[0])(i+1,j,k,n) - (flx_arr[0])(i ,j,k,n) ) * dxInv * mfsq + + ( (flx_arr[1])(i,j+1,k,n) - (flx_arr[1])(i,j ,k,n) ) * dyInv * mfsq + + ( (flx_arr[2])(i,j,k+1,n) - (flx_arr[2])(i,j,k ,n) ) * dzInv ); + }); + +#if 0 + Real fac_for_reflux = 1.0; + + std::array dxD = {{dx, dy, dz}}; + const Real* dxDp = &(dxD[0]); + + if (two_way_coupling) { + if (level < finest_level) { + getAdvFluxReg(level + 1).CrseAdd(mfi, + {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}}, + dxDp, fac_for_reflux*dt, amrex::RunOn::Device); + } + if (level > 0) { + getAdvFluxReg(level).FineAdd(mfi, + {{AMREX_D_DECL(&(flux[0]), &(flux[1]), &(flux[2]))}}, + dxDp, fac_for_reflux*dt, amrex::RunOn::Device); + } + } // two-way coupling +#endif } diff --git a/Source/Advection/AdvectionSrcForState_N.H b/Source/Advection/AdvectionSrcForState_N.H index e594fa1af..3758df06f 100644 --- a/Source/Advection/AdvectionSrcForState_N.H +++ b/Source/Advection/AdvectionSrcForState_N.H @@ -58,11 +58,14 @@ AdvectionSrcForRhoThetaWrapper_N(const amrex::Box& bx, amrex::Real interpy_hi(0.), interpy_lo(0.); amrex::Real interpz_hi(0.), interpz_lo(0.); - interp_prim_h.InterpolateInX(i,j,k,prim_index,interpx_hi,interpx_lo,rho_u(i+1,j ,k ),rho_u(i ,j ,k )); - interp_prim_h.InterpolateInY(i,j,k,prim_index,interpy_hi,interpy_lo,rho_v(i ,j+1,k ),rho_v(i ,j ,k )); + interp_prim_h.InterpolateInX(i+1,j,k,prim_index,interpx_hi,rho_u(i+1,j ,k )); + interp_prim_h.InterpolateInX(i ,j,k,prim_index,interpx_lo,rho_u(i ,j ,k )); - interp_prim_v.InterpolateInZ_hi(i,j,k,prim_index,interpz_hi,rho_w(i ,j ,k+1)); - interp_prim_v.InterpolateInZ_lo(i,j,k,prim_index,interpz_lo,rho_w(i ,j ,k )); + interp_prim_h.InterpolateInY(i,j+1,k,prim_index,interpy_hi,rho_v(i ,j+1,k )); + interp_prim_h.InterpolateInY(i,j ,k,prim_index,interpy_lo,rho_v(i ,j ,k )); + + interp_prim_v.InterpolateInZ(i,j,k+1,prim_index,interpz_hi,rho_w(i ,j ,k+1)); + interp_prim_v.InterpolateInZ(i,j,k ,prim_index,interpz_lo,rho_w(i ,j ,k )); advectionSrc(i,j,k,1) = -( ( xflux_hi * interpx_hi - xflux_lo * interpx_lo ) * dxInv * mfsq + @@ -93,128 +96,39 @@ AdvectionSrcForRhoThetaVert_N(const amrex::Box& bx, const amrex::Array4& mf_v, const AdvType vert_adv_type) { - if (vert_adv_type == AdvType::Centered_2nd) { + switch(vert_adv_type) { + case AdvType::Centered_2nd: AdvectionSrcForRhoThetaWrapper_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, rho_w, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v); - } else if (vert_adv_type == AdvType::Upwind_3rd) { + break; + case AdvType::Upwind_3rd: AdvectionSrcForRhoThetaWrapper_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, rho_w, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v); - } else if (vert_adv_type == AdvType::Centered_4th) { + break; + case AdvType::Centered_4th: AdvectionSrcForRhoThetaWrapper_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, rho_w, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v); - } else if (vert_adv_type == AdvType::Upwind_5th) { + break; + case AdvType::Upwind_5th: AdvectionSrcForRhoThetaWrapper_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, rho_w, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v); - } else if (vert_adv_type == AdvType::Centered_6th) { + break; + case AdvType::Centered_6th: AdvectionSrcForRhoThetaWrapper_N(bx, vbx_hi, fac, advectionSrc, cell_prim, rho_u, rho_v, rho_w, avg_xmom, avg_ymom, avg_zmom, cellSizeInv, mf_m, mf_u, mf_v); - } else { + break; + default: AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } } -/** - * Wrapper function for computing the advective tendency w/ spatial order > 2. - */ -template -void -AdvectionSrcForScalarsWrapper_N(const amrex::Box& bx, - const int& ncomp, const int& icomp, const int& use_terrain, - const amrex::Array4< amrex::Real>& advectionSrc, - const amrex::Array4& cell_prim, - const amrex::Array4& avg_xmom, - const amrex::Array4& avg_ymom, - const amrex::Array4& avg_zmom, - const amrex::Array4& detJ, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_m) -{ - // Instantiate structs for vert/horiz interp - InterpType_H interp_prim_h(cell_prim); - InterpType_V interp_prim_v(cell_prim); - - auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; - amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - amrex::Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; - amrex::Real mfsq = mf_m(i,j,0) * mf_m(i,j,0); - - // NOTE: we don't need to weight avg_xmom, avg_ymom, avg_zmom with terrain metrics - // because that was done when they were constructed in AdvectionSrcForRhoAndTheta - - const int cons_index = icomp + n; - const int prim_index = cons_index - 1; - - amrex::Real interpx_hi(0.), interpx_lo(0.); - amrex::Real interpy_hi(0.), interpy_lo(0.); - amrex::Real interpz_hi(0.), interpz_lo(0.); - - interp_prim_h.InterpolateInX(i,j,k,prim_index,interpx_hi,interpx_lo,avg_xmom(i+1,j ,k ),avg_xmom(i ,j ,k )); - interp_prim_h.InterpolateInY(i,j,k,prim_index,interpy_hi,interpy_lo,avg_ymom(i ,j+1,k ),avg_ymom(i ,j ,k )); - - interp_prim_v.InterpolateInZ_hi(i,j,k,prim_index,interpz_hi,avg_zmom(i ,j ,k+1)); - interp_prim_v.InterpolateInZ_lo(i,j,k,prim_index,interpz_lo,avg_zmom(i ,j ,k )); - - advectionSrc(i,j,k,cons_index) = - invdetJ * ( - ( avg_xmom(i+1,j ,k ) * interpx_hi - avg_xmom(i ,j ,k ) * interpx_lo ) * dxInv * mfsq + - ( avg_ymom(i ,j+1,k ) * interpy_hi - avg_ymom(i ,j ,k ) * interpy_lo ) * dyInv * mfsq + - ( avg_zmom(i ,j ,k+1) * interpz_hi - avg_zmom(i ,j ,k ) * interpz_lo ) * dzInv); - }); -} - -/** - * Wrapper function for templating the vertical advective tendency w/ spatial order > 2. - */ -template -void -AdvectionSrcForScalarsVert_N(const amrex::Box& bx, - const int& ncomp, const int& icomp, const int& use_terrain, - const amrex::Array4< amrex::Real>& advectionSrc, - const amrex::Array4& cell_prim, - const amrex::Array4& avg_xmom, - const amrex::Array4& avg_ymom, - const amrex::Array4& avg_zmom, - const amrex::Array4& detJ, - const amrex::GpuArray& cellSizeInv, - const amrex::Array4& mf_m, - const AdvType vert_adv_type) -{ - if (vert_adv_type == AdvType::Centered_2nd) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_adv_type == AdvType::Upwind_3rd) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_adv_type == AdvType::Centered_4th) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_adv_type == AdvType::Upwind_5th) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_adv_type == AdvType::Centered_6th) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else { - AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); - } -} diff --git a/Source/Advection/AdvectionSrcForState_T.H b/Source/Advection/AdvectionSrcForState_T.H index b04fb42b8..79ee4e0a9 100644 --- a/Source/Advection/AdvectionSrcForState_T.H +++ b/Source/Advection/AdvectionSrcForState_T.H @@ -73,11 +73,14 @@ AdvectionSrcForRhoThetaWrapper_T(const amrex::Box& bx, amrex::Real interpy_hi(0.), interpy_lo(0.); amrex::Real interpz_hi(0.), interpz_lo(0.); - interp_prim_h.InterpolateInX(i,j,k,prim_index,interpx_hi,interpx_lo,rho_u(i+1,j ,k ),rho_u(i ,j ,k )); - interp_prim_h.InterpolateInY(i,j,k,prim_index,interpy_hi,interpy_lo,rho_v(i ,j+1,k ),rho_v(i ,j ,k )); + interp_prim_h.InterpolateInX(i+1,j,k,prim_index,interpx_hi,rho_u(i+1,j ,k )); + interp_prim_h.InterpolateInX(i ,j,k,prim_index,interpx_lo,rho_u(i ,j ,k )); - interp_prim_v.InterpolateInZ_hi(i,j,k,prim_index,interpz_hi,Omega(i ,j ,k+1)); - interp_prim_v.InterpolateInZ_lo(i,j,k,prim_index,interpz_lo,Omega(i ,j ,k )); + interp_prim_h.InterpolateInY(i,j+1,k,prim_index,interpy_hi,rho_v(i ,j+1,k )); + interp_prim_h.InterpolateInY(i,j ,k,prim_index,interpy_lo,rho_v(i ,j ,k )); + + interp_prim_v.InterpolateInZ(i,j,k+1,prim_index,interpz_hi,Omega(i ,j ,k+1)); + interp_prim_v.InterpolateInZ(i,j,k ,prim_index,interpz_lo,Omega(i ,j ,k )); advectionSrc(i,j,k,1) = - invdetJ * ( ( xflux_hi * interpx_hi - xflux_lo * interpx_lo ) * dxInv * mfsq + diff --git a/Source/Advection/Make.package b/Source/Advection/Make.package index 4a77fb744..d2bcefa7f 100644 --- a/Source/Advection/Make.package +++ b/Source/Advection/Make.package @@ -6,4 +6,5 @@ CEXE_headers += AdvectionSrcForMom_N.H CEXE_headers += AdvectionSrcForMom_T.H CEXE_headers += AdvectionSrcForState_N.H CEXE_headers += AdvectionSrcForState_T.H +CEXE_headers += AdvectionSrcForScalars.H CEXE_headers += Interpolation.H diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index b00ea1ba7..e320c6c42 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -159,7 +159,7 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, const Array4& z_phys_nd, const GpuArray dxInv, int bccomp_u, int bccomp_v, int bccomp_w, - int terrain_type) + TerrainType terrain_type) { BL_PROFILE_VAR("impose_vertical_zvel_bcs()",impose_vertical_zvel_bcs); const auto& dom_lo = amrex::lbound(domain); @@ -198,7 +198,7 @@ void ERFPhysBCFunct::impose_vertical_zvel_bcs (const Array4& dest_arr, const amrex::BCRec* bc_ptr_w = bcrs_w_d.data(); bool l_use_terrain = (m_z_phys_nd != nullptr); - bool l_moving_terrain = (terrain_type == 1); + bool l_moving_terrain = (terrain_type == TerrainType::Moving); GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 7d2a4e4ed..bf750a410 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -73,7 +73,7 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) } // var_idx // Coarse-Fine set region - if (lev>0 && coupling_type=="OneWay" && cf_set_width>0 && fillset) { + if (lev>0 && solverChoice.coupling_type == CouplingType::OneWay && cf_set_width>0 && fillset) { FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type); FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type); FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type); @@ -226,7 +226,7 @@ ERF::FillIntermediatePatch (int lev, Real time, } // var_idx // Coarse-Fine set region - if (lev>0 && coupling_type=="OneWay" && cf_set_width>0) { + if (lev>0 && solverChoice.coupling_type == CouplingType::OneWay && cf_set_width>0) { if (cons_only) { FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type); } else { diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index 4c4a7ad55..c7790e0f1 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -25,7 +25,7 @@ public: ERFPhysBCFunct (const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, - const int& terrain_type, + const TerrainType& terrain_type, amrex::Array,AMREX_SPACEDIM+NVAR> bc_extdir_vals, amrex::Array,AMREX_SPACEDIM+NVAR> bc_neumann_vals, std::unique_ptr& z_phys_nd, @@ -92,7 +92,7 @@ public: const amrex::Array4& z_nd, const amrex::GpuArray dxInv, int bccomp_u, int bccomp_v, int bccomp_w, - int terrain_type); + TerrainType terrain_type); void impose_lateral_cons_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, @@ -108,7 +108,7 @@ private: amrex::Geometry m_geom; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; - int m_terrain_type; + TerrainType m_terrain_type; amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals; amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_neumann_vals; amrex::MultiFab* m_z_phys_nd; diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index e8941c11e..abf87a276 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -20,6 +20,14 @@ enum struct ABLDriverType { None, PressureGradient, GeostrophicWind }; +enum struct CouplingType { + OneWay, TwoWay +}; + +enum struct TerrainType { + Static, Moving +}; + enum struct Coord { x, y, z }; @@ -59,11 +67,20 @@ struct SolverChoice { } // Is the terrain static or moving? - pp.query("terrain_type", terrain_type); + static std::string terrain_type_string = "Static"; + pp.query("terrain_type",terrain_type_string); + if (terrain_type_string == "Moving" || terrain_type_string == "moving") { + terrain_type = TerrainType::Moving; + } else if (terrain_type_string == "Static" || terrain_type_string == "static") { + terrain_type = TerrainType::Static; + } else { + amrex::Abort("terrain_type can be either Moving/moving or Static/static"); + } // Use lagged_delta_rt in the fast integrator? pp.query("use_lagged_delta_rt", use_lagged_delta_rt); - if (!use_lagged_delta_rt && !(terrain_type == 1)) { + + if (!use_lagged_delta_rt && !(terrain_type == TerrainType::Moving)) { amrex::Error("Can't turn off lagged_delta_rt when terrain not moving"); } @@ -171,6 +188,15 @@ struct SolverChoice { } } } + + // Which type of refinement + static std::string coupling_type_string = "OneWay"; + pp.query("coupling_type",coupling_type_string); + if (coupling_type_string == "TwoWay") { + coupling_type = CouplingType::TwoWay; + } else { + coupling_type = CouplingType::OneWay; + } } void display() @@ -183,6 +209,18 @@ struct SolverChoice { amrex::Print() << "use_rayleigh_damping : " << use_rayleigh_damping << std::endl; amrex::Print() << "use_gravity : " << use_gravity << std::endl; + if (coupling_type == CouplingType::TwoWay) { + amrex::Print() << "Using two-way coupling " << std::endl; + } else { + amrex::Print() << "Using one-way coupling " << std::endl; + } + + if (terrain_type == TerrainType::Static) { + amrex::Print() << "Using static terrain " << std::endl; + } else { + amrex::Print() << "Using moving terrain " << std::endl; + } + if (abl_driver_type == ABLDriverType::None) { amrex::Print() << "ABL Driver Type: " << "None" << std::endl; amrex::Print() << "No ABL driver selected " << std::endl; @@ -259,7 +297,6 @@ struct SolverChoice { bool test_mapfactor = false; bool use_terrain = false; - int terrain_type = 0; #ifdef ERF_USE_MOISTURE int buoyancy_type = 2; // uses Tprime #else @@ -306,6 +343,9 @@ struct SolverChoice { bool use_NumDiff{false}; amrex::Real NumDiffCoeff{0.}; + CouplingType coupling_type; + TerrainType terrain_type; + ABLDriverType abl_driver_type; amrex::GpuArray abl_pressure_grad; amrex::GpuArray abl_geo_forcing; diff --git a/Source/ERF.H b/Source/ERF.H index 340a4edbb..e7e5d533d 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -19,7 +19,7 @@ #include #include #include -#include +#include #include #ifdef AMREX_MEM_PROFILING @@ -574,7 +574,7 @@ private: amrex::Vector base_state_new; // array of flux registers - amrex::Vector flux_registers; + amrex::Vector advflux_reg; // A BCRec is essentially a 2*DIM integer array storing the boundary // condition type at each lo/hi walls in each direction. We have one BCRec @@ -680,9 +680,6 @@ private: // flag to turn tracer particle generation on/off static bool use_tracer_particles; - // mesh refinement - static std::string coupling_type; - // Diagnostic output interval static int sum_interval; static amrex::Real sum_per; @@ -820,10 +817,9 @@ private: } AMREX_FORCE_INLINE - amrex::FluxRegister& - get_flux_reg (int lev) + amrex::YAFluxRegister& getAdvFluxReg (int lev) { - return *flux_registers[lev]; + return *advflux_reg[lev]; } AMREX_FORCE_INLINE diff --git a/Source/ERF.cpp b/Source/ERF.cpp index ac2641b92..6ecdc9384 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -42,9 +42,6 @@ bool ERF::use_tracer_particles = false; amrex::Vector ERF::tracer_particle_varnames = {AMREX_D_DECL("xvel", "yvel", "zvel")}; #endif -// Type of mesh refinement algorithm -std::string ERF::coupling_type = "OneWay"; - // Dictate verbosity in screen output int ERF::verbose = 0; @@ -169,7 +166,7 @@ ERF::ERF () mri_integrator_mem.resize(nlevs_max); physbcs.resize(nlevs_max); - flux_registers.resize(nlevs_max); + advflux_reg.resize(nlevs_max); // Stresses Tau11_lev.resize(nlevs_max); Tau22_lev.resize(nlevs_max); Tau33_lev.resize(nlevs_max); @@ -315,12 +312,12 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) { BL_PROFILE("ERF::post_timestep()"); - if (coupling_type == "TwoWay") + if (solverChoice.coupling_type == CouplingType::TwoWay) { for (int lev = finest_level-1; lev >= 0; lev--) { // This call refluxes from the lev/lev+1 interface onto lev - get_flux_reg(lev+1).Reflux(vars_new[lev][Vars::cons],1.0, 0, 0, NVAR, geom[lev]); + getAdvFluxReg(lev+1).Reflux(vars_new[lev][Vars::cons], 0, 0, NVAR); // We need to do this before anything else because refluxing changes the // values of coarse cells underneath fine grids with the assumption they'll @@ -367,7 +364,7 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) } // Moving terrain - if ( solverChoice.use_terrain && (solverChoice.terrain_type == 1) ) + if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) { for (int lev = finest_level; lev >= 0; lev--) { @@ -403,8 +400,8 @@ ERF::InitData () } } - if (!solverChoice.use_terrain && solverChoice.terrain_type != 0) { - amrex::Abort("We do not allow terrain_type != 0 with use_terrain = false"); + if (!solverChoice.use_terrain && solverChoice.terrain_type != TerrainType::Static) { + amrex::Abort("We do not allow non-static terrain_type with use_terrain = false"); } last_plot_file_step_1 = -1; @@ -485,7 +482,7 @@ ERF::InitData () } - if (coupling_type == "TwoWay") { + if (solverChoice.coupling_type == CouplingType::TwoWay) { AverageDown(); } @@ -525,11 +522,14 @@ ERF::InitData () } // Initialize flux registers (whether we start from scratch or restart) - if (coupling_type == "TwoWay") { - flux_registers[0] = 0; + if (solverChoice.coupling_type == CouplingType::TwoWay) { + advflux_reg[0] = nullptr; for (int lev = 1; lev <= finest_level; lev++) { - flux_registers[lev] = new FluxRegister(grids[lev], dmap[lev], ref_ratio[lev-1], lev, NVAR); + advflux_reg[lev] = new YAFluxRegister(grids[lev], grids[lev-1], + dmap[lev], dmap[lev-1], + geom[lev], geom[lev-1], + ref_ratio[lev-1], lev, NVAR); } } @@ -634,7 +634,7 @@ ERF::InitData () { // NOTE: we must set up the FillPatcher object before calling // FillPatch at a fine level - if (coupling_type=="OneWay" && cf_width>0 && lev>0) { + if (solverChoice.coupling_type== CouplingType::OneWay && cf_width>0 && lev>0) { Construct_ERFFillPatchers(lev); Register_ERFFillPatchers(lev); } @@ -649,7 +649,7 @@ ERF::InitData () base_state[lev].FillBoundary(geom[lev].periodicity()); // For moving terrain only - if (solverChoice.terrain_type > 0) { + if (solverChoice.terrain_type != TerrainType::Static) { MultiFab::Copy(base_state_new[lev],base_state[lev],0,0,3,1); base_state_new[lev].FillBoundary(geom[lev].periodicity()); } @@ -936,10 +936,6 @@ ERF::ReadParameters () AMREX_ALWAYS_ASSERT(cfl > 0. || fixed_dt > 0.); - // Mesh refinement - pp.query("coupling_type",coupling_type); - AMREX_ALWAYS_ASSERT( (coupling_type == "OneWay") || (coupling_type == "TwoWay") ); - // How to initialize pp.query("init_type",init_type); if (!init_type.empty() && @@ -953,12 +949,8 @@ ERF::ReadParameters () } // No moving terrain with init real - if (init_type == "real") { - int terr_type(0); - pp.query("terrain_type",terr_type); - if (terr_type) { - amrex::Abort("Moving terrain is not supported with init real"); - } + if (init_type == "real" && solverChoice.terrain_type != TerrainType::Static) { + amrex::Abort("Moving terrain is not supported with init real"); } // We use this to keep track of how many boxes we read in from WRF initialization @@ -1072,7 +1064,7 @@ ERF::MakeHorizontalAverages () int lev = 0; // First, average down all levels (if doing two-way coupling) - if (coupling_type == "TwoWay") { + if (solverChoice.coupling_type == CouplingType::TwoWay) { AverageDown(); } @@ -1204,7 +1196,7 @@ ERF::MakeDiagnosticAverage (Vector& h_havg, MultiFab& S, int n) void ERF::AverageDown () { - AMREX_ALWAYS_ASSERT(coupling_type == "TwoWay"); + AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); for (int lev = finest_level-1; lev >= 0; --lev) { AverageDownTo(lev); @@ -1215,7 +1207,7 @@ ERF::AverageDown () void ERF::AverageDownTo (int crse_lev) // NOLINT { - AMREX_ALWAYS_ASSERT(coupling_type == "TwoWay"); + AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); for (int var_idx = 0; var_idx < Vars::NumTypes; ++var_idx) { const BoxArray& ba(vars_new[crse_lev][var_idx].boxArray()); if (ba[0].type() == IntVect::TheZeroVector()) @@ -1376,7 +1368,7 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, nbx = convert(domain_p[0],IntVect(0,0,1)); domain_p.push_back(nbx); - flux_registers.resize(nlevs_max); + advflux_reg.resize(nlevs_max); // Stresses Tau11_lev.resize(nlevs_max); Tau22_lev.resize(nlevs_max); Tau33_lev.resize(nlevs_max); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 2ada7d831..21e635c46 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -105,7 +105,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, base_state[lev].define(ba,dm,3,1); base_state[lev].setVal(0.); - if (solverChoice.use_terrain && solverChoice.terrain_type > 0) { + if (solverChoice.use_terrain && solverChoice.terrain_type != TerrainType::Static) { base_state_new[lev].define(ba,dm,3,1); base_state_new[lev].setVal(0.); } @@ -114,7 +114,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, z_phys_cc[lev] = std::make_unique(ba,dm,1,1); detJ_cc[lev] = std::make_unique(ba,dm,1,1); - if (solverChoice.terrain_type > 0) { + if (solverChoice.terrain_type != TerrainType::Static) { detJ_cc_new[lev] = std::make_unique(ba,dm,1,1); detJ_cc_src[lev] = std::make_unique(ba,dm,1,1); z_t_rk[lev] = std::make_unique( convert(ba, IntVect(0,0,1)), dm, 1, 1 ); @@ -126,7 +126,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, // We need this to be one greater than the ghost cells to handle levels > 0 int ngrow = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 2; z_phys_nd[lev] = std::make_unique(ba_nd,dm,1,IntVect(ngrow,ngrow,1)); - if (solverChoice.terrain_type > 0) { + if (solverChoice.terrain_type != TerrainType::Static) { z_phys_nd_new[lev] = std::make_unique(ba_nd,dm,1,IntVect(ngrow,ngrow,1)); z_phys_nd_src[lev] = std::make_unique(ba_nd,dm,1,IntVect(ngrow,ngrow,1)); } diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 207440853..621bc96d7 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -119,14 +119,14 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ // Update the dycore advance_dycore(lev, - cons_mf, S_new, - U_old, V_old, W_old, - U_new, V_new, W_new, - rU_old[lev], rV_old[lev], rW_old[lev], - rU_new[lev], rV_new[lev], rW_new[lev], - rU_crse, rV_crse, rW_crse, - source, buoyancy, - Geom(lev), dt_lev, time, &ifr); + cons_mf, S_new, + U_old, V_old, W_old, + U_new, V_new, W_new, + rU_old[lev], rV_old[lev], rW_old[lev], + rU_new[lev], rV_new[lev], rW_new[lev], + rU_crse, rV_crse, rW_crse, + source, buoyancy, + Geom(lev), dt_lev, time, &ifr); #if defined(ERF_USE_MOISTURE) // Update the microphysics diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 5d5361c20..5244dcefb 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -36,7 +36,7 @@ ERF::timeStep (int lev, Real time, int iteration) // NOTE: Def & Reg index lev backwards (so we add 1 here) // Redefine & register the ERFFillpatcher objects - if (coupling_type=="OneWay" && cf_width>0) { + if (solverChoice.coupling_type == CouplingType::OneWay && cf_width>0) { Define_ERFFillPatchers(lev+1); Register_ERFFillPatchers(lev+1); if (lev < max_level-1) { @@ -87,7 +87,7 @@ ERF::timeStep (int lev, Real time, int iteration) timeStep(lev+1, time+(i-1)*dt[lev+1], i); } - if (coupling_type == "TwoWay") { + if (solverChoice.coupling_type == CouplingType::TwoWay) { AverageDownTo(lev); // average lev+1 down to lev } } diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index e07bc1881..25a080309 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -314,7 +314,7 @@ void ERF::advance_dycore(int level, mri_integrator.advance(state_old, state_new, old_time, dt_advance); // Register coarse data for coarse-fine fill - if (level0) { + if (level0) { FPr_c[level].RegisterCoarseData({&cons_old, &cons_new}, {old_time, old_time + dt_advance}); FPr_u[level].RegisterCoarseData({&xvel_old, &xvel_new}, {old_time, old_time + dt_advance}); FPr_v[level].RegisterCoarseData({&yvel_old, &yvel_new}, {old_time, old_time + dt_advance}); diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 32579dbfb..51c10dd7b 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -44,7 +44,7 @@ using namespace amrex; * @param[in] mapfac_v map factor at y-faces */ -void erf_slow_rhs_post (int level, +void erf_slow_rhs_post (int level, int finest_level, int nrk, Real dt, Vector& S_rhs, @@ -93,8 +93,9 @@ void erf_slow_rhs_post (int level, const MultiFab* t_mean_mf = nullptr; if (most) t_mean_mf = most->get_mac_avg(0,2); - const bool l_use_terrain = solverChoice.use_terrain; - const bool l_moving_terrain = (solverChoice.terrain_type == 1); + const bool l_use_terrain = solverChoice.use_terrain; + const bool l_two_way_coupling = (solverChoice.coupling_type == CouplingType::TwoWay); + const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving); if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain); const bool l_use_ndiff = solverChoice.use_NumDiff; @@ -250,27 +251,27 @@ void erf_slow_rhs_post (int level, if (l_use_deardorff) { start_comp = RhoKE_comp; num_comp = 1; - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + AdvectionSrcForScalars(level, finest_level, mfi, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dt, dxInv, mf_m, horiz_adv_type, vert_adv_type, - l_use_terrain); + l_use_terrain, l_two_way_coupling); } if (l_use_QKE) { start_comp = RhoQKE_comp; num_comp = 1; - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + AdvectionSrcForScalars(level, finest_level, mfi, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dt, dxInv, mf_m, horiz_adv_type, vert_adv_type, - l_use_terrain); + l_use_terrain, l_two_way_coupling); } // This is simply an advected scalar for convenience start_comp = RhoScalar_comp; num_comp = 1; - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, - horiz_adv_type, vert_adv_type, - l_use_terrain); + AdvectionSrcForScalars(level, finest_level, mfi, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dt, dxInv, mf_m, + horiz_adv_type, vert_adv_type, + l_use_terrain, l_two_way_coupling); #ifdef ERF_USE_MOISTURE start_comp = RhoQt_comp; @@ -283,10 +284,10 @@ void erf_slow_rhs_post (int level, moist_horiz_adv_type = EfficientAdvType(nrk,ac.moistscal_horiz_adv_type); moist_vert_adv_type = EfficientAdvType(nrk,ac.moistscal_vert_adv_type); } - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + AdvectionSrcForScalars(level, finest_level, mfi, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dt, dxInv, mf_m, moist_horiz_adv_type, moist_vert_adv_type, - l_use_terrain); + l_use_terrain, l_two_way_coupling); #elif defined(ERF_USE_WARM_NO_PRECIP) start_comp = RhoQv_comp; @@ -300,10 +301,10 @@ void erf_slow_rhs_post (int level, moist_vert_adv_type = EfficientAdvType(nrk,solverChoice.moistscal_vert_adv_type); } - AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + AdvectionSrcForScalars(level, finest_level, mfi, tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dt, dxInv, mf_m, moist_horiz_adv_type, moist_vert_adv_type, - l_use_terrain); + l_use_terrain, l_two_way_coupling); #endif if (l_use_diff) { diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 58a768be7..50b509895 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -119,7 +119,7 @@ void erf_slow_rhs_pre (int level, const AdvType l_horiz_adv_type = solverChoice.advChoice.dycore_horiz_adv_type; const AdvType l_vert_adv_type = solverChoice.advChoice.dycore_vert_adv_type; const bool l_use_terrain = solverChoice.use_terrain; - const bool l_moving_terrain = (solverChoice.terrain_type == 1); + const bool l_moving_terrain = (solverChoice.terrain_type == TerrainType::Moving); if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain); const bool l_use_ndiff = solverChoice.use_NumDiff; @@ -801,8 +801,7 @@ void erf_slow_rhs_pre (int level, q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i-1,j,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); #endif - rho_u_rhs(i, j, k) += -gpx / (1.0 + q) - - abl_pressure_grad[0] + rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i-1,j,k,Rho_comp)) * abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) @@ -810,8 +809,7 @@ void erf_slow_rhs_pre (int level, { Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k)); Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k)); - rho_u_rhs(i, j, k) += coriolis_factor * - (rho_v_loc * sinphi - rho_w_loc * cosphi); + rho_u_rhs(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi); } // Add Rayleigh damping @@ -846,8 +844,7 @@ void erf_slow_rhs_pre (int level, q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i-1,j,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); #endif - rho_u_rhs(i, j, k) += -gpx / (1.0 + q) - - abl_pressure_grad[0] + rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i-1,j,k,Rho_comp)) * abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) @@ -855,8 +852,7 @@ void erf_slow_rhs_pre (int level, { Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k)); Real rho_w_loc = 0.25 * (rho_w(i,j,k+1) + rho_w(i,j,k) + rho_w(i,j-1,k+1) + rho_w(i,j-1,k)); - rho_u_rhs(i, j, k) += coriolis_factor * - (rho_v_loc * sinphi - rho_w_loc * cosphi); + rho_u_rhs(i, j, k) += coriolis_factor * (rho_v_loc * sinphi - rho_w_loc * cosphi); } // Add Rayleigh damping @@ -911,8 +907,7 @@ void erf_slow_rhs_pre (int level, q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j-1,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); #endif - rho_v_rhs(i, j, k) += -gpy / (1.0_rt + q) - - abl_pressure_grad[1] + rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j-1,k,Rho_comp)) * abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -953,8 +948,7 @@ void erf_slow_rhs_pre (int level, q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j-1,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); #endif - rho_v_rhs(i, j, k) += -gpy / (1.0_rt + q) - - abl_pressure_grad[1] + rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j-1,k,Rho_comp)) * abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) @@ -1007,8 +1001,7 @@ void erf_slow_rhs_pre (int level, q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j,k-1,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); #endif - rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz) / (1.0_rt + q) - - abl_pressure_grad[2] + rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j,k-1,Rho_comp)) * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) @@ -1048,8 +1041,7 @@ void erf_slow_rhs_pre (int level, q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j,k-1,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); #endif - rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz) / (1.0_rt + q) - - abl_pressure_grad[2] + rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j,k-1,Rho_comp)) * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) diff --git a/Source/TimeIntegration/TI_fast_rhs_fun.H b/Source/TimeIntegration/TI_fast_rhs_fun.H index 6ee23f68c..f27ea13ad 100644 --- a/Source/TimeIntegration/TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/TI_fast_rhs_fun.H @@ -24,7 +24,7 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, // Moving terrain std::unique_ptr z_t_pert; - if ( solverChoice.use_terrain && (solverChoice.terrain_type == 1) ) + if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) { // Make "old" fast geom -- store in z_phys_nd for convenience if (verbose) Print() << "Making geometry at start of substep time: " << old_substep_time << std::endl; @@ -94,7 +94,7 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub, dtau, beta_s, inv_fac, mapfac_m[level], mapfac_u[level], mapfac_v[level]); } - } else if (solverChoice.use_terrain && solverChoice.terrain_type == 0) { + } else if (solverChoice.use_terrain && solverChoice.terrain_type == TerrainType::Static) { if (fast_step == 0) { // If this is the first substep we make the coefficients since they are based only on stage data diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H index f1eb5fbaf..4c226be3b 100644 --- a/Source/TimeIntegration/TI_headers.H +++ b/Source/TimeIntegration/TI_headers.H @@ -62,7 +62,7 @@ void erf_slow_rhs_pre (int level, int nrk, * Function for computing the slow RHS for the evolution equations for the scalars other than density or potential temperature * */ -void erf_slow_rhs_post (int level, int nrk, +void erf_slow_rhs_post (int level, int finest_level, int nrk, amrex::Real dt, amrex::Vector& S_rhs, amrex::Vector& S_old, diff --git a/Source/TimeIntegration/TI_no_substep_fun.H b/Source/TimeIntegration/TI_no_substep_fun.H index 3e0f57383..14bdc3d37 100644 --- a/Source/TimeIntegration/TI_no_substep_fun.H +++ b/Source/TimeIntegration/TI_no_substep_fun.H @@ -49,7 +49,7 @@ Array4* fslow = fslow_d.dataPtr(); // Moving terrain - if ( solverChoice.use_terrain && solverChoice.terrain_type == 1 ) + if ( solverChoice.use_terrain && solverChoice.terrain_type == TerrainType::Moving ) { const Array4& dJ_old = detJ_cc[level]->const_array(mfi); const Array4& dJ_new = detJ_cc_new[level]->const_array(mfi); diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index f30077ddc..fd599247b 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -18,7 +18,7 @@ Real slow_dt = new_stage_time - old_step_time; // Moving terrain - if ( solverChoice.use_terrain && (solverChoice.terrain_type == 1) ) + if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) { // Note that the "old" and "new" metric terms correspond to // t^n and the RK stage (either t^*, t^** or t^{n+1} that this source @@ -203,7 +203,7 @@ #endif // Compute RHS for fine interior ghost - if (level>0 && coupling_type=="OneWay" && cf_width>0) { + if (level>0 && solverChoice.coupling_type == CouplingType::OneWay && cf_width>0) { fine_compute_interior_ghost_rhs(new_stage_time, slow_dt, cf_width, cf_set_width, fine_geom, &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1], @@ -267,8 +267,8 @@ #endif // Moving terrain - if ( solverChoice.use_terrain && (solverChoice.terrain_type == 1) ) { - erf_slow_rhs_post(level, nrk, slow_dt, + if ( solverChoice.use_terrain && (solverChoice.terrain_type == TerrainType::Moving) ) { + erf_slow_rhs_post(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_new, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, source, SmnSmn, eddyDiffs, @@ -283,7 +283,7 @@ #endif ); } else { - erf_slow_rhs_post(level, nrk, slow_dt, + erf_slow_rhs_post(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_new, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, source, SmnSmn, eddyDiffs, diff --git a/Source/Utils/Interpolation_UPW.H b/Source/Utils/Interpolation_UPW.H index 89586d5e2..e5f061dea 100644 --- a/Source/Utils/Interpolation_UPW.H +++ b/Source/Utils/Interpolation_UPW.H @@ -18,19 +18,13 @@ struct CENTERED2 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real /*upw_hi*/, amrex::Real /*upw_lo*/) const { // Data to interpolate on - amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i-1, j , k , qty_index); - // Interpolate hi - val_hi = Evaluate(sp1,s); - // Interpolate lo val_lo = Evaluate(s,sm1); } @@ -42,19 +36,13 @@ struct CENTERED2 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real /*upw_hi*/, amrex::Real /*upw_lo*/) const { // Data to interpolate on - amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i , j-1, k , qty_index); - // Interpolate hi - val_hi = Evaluate(sp1,s); - // Interpolate lo val_lo = Evaluate(s,sm1); } @@ -62,7 +50,7 @@ struct CENTERED2 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -77,24 +65,6 @@ struct CENTERED2 val_lo = Evaluate(s,sm1); } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real /*upw_hi*/) const - { - // Data to interpolate on - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - - // Interpolate hi - val_hi = Evaluate(sp1,s); - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -128,25 +98,18 @@ struct UPWIND3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i-1, j , k , qty_index); amrex::Real sm2 = m_phi(i-2, j , k , qty_index); // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; if (upw_lo != 0.) upw_lo = (upw_lo > 0) ? 1. : -1.; - // Interpolate hi - val_hi = Evaluate(sp2,sp1,s,sm1,upw_hi); - // Interpolate lo val_lo = Evaluate(sp1,s,sm1,sm2,upw_lo); } @@ -158,25 +121,18 @@ struct UPWIND3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i , j-1, k , qty_index); amrex::Real sm2 = m_phi(i , j-2, k , qty_index); // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; if (upw_lo != 0.) upw_lo = (upw_lo > 0) ? 1. : -1.; - // Interpolate hi - val_hi = Evaluate(sp2,sp1,s,sm1,upw_hi); - // Interpolate lo val_lo = Evaluate(sp1,s,sm1,sm2,upw_lo); } @@ -184,7 +140,7 @@ struct UPWIND3 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -204,29 +160,6 @@ struct UPWIND3 val_lo = Evaluate(sp1,s,sm1,sm2,upw_lo); } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - - // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; - - // Interpolate lo - val_hi = Evaluate(sp2,sp1,s,sm1,upw_hi); - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -268,21 +201,15 @@ struct CENTERED4 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real /*upw_hi*/, amrex::Real /*upw_lo*/) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i-1, j , k , qty_index); amrex::Real sm2 = m_phi(i-2, j , k , qty_index); - // Interpolate hi - val_hi = Evaluate(sp2,sp1,s,sm1); - // Interpolate lo val_lo = Evaluate(sp1,s,sm1,sm2); } @@ -294,21 +221,15 @@ struct CENTERED4 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real /*upw_hi*/, amrex::Real /*upw_lo*/) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i , j-1, k , qty_index); amrex::Real sm2 = m_phi(i , j-2, k , qty_index); - // Interpolate hi - val_hi = Evaluate(sp2,sp1,s,sm1); - // Interpolate lo val_lo = Evaluate(sp1,s,sm1,sm2); } @@ -316,7 +237,7 @@ struct CENTERED4 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -333,26 +254,6 @@ struct CENTERED4 val_lo = Evaluate(sp1,s,sm1,sm2); } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real /*upw_hi*/) const - { - // Data to interpolate on - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - - // Interpolate hi - val_hi = Evaluate(sp2,sp1,s,sm1); - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -390,13 +291,10 @@ struct UPWIND5 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i+3, j , k , qty_index); amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -405,12 +303,8 @@ struct UPWIND5 amrex::Real sm3 = m_phi(i-3, j , k , qty_index); // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; if (upw_lo != 0.) upw_lo = (upw_lo > 0) ? 1. : -1.; - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2,upw_hi); - // Interpolate lo val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3,upw_lo); } @@ -422,13 +316,10 @@ struct UPWIND5 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i , j+3, k , qty_index); amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -437,12 +328,8 @@ struct UPWIND5 amrex::Real sm3 = m_phi(i , j-3, k , qty_index); // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; if (upw_lo != 0.) upw_lo = (upw_lo > 0) ? 1. : -1.; - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2,upw_hi); - // Interpolate lo val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3,upw_lo); } @@ -450,7 +337,7 @@ struct UPWIND5 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -472,31 +359,6 @@ struct UPWIND5 val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3,upw_lo); } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp3 = m_phi(i , j , k+3, qty_index); - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - amrex::Real sm2 = m_phi(i , j , k-2, qty_index); - - // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; - - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2,upw_hi); - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -542,13 +404,10 @@ struct CENTERED6 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real /*upw_hi*/, amrex::Real /*upw_lo*/) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i+3, j , k , qty_index); amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -556,9 +415,6 @@ struct CENTERED6 amrex::Real sm2 = m_phi(i-2, j , k , qty_index); amrex::Real sm3 = m_phi(i-3, j , k , qty_index); - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2); - // Interpolate lo val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3); } @@ -570,13 +426,10 @@ struct CENTERED6 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real /*upw_hi*/, amrex::Real /*upw_lo*/) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i , j+3, k , qty_index); amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -584,9 +437,6 @@ struct CENTERED6 amrex::Real sm2 = m_phi(i , j-2, k , qty_index); amrex::Real sm3 = m_phi(i , j-3, k , qty_index); - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2); - // Interpolate lo val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3); } @@ -594,7 +444,7 @@ struct CENTERED6 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -613,28 +463,6 @@ struct CENTERED6 val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3); } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real /*upw_hi*/) const - { - // Data to interpolate on - amrex::Real sp3 = m_phi(i , j , k+3, qty_index); - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - amrex::Real sm2 = m_phi(i , j , k-2, qty_index); - - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2); - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -672,7 +500,7 @@ struct UPWINDALL AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -700,37 +528,6 @@ struct UPWINDALL } } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi, - const AdvType adv_type) const - { - if (adv_type == AdvType::Centered_2nd) { - val_hi = 0.5 * ( m_phi(i,j,k,qty_index) + m_phi(i,j,k+1,qty_index) ); - return; - } else { - // Data to interpolate on - amrex::Real sp3 = (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th ) ? m_phi(i , j , k+3, qty_index): 0.; - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - amrex::Real sm2 = (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th ) ? m_phi(i , j , k-2, qty_index) : 0.; - - // Upwinding flags - if (upw_hi != 0.) upw_hi = (upw_hi > 0) ? 1. : -1.; - - // Interpolate hi - val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2,upw_hi,adv_type); - } - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real diff --git a/Source/Utils/Interpolation_WENO.H b/Source/Utils/Interpolation_WENO.H index bf9a24f40..01a4ca68d 100644 --- a/Source/Utils/Interpolation_WENO.H +++ b/Source/Utils/Interpolation_WENO.H @@ -18,26 +18,15 @@ struct WENO3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i-1, j , k , qty_index); amrex::Real sm2 = m_phi(i-2, j , k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm2,sm1,s ); } else if (upw_lo < -tol) { @@ -54,26 +43,15 @@ struct WENO3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i , j-1, k , qty_index); amrex::Real sm2 = m_phi(i , j-2, k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm2,sm1,s ); } else if (upw_lo < -tol) { @@ -86,7 +64,7 @@ struct WENO3 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -108,31 +86,6 @@ struct WENO3 } } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -182,13 +135,10 @@ struct WENO5 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i+3, j , k , qty_index); amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -196,14 +146,6 @@ struct WENO5 amrex::Real sm2 = m_phi(i-2, j , k , qty_index); amrex::Real sm3 = m_phi(i-3, j , k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm2,sm1,s,sp1,sp2); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp3,sp2,sp1,s ,sm1); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm3,sm2,sm1,s ,sp1); } else if (upw_lo < -tol) { @@ -220,13 +162,10 @@ struct WENO5 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i , j+3, k , qty_index); amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -234,14 +173,6 @@ struct WENO5 amrex::Real sm2 = m_phi(i , j-2, k , qty_index); amrex::Real sm3 = m_phi(i , j-3, k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm2,sm1,s,sp1,sp2); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp3,sp2,sp1,s ,sm1); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm3,sm2,sm1,s ,sp1); } else if (upw_lo < -tol) { @@ -254,7 +185,7 @@ struct WENO5 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -278,33 +209,6 @@ struct WENO5 } } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp3 = m_phi(i , j , k+3, qty_index); - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - amrex::Real sm2 = m_phi(i , j , k-2, qty_index); - - if (upw_hi > tol){ - val_hi = Evaluate(sm2,sm1,s,sp1,sp2); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp3,sp2,sp1,s ,sm1); - } else { - val_hi = 0.5 * (s + sp1); - } - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real diff --git a/Source/Utils/Interpolation_WENO_Z.H b/Source/Utils/Interpolation_WENO_Z.H index 0d77bec86..dc65811b0 100644 --- a/Source/Utils/Interpolation_WENO_Z.H +++ b/Source/Utils/Interpolation_WENO_Z.H @@ -18,26 +18,15 @@ struct WENO_Z3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i-1, j , k , qty_index); amrex::Real sm2 = m_phi(i-2, j , k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm2,sm1,s ); } else if (upw_lo < -tol) { @@ -54,26 +43,15 @@ struct WENO_Z3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i , j-1, k , qty_index); amrex::Real sm2 = m_phi(i , j-2, k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm2,sm1,s ); } else if (upw_lo < -tol) { @@ -86,7 +64,7 @@ struct WENO_Z3 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -108,31 +86,6 @@ struct WENO_Z3 } } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -183,26 +136,15 @@ struct WENO_MZQ3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i-1, j , k , qty_index); amrex::Real sm2 = m_phi(i-2, j , k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm2,sm1,s ); } else if (upw_lo < -tol) { @@ -219,26 +161,15 @@ struct WENO_MZQ3 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); amrex::Real sm1 = m_phi(i , j-1, k , qty_index); amrex::Real sm2 = m_phi(i , j-2, k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm2,sm1,s ); } else if (upw_lo < -tol) { @@ -251,7 +182,7 @@ struct WENO_MZQ3 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -273,31 +204,6 @@ struct WENO_MZQ3 } } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - - if (upw_hi > tol){ - val_hi = Evaluate(sm1,s ,sp1); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp2,sp1,s ); - } else { - val_hi = 0.5 * (s + sp1); - } - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -353,13 +259,10 @@ struct WENO_Z5 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i+3, j , k , qty_index); amrex::Real sp2 = m_phi(i+2, j , k , qty_index); amrex::Real sp1 = m_phi(i+1, j , k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -367,14 +270,6 @@ struct WENO_Z5 amrex::Real sm2 = m_phi(i-2, j , k , qty_index); amrex::Real sm3 = m_phi(i-3, j , k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm2,sm1,s,sp1,sp2); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp3,sp2,sp1,s ,sm1); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm3,sm2,sm1,s ,sp1); } else if (upw_lo < -tol) { @@ -391,13 +286,10 @@ struct WENO_Z5 const int& j, const int& k, const int& qty_index, - amrex::Real& val_hi, amrex::Real& val_lo, - amrex::Real upw_hi, amrex::Real upw_lo) const { // Data to interpolate on - amrex::Real sp3 = m_phi(i , j+3, k , qty_index); amrex::Real sp2 = m_phi(i , j+2, k , qty_index); amrex::Real sp1 = m_phi(i , j+1, k , qty_index); amrex::Real s = m_phi(i , j , k , qty_index); @@ -405,14 +297,6 @@ struct WENO_Z5 amrex::Real sm2 = m_phi(i , j-2, k , qty_index); amrex::Real sm3 = m_phi(i , j-3, k , qty_index); - if (upw_hi > tol){ - val_hi = Evaluate(sm2,sm1,s,sp1,sp2); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp3,sp2,sp1,s ,sm1); - } else { - val_hi = 0.5 * (s + sp1); - } - if (upw_lo > tol) { val_lo = Evaluate(sm3,sm2,sm1,s ,sp1); } else if (upw_lo < -tol) { @@ -425,7 +309,7 @@ struct WENO_Z5 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void - InterpolateInZ_lo (const int& i, + InterpolateInZ (const int& i, const int& j, const int& k, const int& qty_index, @@ -449,33 +333,6 @@ struct WENO_Z5 } } - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE - void - InterpolateInZ_hi (const int& i, - const int& j, - const int& k, - const int& qty_index, - amrex::Real& val_hi, - amrex::Real upw_hi) const - { - // Data to interpolate on - amrex::Real sp3 = m_phi(i , j , k+3, qty_index); - amrex::Real sp2 = m_phi(i , j , k+2, qty_index); - amrex::Real sp1 = m_phi(i , j , k+1, qty_index); - amrex::Real s = m_phi(i , j , k , qty_index); - amrex::Real sm1 = m_phi(i , j , k-1, qty_index); - amrex::Real sm2 = m_phi(i , j , k-2, qty_index); - - if (upw_hi > tol){ - val_hi = Evaluate(sm2,sm1,s,sp1,sp2); - } else if (upw_hi < -tol) { - val_hi = Evaluate(sp3,sp2,sp1,s ,sm1); - } else { - val_hi = 0.5 * (s + sp1); - } - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real diff --git a/Source/Utils/TerrainMetrics.cpp b/Source/Utils/TerrainMetrics.cpp index 4182a92b5..175ccd2c3 100644 --- a/Source/Utils/TerrainMetrics.cpp +++ b/Source/Utils/TerrainMetrics.cpp @@ -69,7 +69,6 @@ init_zlevels (amrex::Vector& zlevels_stag, void init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd, amrex::Vector const& z_levels_h) { - auto dx = geom.CellSizeArray(); auto ProbHiArr = geom.ProbHiArray(); // z_nd is nodal in all directions @@ -84,9 +83,6 @@ init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd, amrex::Vector z_levels_d; z_levels_d.resize(nz); #ifdef AMREX_USE_GPU diff --git a/Tests/test_files/DensityCurrent_detJ2_MT/DensityCurrent_detJ2_MT.i b/Tests/test_files/DensityCurrent_detJ2_MT/DensityCurrent_detJ2_MT.i index 1382e7ba3..0498ffd33 100644 --- a/Tests/test_files/DensityCurrent_detJ2_MT/DensityCurrent_detJ2_MT.i +++ b/Tests/test_files/DensityCurrent_detJ2_MT/DensityCurrent_detJ2_MT.i @@ -67,6 +67,6 @@ erf.init_shrink = 1.0 # scale back initial timestep erf.terrain_z_levels = 0. 100. 200. 300. 400. 500. 600. 700. 800. 900. 1000. 1100. 1200. 1300. 1400. 1500. 1600. 1700. 1800. 1900. 2000. 2100. 2200. 2300. 2400. 2500. 2600. 2700. 2800. 2900. 3000. 3100. 3200. 3300. 3400. 3500. 3600. 3700. 3800. 3900. 4000. 4100. 4200. 4300. 4400. 4500. 4600. 4700. 4800. 4900. 5000. 5100. 5200. 5300. 5400. 5500. 5600. 5700. 5800. 5900. 6000. 6100. 6200. 6300. 6400. # TERRRAIN GRID TYPE -erf.use_terrain = 1 -erf.terrain_type = 1 +erf.use_terrain = true +erf.terrain_type = Moving erf.terrain_smoothing = 2 // Sullivan 2004 approach diff --git a/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i b/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i index 0642c91fb..e100c34f4 100644 --- a/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i +++ b/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i @@ -15,8 +15,8 @@ zlo.type = "SlipWall" zhi.type = "SlipWall" # TERRRAIN GRID TYPE -erf.use_terrain = 1 # enable terrain stencils -erf.terrain_type = 1 # moving terrain +erf.use_terrain = true # enable terrain stencils +erf.terrain_type = Moving # moving terrain erf.terrain_smoothing = 2 # Sullivan 2004 approach erf.no_substepping = 1 diff --git a/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i b/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i index bb9e799c2..9186f4556 100644 --- a/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i +++ b/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i @@ -15,8 +15,8 @@ zlo.type = "SlipWall" zhi.type = "SlipWall" # TERRRAIN GRID TYPE -erf.use_terrain = 1 # enable terrain stencils -erf.terrain_type = 1 # moving terrain +erf.use_terrain = true # enable terrain stencils +erf.terrain_type = Moving # moving terrain erf.terrain_smoothing = 2 # Sullivan 2004 approach erf.use_lagged_delta_rt = false From 80397f62e529e52d1171563ed0554e40f97d525a Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 7 Nov 2023 09:29:34 -0800 Subject: [PATCH 2/2] Fix for extrapolation of k_turb (#1285) Co-authored-by: Aaron Lattanzi --- Source/Diffusion/ComputeTurbulentViscosity.cpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index bf3709cf4..2cc5d5946 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -61,7 +61,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu #endif for (amrex::MFIter mfi(eddyViscosity,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - // NOTE: This gets us the lateral ghost cells for lev>1; which + // NOTE: This gets us the lateral ghost cells for lev>0; which // have been filled from FP Two Levels. Box bxcc = mfi.growntilebox() & domain; @@ -164,7 +164,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu } } - // Extrapolate Kturb in extrap x/y, fill remaining elements + // Extrapolate Kturb in x/y, fill remaining elements (relevent to lev==0) //*********************************************************************************** int ngc(1); Real inv_Pr_t = turbChoice.Pr_t_inv; @@ -194,7 +194,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu { Box bxcc = mfi.tilebox(); Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1); - Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); + Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); planey.grow(0,1); int i_lo = bxcc.smallEnd(0); int i_hi = bxcc.bigEnd(0); int j_lo = bxcc.smallEnd(1); int j_hi = bxcc.bigEnd(1); bxcc.growLo(0,ngc); bxcc.growHi(0,ngc); @@ -222,15 +222,17 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu if (j_lo == domain.smallEnd(1)) { amrex::ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - mu_turb(i, j_lo-j, k, EddyDiff::Mom_h) = mu_turb(i, j_lo, k, EddyDiff::Mom_h); - mu_turb(i, j_lo-j, k, EddyDiff::Mom_v) = mu_turb(i, j_lo, k, EddyDiff::Mom_v); + int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0)); + mu_turb(i, j_lo-j, k, EddyDiff::Mom_h) = mu_turb(li, j_lo, k, EddyDiff::Mom_h); + mu_turb(i, j_lo-j, k, EddyDiff::Mom_v) = mu_turb(li, j_lo, k, EddyDiff::Mom_v); }); } if (j_hi == domain.bigEnd(1)) { amrex::ParallelFor(planey, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - mu_turb(i, j_hi+j, k, EddyDiff::Mom_h) = mu_turb(i, j_hi, k, EddyDiff::Mom_h); - mu_turb(i, j_hi+j, k, EddyDiff::Mom_v) = mu_turb(i, j_hi, k, EddyDiff::Mom_v); + int li = amrex::min(amrex::max(i, domain.smallEnd(0)), domain.bigEnd(0)); + mu_turb(i, j_hi+j, k, EddyDiff::Mom_h) = mu_turb(li, j_hi, k, EddyDiff::Mom_h); + mu_turb(i, j_hi+j, k, EddyDiff::Mom_v) = mu_turb(li, j_hi, k, EddyDiff::Mom_v); }); }