From 5f1f91626011d4ee56ad17c3d685eb1cd02cd066 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 24 Jul 2023 11:02:45 -0700 Subject: [PATCH] Update advection flags (#1180) * Update advection_type flags * fix trailing white space * update advection stuff * Remove WENO name conflict from interpolation struct and index defines. --------- Co-authored-by: AMLattanzi --- Docs/sphinx_doc/Discretizations.rst | 9 +- Docs/sphinx_doc/ERFvsWRF.rst | 3 +- Docs/sphinx_doc/Inputs.rst | 130 +++++----- Docs/sphinx_doc/RegressionTests.rst | 10 +- Exec/ABL/inputs_most_msf | 6 +- Exec/ABL/inputs_most_no_msf | 6 +- Exec/ABL_input_sounding/abl_mwe.i | 7 +- Exec/DevTests/ParticlesOverWoA/inputs_zlevels | 1 - Exec/RegTests/Bubble/inputs_squall2d_x | 1 - Exec/RegTests/ScalarAdvDiff/inputs_WENO | 13 +- Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z | 10 +- Exec/RegTests/WitchOfAgnesi/inputs_zlevels | 1 - Source/Advection/Advection.H | 14 +- Source/Advection/AdvectionSrcForMom.cpp | 70 +++--- Source/Advection/AdvectionSrcForMom_N.H | 108 +++++---- Source/Advection/AdvectionSrcForMom_T.H | 104 ++++---- Source/Advection/AdvectionSrcForState.cpp | 227 ++++++++---------- Source/Advection/AdvectionSrcForState_N.H | 86 ++++--- Source/Advection/AdvectionSrcForState_T.H | 42 ++-- Source/DataStruct.H | 213 +++++++++++++--- Source/ERF.H | 62 +++-- Source/ERF_make_new_level.cpp | 20 +- Source/IO/NCCheckpoint.cpp | 9 +- Source/IndexDefines.H | 16 +- Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 8 +- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 47 ++-- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 8 +- Source/TimeIntegration/TI_slow_rhs_fun.H | 6 +- Source/Utils/Interpolation.H | 79 +++--- Source/Utils/Interpolation_UPW.H | 58 ++--- Source/Utils/Interpolation_WENO_Z.H | 165 +++++++++++++ .../Bubble_DensityCurrent.i | 2 - Tests/test_files/CouetteFlow/CouetteFlow.i | 3 - .../Deardorff_stationary.i | 1 - .../DensityCurrent/DensityCurrent.i | 2 - .../DensityCurrent_detJ2.i | 2 - .../DensityCurrent_detJ2_MT.i | 2 - .../DensityCurrent_detJ2_nosub.i | 3 - Tests/test_files/EkmanSpiral/EkmanSpiral.i | 3 - .../IsentropicVortexAdvecting.i | 3 - .../IsentropicVortexStationary.i | 3 - .../MSF_NoSub_IsentropicVortexAdv.i | 3 - .../MSF_Sub_IsentropicVortexAdv.i | 3 - .../MovingTerrain_nosub/MovingTerrain_nosub.i | 3 - .../MovingTerrain_sub/MovingTerrain_sub.i | 3 - .../PoiseuilleFlow/PoiseuilleFlow.i | 3 - .../RayleighDamping/RayleighDamping.i | 3 - .../ScalarAdvDiff_order2.i | 3 - .../ScalarAdvDiff_order3.i | 6 +- .../ScalarAdvDiff_order4.i | 6 +- .../ScalarAdvDiff_order5.i | 6 +- .../ScalarAdvDiff_order6.i | 6 +- .../ScalarAdvectionShearedU.i | 3 - .../ScalarAdvectionUniformU.i | 8 +- .../ScalarDiffusionGaussian.i | 3 - .../ScalarDiffusionSine/ScalarDiffusionSine.i | 3 - .../TaylorGreenAdvecting.i | 3 - .../TaylorGreenAdvectingDiffusing.i | 3 - 58 files changed, 959 insertions(+), 672 deletions(-) diff --git a/Docs/sphinx_doc/Discretizations.rst b/Docs/sphinx_doc/Discretizations.rst index 2f42782af..09abf73f9 100644 --- a/Docs/sphinx_doc/Discretizations.rst +++ b/Docs/sphinx_doc/Discretizations.rst @@ -230,12 +230,15 @@ With the WENO5 scheme, one has :math:`N=3, \; \omega_{1} = 1/10, \; \omega_{2} = q_{m + \frac{1}{2}}^{(3)} = \frac{1}{3} q_{m} + \frac{5}{6} q_{m+1} - \frac{1}{6} q_{m+2} \end{array} -By default, the WENO3 scheme will be employed for moisture variables if the code is compiled with moisture support. However, users may utilize the WENO scheme for all scalar variables, with a specified order, via the following inputs: +By default, the WENO3 scheme will be employed for moisture variables if the code is compiled with moisture support. +However, users may utilize the WENO scheme for dry scalar variables as well. The scheme for each type is specified by :: - erf.all_use_WENO = true # Default is false - erf.spatial_order_WENO = 5 # Default is 3 + erf.dryscal_horiz_adv_type = + erf.dryscal_vert_adv_type = + erf.moistscal_horiz_adv_type = + erf.moistscal_vert_adv_type = Ref: Muñoz-Esparza, D., Sauer, J. A., Jensen, A. A., Xue, L., & Grabowski, W. W. (2022). The FastEddy® resident-GPU accelerated large-eddy simulation framework: Moist dynamics extension, validation and sensitivities of modeling non-precipitating shallow cumulus clouds. Journal of Advances in Modeling Earth Systems, 14, e2021MS002904. `https://onlinelibrary.wiley.com/doi/10.1029/2021MS002904>`_ diff --git a/Docs/sphinx_doc/ERFvsWRF.rst b/Docs/sphinx_doc/ERFvsWRF.rst index 143180860..90d60c704 100644 --- a/Docs/sphinx_doc/ERFvsWRF.rst +++ b/Docs/sphinx_doc/ERFvsWRF.rst @@ -24,7 +24,8 @@ cloud water / ice mixing ratio. **Time Integration**: Time-split integration using 3rd-order Runge-Kutta scheme with smaller time step for acoustic and gravity wave modes. Variable time step capability. -**Spatial Discretization**: 2nd- to 6th-order advection options in horizontal and vertical +**Spatial Discretization**: 2nd- to 6th-order advection options in horizontal and vertical. In addition, several +different WENO schemes are available for scalar variables other than density and potential temperature. **Turbulent Mixing**: Sub-grid scale turbulence formulation. Vertically implicit acoustic step off-centering. diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index dc426734f..3f6f26780 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -491,60 +491,82 @@ Diffusive Physics List of Parameters ------------------ -+----------------------------------+--------------------+---------------------+-------------+ -| Parameter | Definition | Acceptable | Default | -| | | Values | | -+==================================+====================+=====================+=============+ -| **erf.alpha_T** | Diffusion coeff. | Real | 0.0 | -| | for temperature | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.alpha_C** | Diffusion coeff. | Real | 0.0 | -| | for scalar | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.rho0_trans** | Reference density | Real | 1.0 | -| | to compute const. | | | -| | rho*Alpha | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.les_type** | Using an LES | "None", | "None" | -| | model, and if so, | "Smagorinsky", | | -| | which type? | "Deardorff" | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.molec_diff_type** | Using molecular | "None", | "None" | -| | viscosity and | "Constant", or | | -| | diffusivity? | "ConstantAlpha" | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.dynamicViscosity** | Viscous coeff. if | Real | 0.0 | -| | DNS | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.Cs** | Constant | Real | 0.0 | -| | Smagorinsky coeff. | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.Pr_t** | Turbulent Prandtl | Real | 1.0 | -| | Number | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.Sc_t** | Turbulent Schmidt | Real | 1.0 | -| | Number | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.horiz/vert_spatial_order** | Horizontal and | 2 / 3 / 4 / 5 / 6 | 2 | -| | vertical spatial | | | -| | order | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.all_use_WENO** | Use WENO for all | "True", | "False" | -| | scalar variables | "False" | | -| | | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.spatial_order_WENO** | Spatial order | 3/5 | 3 | -| | for WENO | | | -| | | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.use_NumDiff** | Use 6th order | "True", | "False" | -| | numerical diffusion| "False" | | -| | | | | -+----------------------------------+--------------------+---------------------+-------------+ -| **erf.NumDiffCoeff** | Coefficient for | Real | 0.0 | -| | 6th order | [0.0, 1.0] | | -| | numerical diffusion| | | -+----------------------------------+--------------------+---------------------+-------------+ ++----------------------------------+--------------------+---------------------+--------------+ +| Parameter | Definition | Acceptable | Default | +| | | Values | | ++==================================+====================+=====================+==============+ +| **erf.alpha_T** | Diffusion coeff. | Real | 0.0 | +| | for temperature | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.alpha_C** | Diffusion coeff. | Real | 0.0 | +| | for scalar | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.rho0_trans** | Reference density | Real | 1.0 | +| | to compute const. | | | +| | rho*Alpha | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.les_type** | Using an LES | "None", | "None" | +| | model, and if so, | "Smagorinsky", | | +| | which type? | "Deardorff" | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.molec_diff_type** | Using molecular | "None", | "None" | +| | viscosity and | "Constant", or | | +| | diffusivity? | "ConstantAlpha" | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.dynamicViscosity** | Viscous coeff. if | Real | 0.0 | +| | DNS | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.Cs** | Constant | Real | 0.0 | +| | Smagorinsky coeff. | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.Pr_t** | Turbulent Prandtl | Real | 1.0 | +| | Number | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.Sc_t** | Turbulent Schmidt | Real | 1.0 | +| | Number | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.dycore_horiz_adv_type** | Horizontal | see below | Centered_2nd | +| | advection type | | | +| | for dycore vars | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.dycore_vert_adv_type** | Vertical | see below | Centered_2nd | +| | advection type | | | +| | for dycore vars | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.dryscal_horiz_adv_type** | Horizontal | see below | Centered_2nd | +| | advection type | | | +| | for dry scalars | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.dryscal_vert_adv_type** | Vertical | see below | Centered_2nd | +| | advection type | | | +| | for dry scalars | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.moistscal_horiz_adv_type** | Horizontal | see below | WENO3 | +| | advection type | | | +| | for dry scalars | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.moistscal_vert_adv_type** | Vertical | see below | WENO3 | +| | advection type | | | +| | for dry scalars | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.use_NumDiff** | Use 6th order | "True", | "False" | +| | numerical diffusion| "False" | | +| | | | | ++----------------------------------+--------------------+---------------------+--------------+ +| **erf.NumDiffCoeff** | Coefficient for | Real | 0.0 | +| | 6th order | [0.0, 1.0] | | +| | numerical diffusion| | | ++----------------------------------+--------------------+---------------------+--------------+ + +The allowed advection types for the dycore variables are +"Centered_2nd", "Upwind_3rd", "Centered_4th", "Upwind_5th" and "Centered_6th". + +The allowed advection types for the dry and moist scalars are +"Centered_2nd", "Upwind_3rd", "Centered_4th", "Upwind_5th", "Centered_6th" and in addition, +"WENO3", "WENOZ3", "WENOMZQ3", "WENO5", and "WENOZ5." + +Note: if using WENO schemes, the horizontal and vertical advection types must be set to +the same string. Note: in the equations for the evolution of momentum, potential temperature and advected scalars, the diffusion coefficients are written as :math:`\mu`, :math:`\rho \alpha_T` and :math:`\rho \alpha_C`, respectively. diff --git a/Docs/sphinx_doc/RegressionTests.rst b/Docs/sphinx_doc/RegressionTests.rst index 0eb1a42ca..d5c858677 100644 --- a/Docs/sphinx_doc/RegressionTests.rst +++ b/Docs/sphinx_doc/RegressionTests.rst @@ -75,19 +75,19 @@ The following problems are currently tested in the CI: | | | | | SlipWall | | | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | ScalarAdvDiff_order2 | 32 32 32 | Periodic | Periodic | SlipWall | None | advection + diffusion | -| | | | | SlipWall | | spatial_order = 2 | +| | | | | SlipWall | | "Centered_2nd" | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | ScalarAdvDiff_order3 | 32 32 32 | Periodic | Periodic | SlipWall | None | advection + diffusion | -| | | | | SlipWall | | spatial_order = 3 | +| | | | | SlipWall | | "Upwind_3rd" | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | ScalarAdvDiff_order4 | 32 32 32 | Periodic | Periodic | SlipWall | None | advection + diffusion | -| | | | | SlipWall | | spatial_order = 4 | +| | | | | SlipWall | | "Centered_4th" | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | ScalarAdvDiff_order5 | 32 32 32 | Periodic | Periodic | SlipWall | None | advection + diffusion | -| | | | | SlipWall | | spatial_order = 4 | +| | | | | SlipWall | | "Upwind_5th" | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | ScalarAdvDiff_order6 | 32 32 32 | Periodic | Periodic | SlipWall | None | advection + diffusion | -| | | | | SlipWall | | spatial_order = 6 | +| | | | | SlipWall | | "Centered_6th" | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | ScalarDiffusionGaussian | 16 16 16 | Periodic | Periodic | SlipWall | None | | | | | | | SlipWall | | | diff --git a/Exec/ABL/inputs_most_msf b/Exec/ABL/inputs_most_msf index 1558c131a..bba676fa2 100644 --- a/Exec/ABL/inputs_most_msf +++ b/Exec/ABL/inputs_most_msf @@ -55,8 +55,10 @@ erf.molec_diff_type = "Constant" erf.les_type = "Smagorinsky" erf.Cs = 0.1 -erf.horiz_spatial_order = 3 -erf.vert_spatial_order = 3 +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Exec/ABL/inputs_most_no_msf b/Exec/ABL/inputs_most_no_msf index cfb1d8c09..a2bc1965e 100644 --- a/Exec/ABL/inputs_most_no_msf +++ b/Exec/ABL/inputs_most_no_msf @@ -53,8 +53,10 @@ erf.molec_diff_type = "Constant" erf.les_type = "Smagorinsky" erf.Cs = 0.1 -erf.horiz_spatial_order = 3 -erf.vert_spatial_order = 3 +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Exec/ABL_input_sounding/abl_mwe.i b/Exec/ABL_input_sounding/abl_mwe.i index 9c09cd6f7..e7b21d6ec 100644 --- a/Exec/ABL_input_sounding/abl_mwe.i +++ b/Exec/ABL_input_sounding/abl_mwe.i @@ -53,8 +53,11 @@ erf.plot_int_1 = 100 # number of timesteps between plotfiles erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoKE pres_hse dens_hse pert_pres pert_dens # SOLVER CHOICES -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 +erf.dycore_horiz_adv_type = "Centered_2nd" +erf.dycore_vert_adv_type = "Centered_2nd" +erf.dryscal_horiz_adv_type = "Centered_2nd" +erf.dryscal_vert_adv_type = "Centered_2nd" + erf.use_gravity = true erf.use_gravity = false diff --git a/Exec/DevTests/ParticlesOverWoA/inputs_zlevels b/Exec/DevTests/ParticlesOverWoA/inputs_zlevels index e14a6e1c9..2a1758963 100644 --- a/Exec/DevTests/ParticlesOverWoA/inputs_zlevels +++ b/Exec/DevTests/ParticlesOverWoA/inputs_zlevels @@ -48,7 +48,6 @@ erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pr erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.spatial_order = 2 erf.les_type = "None" # MULTILEVEL diff --git a/Exec/RegTests/Bubble/inputs_squall2d_x b/Exec/RegTests/Bubble/inputs_squall2d_x index 656d95a66..3328ad678 100644 --- a/Exec/RegTests/Bubble/inputs_squall2d_x +++ b/Exec/RegTests/Bubble/inputs_squall2d_x @@ -37,7 +37,6 @@ erf.plot_int_1 = 10 # number of timesteps between plotfiles erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pres_hse # SOLVER CHOICES -erf.spatial_order = 2 erf.use_gravity = true erf.use_coriolis = true erf.latitude = 40.0 # hard-coded in dyn_em/module_initialize_ideal.F diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_WENO b/Exec/RegTests/ScalarAdvDiff/inputs_WENO index 7e26ec756..f59391061 100644 --- a/Exec/RegTests/ScalarAdvDiff/inputs_WENO +++ b/Exec/RegTests/ScalarAdvDiff/inputs_WENO @@ -1,4 +1,4 @@ -# ------------------ INPUTS TO MAIN PROGRAM ------------------- +p ------------------ INPUTS TO MAIN PROGRAM ------------------- max_step = 20 amrex.fpe_trap_invalid = 1 @@ -18,9 +18,12 @@ zhi.type = "SlipWall" erf.use_lowM_dt = 1 erf.cfl = 0.5 # cfl number for hyperbolic system -# TEST WENO SCHEMES -erf.all_use_WENO = true -erf.spatial_order_WENO = 5 # MUST BE 3 or 5 +erf.dycore_horiz_adv_type = "Centered_2nd" +erf.dycore_vert_adv_type = "Centered_2nd" +erf.dryscal_horiz_adv_type = "WENO5" +erf.dryscal_vert_adv_type = "WENO5" +erf.moistscal_horiz_adv_type = "WENO5" +erf.moistscal_vert_adv_type = "WENO5" # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -64,4 +67,4 @@ prob.prob_type = 11 #prob.rad_0 = 0.25 #prob.xc_frac = 0.5 #prob.yc_frac = 0.5 -#prob.zc_frac = 0.5 \ No newline at end of file +#prob.zc_frac = 0.5 diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z b/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z index 9017ef3a5..6b72bc8fa 100644 --- a/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z +++ b/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z @@ -19,8 +19,12 @@ erf.use_lowM_dt = 1 erf.cfl = 0.5 # cfl number for hyperbolic system # TEST WENO SCHEMES -erf.all_use_WENO_Z = true -erf.spatial_order_WENO = 3 # MUST BE 3 or 5 +erf.dycore_horiz_adv_type = "Centered_2nd" +erf.dycore_horiz_adv_type = "Centered_2nd" +erf.dryscal_vert_adv_type = "WENO3" +erf.dryscal_vert_adv_type = "WENO3" +erf.moistscal_vert_adv_type = "WENO3" +erf.moistscal_vert_adv_type = "WENO3" # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -64,4 +68,4 @@ prob.prob_type = 11 #prob.rad_0 = 0.25 #prob.xc_frac = 0.5 #prob.yc_frac = 0.5 -#prob.zc_frac = 0.5 \ No newline at end of file +#prob.zc_frac = 0.5 diff --git a/Exec/RegTests/WitchOfAgnesi/inputs_zlevels b/Exec/RegTests/WitchOfAgnesi/inputs_zlevels index e14a6e1c9..2a1758963 100644 --- a/Exec/RegTests/WitchOfAgnesi/inputs_zlevels +++ b/Exec/RegTests/WitchOfAgnesi/inputs_zlevels @@ -48,7 +48,6 @@ erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pr erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.spatial_order = 2 erf.les_type = "None" # MULTILEVEL diff --git a/Source/Advection/Advection.H b/Source/Advection/Advection.H index 76f0884d1..1bb262e94 100644 --- a/Source/Advection/Advection.H +++ b/Source/Advection/Advection.H @@ -27,12 +27,12 @@ void AdvectionSrcForRhoAndTheta (const amrex::Box& bx, const amrex::Box& valid_b const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int horiz_spatial_order, const int vert_spatial_order, + const int horiz_adv_type, const int vert_adv_type, const int use_terrain); /** Compute advection tendency for all scalars other than density and potential temperature */ void AdvectionSrcForScalars (const amrex::Box& bx, - const int &icomp, const int &ncomp, + const int icomp, const int ncomp, const amrex::Array4& avg_xmom, const amrex::Array4& avg_ymom, const amrex::Array4& avg_zmom, @@ -41,12 +41,8 @@ void AdvectionSrcForScalars (const amrex::Box& bx, const amrex::Array4& detJ, const amrex::GpuArray& cellSize, const amrex::Array4& mf_m, - const bool all_use_WENO, - const bool moist_use_WENO, - const bool all_use_WENO_Z, - const bool moist_use_WENO_Z, - const int spatial_order_WENO, - const int horiz_spatial_order, const int vert_spatial_order, const int use_terrain); + const int horiz_adv_type, const int vert_adv_type, + const int use_terrain); /** Compute advection tendencies for all components of momentum */ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz, @@ -61,7 +57,7 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int horiz_spatial_order, const int vert_spatial_order, + const int horiz_adv_type, const int vert_adv_type, const int use_terrain, const int domhi_z); #endif diff --git a/Source/Advection/AdvectionSrcForMom.cpp b/Source/Advection/AdvectionSrcForMom.cpp index abad50f45..35eb9b1f7 100644 --- a/Source/Advection/AdvectionSrcForMom.cpp +++ b/Source/Advection/AdvectionSrcForMom.cpp @@ -27,8 +27,8 @@ using namespace amrex; * @param[in] mf_m map factor at cell centers * @param[in] mf_u map factor at x-faces * @param[in] mf_v map factor at y-faces - * @param[in] horiz_spatial_order sets the spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order sets the spatial order to be used for vertical derivatives + * @param[in] horiz_adv_type sets the spatial order to be used for lateral derivatives + * @param[in] vert_adv_type sets the spatial order to be used for vertical derivatives * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms) * @param[in] domhi_z maximum k value in the domain */ @@ -49,8 +49,8 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, const Array4& mf_m, const Array4& mf_u, const Array4& mf_v, - const int horiz_spatial_order, - const int vert_spatial_order, + const int horiz_adv_type, + const int vert_adv_type, const int use_terrain, const int domhi_z) { @@ -62,7 +62,8 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, if (!use_terrain) { // Inline with 2nd order for efficiency - if (std::max(horiz_spatial_order,vert_spatial_order) == 2) { + if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) + { ParallelFor(bxx, bxy, bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -135,43 +136,46 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, }); // Template higher order methods } else { - if (horiz_spatial_order == 2) { - AdvectionSrcForMomVert_N(bxx, bxy, bxz, + if (horiz_adv_type == AdvType::Centered_2nd) { + AdvectionSrcForMomVert_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 3) { + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Upwind_3rd) { AdvectionSrcForMomVert_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 4) { - AdvectionSrcForMomVert_N(bxx, bxy, bxz, + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Centered_4th) { + AdvectionSrcForMomVert_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 5) { + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Upwind_5th) { AdvectionSrcForMomVert_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 6) { - AdvectionSrcForMomVert_N(bxx, bxy, bxz, + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Centered_6th) { + AdvectionSrcForMomVert_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); + vert_adv_type, domhi_z); } else { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } } - } else { + } // end of use_terrain == false + else + { // now do use_terrain == true // Inline with 2nd order for efficiency - if (std::max(horiz_spatial_order,vert_spatial_order) == 2) { + if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) + { ParallelFor(bxx, bxy, bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -264,36 +268,36 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, }); // Template higher order methods } else { - if (horiz_spatial_order == 2) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, + if (horiz_adv_type == AdvType::Centered_2nd) { + AdvectionSrcForMomVert_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 3) { + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Upwind_3rd) { AdvectionSrcForMomVert_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 4) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Centered_4th) { + AdvectionSrcForMomVert_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 5) { + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Upwind_5th) { AdvectionSrcForMomVert_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (horiz_spatial_order == 6) { - AdvectionSrcForMomVert_T(bxx, bxy, bxz, + vert_adv_type, domhi_z); + } else if (horiz_adv_type == AdvType::Centered_6th) { + AdvectionSrcForMomVert_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); + vert_adv_type, domhi_z); } else { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } diff --git a/Source/Advection/AdvectionSrcForMom_N.H b/Source/Advection/AdvectionSrcForMom_N.H index 691db3835..0bdbaf9d7 100644 --- a/Source/Advection/AdvectionSrcForMom_N.H +++ b/Source/Advection/AdvectionSrcForMom_N.H @@ -3,7 +3,7 @@ /** * Function for computing the advective tendency for the x-component of momentum - * without metric terms and for spatial order > 2 + * without metric terms and for higher-order stencils * * @param[in] i,j,k indices of x-face at which to create tendency * @param[in] rho_u x-component of momentum @@ -13,8 +13,6 @@ * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_u map factor on x-faces * @param[in] mf_v map factor on y-faces - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives */ template AMREX_GPU_DEVICE @@ -30,8 +28,6 @@ AdvectionSrcForXMom_N (int i, int j, int k, const amrex::Array4& mf_u, const amrex::Array4& mf_v) { - //used when spatial order > 2 - amrex::Real advectionSrc; auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; @@ -43,8 +39,10 @@ AdvectionSrcForXMom_N (int i, int j, int k, amrex::Real yflux_hi; amrex::Real yflux_lo; amrex::Real zflux_hi; amrex::Real zflux_lo; - amrex::Real mf_u_inv_hi = 1. / mf_u(i+1,j ,0); amrex::Real mf_u_inv_mid = 1. / mf_u(i ,j ,0); amrex::Real mf_u_inv_lo = 1. / mf_u(i-1,j ,0); - amrex::Real mf_v_inv_1 = 1. / mf_v(i ,j+1,0); amrex::Real mf_v_inv_2 = 1. / mf_v(i-1,j+1,0); amrex::Real mf_v_inv_3 = 1. / mf_v(i ,j ,0); amrex::Real mf_v_inv_4 = 1. / mf_v(i-1,j ,0); + amrex::Real mf_u_inv_hi = 1. / mf_u(i+1,j ,0); amrex::Real mf_u_inv_mid = 1. / mf_u(i ,j ,0); + amrex::Real mf_u_inv_lo = 1. / mf_u(i-1,j ,0); + amrex::Real mf_v_inv_1 = 1. / mf_v(i ,j+1,0); amrex::Real mf_v_inv_2 = 1. / mf_v(i-1,j+1,0); + amrex::Real mf_v_inv_3 = 1. / mf_v(i ,j ,0); amrex::Real mf_v_inv_4 = 1. / mf_v(i-1,j ,0); amrex::Real interp_hi(0.), interp_lo(0.); @@ -78,7 +76,7 @@ AdvectionSrcForXMom_N (int i, int j, int k, /** * Function for computing the advective tendency for the y-component of momentum - * without metric terms and for spatial order > 2 + * without metric terms and for higher-order stencils * * @param[in] i,j,k indices of y-face at which to create tendency * @param[in] rho_u x-component of momentum @@ -88,8 +86,6 @@ AdvectionSrcForXMom_N (int i, int j, int k, * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_u map factor on x-faces * @param[in] mf_v map factor on y-faces - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives */ template AMREX_GPU_DEVICE @@ -116,8 +112,10 @@ AdvectionSrcForYMom_N (int i, int j, int k, amrex::Real yflux_hi; amrex::Real yflux_lo; amrex::Real zflux_hi; amrex::Real zflux_lo; - amrex::Real mf_v_inv_hi = 1. / mf_v(i ,j+1,0); amrex::Real mf_v_inv_mid = 1. / mf_v(i ,j ,0); amrex::Real mf_v_inv_lo = 1. / mf_v(i ,j-1,0); - amrex::Real mf_u_inv_1 = 1. / mf_u(i+1,j ,0); amrex::Real mf_u_inv_2 = 1. / mf_u(i+1,j-1,0); amrex::Real mf_u_inv_3 = 1. / mf_u(i ,j ,0); amrex::Real mf_u_inv_4 = 1. / mf_u(i ,j-1,0); + amrex::Real mf_v_inv_hi = 1. / mf_v(i ,j+1,0); amrex::Real mf_v_inv_mid = 1. / mf_v(i ,j ,0); + amrex::Real mf_v_inv_lo = 1. / mf_v(i ,j-1,0); + amrex::Real mf_u_inv_1 = 1. / mf_u(i+1,j ,0); amrex::Real mf_u_inv_2 = 1. / mf_u(i+1,j-1,0); + amrex::Real mf_u_inv_3 = 1. / mf_u(i ,j ,0); amrex::Real mf_u_inv_4 = 1. / mf_u(i ,j-1,0); amrex::Real interp_hi(0.), interp_lo(0.); @@ -151,7 +149,7 @@ AdvectionSrcForYMom_N (int i, int j, int k, /** * Function for computing the advective tendency for the z-component of momentum - * without metric terms and for spatial order > 2 + * without metric terms and for higher-order stencils * * @param[in] i,j,k indices of z-face at which to create tendency * @param[in] rho_u x-component of momentum @@ -162,8 +160,6 @@ AdvectionSrcForYMom_N (int i, int j, int k, * @param[in] mf_m map factor on cell centers * @param[in] mf_u map factor on x-faces * @param[in] mf_v map factor on y-faces - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives * @param[in] domhi_z maximum k value in the domain */ template @@ -182,8 +178,7 @@ AdvectionSrcForZMom_N (int i, int j, int k, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int vert_spatial_order, - const int domhi_z) + const int vert_adv_type, const int domhi_z) { amrex::Real advectionSrc; @@ -214,35 +209,50 @@ AdvectionSrcForZMom_N (int i, int j, int k, yflux_hi = rho_v_avg_hi * interp_hi; yflux_lo = rho_v_avg_lo * interp_lo; - // If k+1 == domhi_z and spatial_order >= 3 we would reach to k = domhi_z+2 so we set to spatial_order = min(spatial_order,2) - // If k+1 == domhi_z-1 and spatial_order >= 5 we would reach to k = domhi_z+2 so we set to spatial_order = min(spatial_order,4) - int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1)); + // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1)); + // If k == domhi_z , l_spatial_order_hi = 2 + // If k == domhi_z-1 , l_spatial_order_hi = std::min(vert_spatial_order, 4); + // If k == 1 , l_spatial_order_hi = std::min(vert_spatial_order, 4); if (k == domhi_z+1) { zflux_hi = rho_w(i,j,k) * w(i,j,k); - } else if (l_spatial_order_hi < vert_spatial_order) { - rho_w_avg_hi = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k+1)); - interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,l_spatial_order_hi); - zflux_hi = rho_w_avg_hi * interp_hi; } else { rho_w_avg_hi = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k+1)); - interp_w_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi); + if (k == domhi_z) + { + interp_w_wall.InterpolateInZ_hi(i,j,k,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); + } else { + interp_w_wall.InterpolateInZ_hi(i,j,k,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); + } zflux_hi = rho_w_avg_hi * interp_hi; } - // If k == 1 and spatial_order >= 3 we would reach to k = -1 so we set to spatial_order = min(spatial_order,2) - // If k == 2 and spatial_order >= 5 we would reach to k = -1 so we set to spatial_order = min(spatial_order,4) - int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(domhi_z+2-k)), 2*k); + // int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(domhi_z+2-k)), 2*k); + // If k == 1 , l_spatial_order_hi = 2 + // If k == 2 , l_spatial_order_hi = std::min(vert_spatial_order, 4); + // If k == domhi_z , l_spatial_order_hi = std::min(vert_spatial_order, 4); if (k == 0) { zflux_lo = rho_w(i,j,k) * w(i,j,k); - } else if (l_spatial_order_lo < vert_spatial_order) { - rho_w_avg_lo = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k-1)); - interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,l_spatial_order_lo); - zflux_lo = rho_w_avg_lo * interp_lo; } else { rho_w_avg_lo = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k-1)); - interp_w_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo); + if (k == 1) { + interp_w_wall.InterpolateInZ_lo(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); + } else { + interp_w_wall.InterpolateInZ_lo(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); + } zflux_lo = rho_w_avg_lo * interp_lo; } @@ -274,7 +284,7 @@ AdvectionSrcForMomWrapper_N(const amrex::Box& bxx, const amrex::Box& bxy, const const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int vert_spatial_order, + const int vert_adv_type, const int domhi_z) { // Instantiate the appropriate structs @@ -299,7 +309,7 @@ AdvectionSrcForMomWrapper_N(const amrex::Box& bxx, const amrex::Box& bxy, const rho_w_rhs(i, j, k) = -AdvectionSrcForZMom_N(i, j, k, rho_u, rho_v, rho_w, w, interp_w_h, interp_w_v, interp_w_wall, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); + vert_adv_type, domhi_z); }); } @@ -322,39 +332,39 @@ AdvectionSrcForMomVert_N(const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int vert_spatial_order, + const int vert_adv_type, const int domhi_z) { - if (vert_spatial_order == 2) { - AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, + if (vert_adv_type == AdvType::Centered_2nd) { + AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, rho_w, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 3) { + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Upwind_3rd) { AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, rho_w, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 4) { - AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Centered_4th) { + AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, rho_w, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 5) { + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Upwind_5th) { AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, rho_w, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 6) { - AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Centered_6th) { + AdvectionSrcForMomWrapper_N(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, rho_w, u, v, w, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); + vert_adv_type, domhi_z); } else { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } diff --git a/Source/Advection/AdvectionSrcForMom_T.H b/Source/Advection/AdvectionSrcForMom_T.H index 301e5d390..a09f2f2f7 100644 --- a/Source/Advection/AdvectionSrcForMom_T.H +++ b/Source/Advection/AdvectionSrcForMom_T.H @@ -16,8 +16,6 @@ * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_u map factor on x-faces * @param[in] mf_v map factor on y-faces - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives */ template AMREX_GPU_DEVICE @@ -110,8 +108,6 @@ AdvectionSrcForXMom_T (int i, int j, int k, * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_u map factor on x-faces * @param[in] mf_v map factor on y-faces - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives */ template AMREX_GPU_DEVICE @@ -204,8 +200,7 @@ AdvectionSrcForYMom_T (int i, int j, int k, * @param[in] mf_m map factor on cell centers * @param[in] mf_u map factor on x-faces * @param[in] mf_v map factor on y-faces - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives + * @param[in] vert_adv_type int that defines advection stencil * @param[in] domhi_z maximum k value in the domain */ template @@ -226,7 +221,7 @@ AdvectionSrcForZMom_T (int i, int j, int k, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int vert_spatial_order, + const int vert_adv_type, const int domhi_z) { amrex::Real advectionSrc; @@ -270,16 +265,24 @@ AdvectionSrcForZMom_T (int i, int j, int k, Omega_avg_hi = (k == domhi_z+1) ? Omega(i,j,k) : 0.5 * (Omega(i,j,k) + Omega(i,j,k+1)); amrex::Real centFluxZZNext = Omega_avg_hi; - // If k == domhi_z and spatial_order >= 3 we would reach to k = domhi_z+2 so we set to spatial_order = min(spatial_order,2) - // If k == domhi_z-1 and spatial_order >= 5 we would reach to k = domhi_z+2 so we set to spatial_order = min(spatial_order,4) - int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1)); + // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1)); + // If k == domhi_z , l_spatial_order_hi = 2 + // If k == domhi_z-1 , l_spatial_order_hi = std::min(vert_spatial_order, 4); + // If k == 1 , l_spatial_order_hi = std::min(vert_spatial_order, 4); if (k == domhi_z+1) { centFluxZZNext *= w(i,j,k); - } else if (l_spatial_order_hi < vert_spatial_order) { - interp_omega_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi,l_spatial_order_hi); - centFluxZZNext *= interp_hi; } else { - interp_omega_v.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi); + if (k == domhi_z) { + interp_omega_wall.InterpolateInZ_hi(i,j,k,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); + } else { + interp_omega_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi,vert_adv_type); + } + } else { + interp_omega_v.InterpolateInZ_hi(i,j,k,0,interp_hi,Omega_avg_hi); + } centFluxZZNext *= interp_hi; } @@ -288,16 +291,24 @@ AdvectionSrcForZMom_T (int i, int j, int k, Omega_avg_lo = (k == 0) ? Omega(i,j,k) : 0.5 * (Omega(i,j,k) + Omega(i,j,k-1)); amrex::Real centFluxZZPrev = Omega_avg_lo; - // If k == 1 and spatial_order >= 3 we would reach to k = -1 so we set to spatial_order = min(spatial_order,2) - // If k == 2 and spatial_order >= 5 we would reach to k = -1 so we set to spatial_order = min(spatial_order,4) - int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(domhi_z+2-k)), 2*k); + // int l_spatial_order_lo = std::min(std::min(vert_spatial_order, 2*(domhi_z+2-k)), 2*k); + // If k == domhi_z , l_spatial_order_hi = 2 + // If k == domhi_z-1 , l_spatial_order_hi = std::min(vert_spatial_order, 4); + // If k == 1 , l_spatial_order_hi = std::min(vert_spatial_order, 4); if (k == 0) { centFluxZZPrev *= w(i,j,k); - } else if (l_spatial_order_lo < vert_spatial_order) { - interp_omega_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo,l_spatial_order_lo); - centFluxZZPrev *= interp_lo; } else { - interp_omega_v.InterpolateInZ_lo(i,j,k,0,interp_lo,Omega_avg_lo); + if (k == 1) { + interp_omega_wall.InterpolateInZ_lo(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); + } else { + interp_omega_wall.InterpolateInZ_lo(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); + } centFluxZZPrev *= interp_lo; } @@ -336,7 +347,7 @@ AdvectionSrcForMomWrapper_T(const amrex::Box& bxx, const amrex::Box& bxy, const const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int vert_spatial_order, + const int vert_adv_type, const int domhi_z) { // Instantiate the appropriate structs @@ -361,7 +372,7 @@ AdvectionSrcForMomWrapper_T(const amrex::Box& bxx, const amrex::Box& bxy, const rho_w_rhs(i, j, k) = -AdvectionSrcForZMom_T(i, j, k, rho_u, rho_v, Omega, w, z_nd, detJ, interp_w_h, interp_w_v, interp_w_wall, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); + vert_adv_type, domhi_z); }); } @@ -386,39 +397,38 @@ AdvectionSrcForMomVert_T(const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int vert_spatial_order, - const int domhi_z) + const int vert_adv_type, const int domhi_z) { - if (vert_spatial_order == 2) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 3) { + if (vert_adv_type == AdvType::Centered_2nd) { + AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, z_nd, detJ, + cellSizeInv, mf_m, mf_u, mf_v, + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Upwind_3rd) { AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 4) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 5) { + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Centered_4th) { + AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, z_nd, detJ, + cellSizeInv, mf_m, mf_u, mf_v, + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Upwind_5th) { AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, rho_u, rho_v, Omega, u, v, w, z_nd, detJ, cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); - } else if (vert_spatial_order == 6) { - AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, - rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, z_nd, detJ, - cellSizeInv, mf_m, mf_u, mf_v, - vert_spatial_order, domhi_z); + vert_adv_type, domhi_z); + } else if (vert_adv_type == AdvType::Centered_6th) { + AdvectionSrcForMomWrapper_T(bxx, bxy, bxz, + rho_u_rhs, rho_v_rhs, rho_w_rhs, + rho_u, rho_v, Omega, u, v, w, z_nd, detJ, + cellSizeInv, mf_m, mf_u, mf_v, + vert_adv_type, domhi_z); } else { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index 530099052..63d82adf1 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -29,10 +29,8 @@ using namespace amrex; * @param[in] mf_m map factor at cell centers * @param[in] mf_u map factor at x-faces * @param[in] mf_v map factor at y-faces - * @param[in] all_use_WENO whether all variables (or just moisture variables) use WENO advection scheme - * @param[in] spatial_order_WENO spatial order if using WENO (3,5, or 7) - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives if not using WENO (2-6) - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives if not using WENO (2-6) + * @param[in] horiz_adv_type advection scheme to be used in horiz. directions for dry scalars + * @param[in] vert_adv_type advection scheme to be used in horiz. directions for dry scalars * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms) */ @@ -51,8 +49,8 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, const Array4& mf_m, const Array4& mf_u, const Array4& mf_v, - const int horiz_spatial_order, - const int vert_spatial_order, + const int horiz_adv_type, + const int vert_adv_type, const int use_terrain) { BL_PROFILE_VAR("AdvectionSrcForRhoAndTheta", AdvectionSrcForRhoAndTheta); @@ -63,7 +61,8 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, if (!use_terrain) { // Inline with 2nd order for efficiency - if (std::max(horiz_spatial_order,vert_spatial_order) == 2) { + if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) + { ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real xflux_lo = rho_u(i ,j,k) / mf_u(i ,j ,0); @@ -101,36 +100,36 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, }); // Template higher order methods } else { - if (horiz_spatial_order == 2) { - 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_spatial_order); - } else if (horiz_spatial_order == 3) { + if (horiz_adv_type == 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) { 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_spatial_order); - } else if (horiz_spatial_order == 4) { - 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_spatial_order); - } else if (horiz_spatial_order == 5) { + vert_adv_type); + } else if (horiz_adv_type == 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) { 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_spatial_order); - } else if (horiz_spatial_order == 6) { - 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_spatial_order); + vert_adv_type); + } else if (horiz_adv_type == 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 { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } @@ -138,7 +137,8 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, } else { // Inline with 2nd order for efficiency - if (std::max(horiz_spatial_order,vert_spatial_order) == 2) { + if (horiz_adv_type == AdvType::Centered_2nd && vert_adv_type == AdvType::Centered_2nd) + { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real invdetJ = 1./ detJ(i,j,k); @@ -188,36 +188,36 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, }); // Template higher order methods (horizontal first) } else { - if (horiz_spatial_order == 2) { - 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_spatial_order); - } else if (horiz_spatial_order == 3) { + if (horiz_adv_type == 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) { 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_spatial_order); - } else if (horiz_spatial_order == 4) { - 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_spatial_order); - } else if (horiz_spatial_order == 5) { + mf_u, mf_v, vert_adv_type); + } else if (horiz_adv_type == 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) { 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_spatial_order); - } else if (horiz_spatial_order == 6) { - 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_spatial_order); + mf_u, mf_v, vert_adv_type); + } else if (horiz_adv_type == 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 { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } @@ -242,16 +242,13 @@ AdvectionSrcForRhoAndTheta (const Box& bx, const Box& valid_bx, * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing * @param[in] mf_m map factor at cell centers - * @param[in] all_use_WENO whether all variables use WENO advection scheme - * @param[in] moist_use_WENO whether moisture variables use WENO advection scheme - * @param[in] spatial_order_WENO spatial order if using WENO (3,5, or 7) - * @param[in] horiz_spatial_order spatial order to be used for lateral derivatives if not using WENO (2-6) - * @param[in] vert_spatial_order spatial order to be used for vertical derivatives if not using WENO (2-6) + * @param[in] horiz_adv_type advection scheme to be used in horiz. directions for dry scalars + * @param[in] vert_adv_type advection scheme to be used in horiz. directions for dry scalars * @param[in] use_terrain if true, use the terrain-aware derivatives (with metric terms) */ void -AdvectionSrcForScalars (const Box& bx, const int &icomp, const int &ncomp, +AdvectionSrcForScalars (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, @@ -259,58 +256,17 @@ AdvectionSrcForScalars (const Box& bx, const int &icomp, const int &ncomp, const Array4& detJ, const GpuArray& cellSizeInv, const Array4& mf_m, - const bool all_use_WENO, - const bool moist_use_WENO, - const bool all_use_WENO_Z, - const bool moist_use_WENO_Z, - const int spatial_order_WENO, - const int horiz_spatial_order, - const int vert_spatial_order, + const int horiz_adv_type, + const int vert_adv_type, const int use_terrain) { BL_PROFILE_VAR("AdvectionSrcForScalars", AdvectionSrcForScalars); auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; - int ncomp_end{ncomp}; - int moist_off = RhoScalar_comp; -#if defined(ERF_USE_MOISTURE) - moist_off = RhoQt_comp; -#elif defined(ERF_USE_WARM_NO_PRECIP) - moist_off = RhoQv_comp; -#endif - - // Running with WENO for moisture but not for other vars - if((moist_use_WENO || moist_use_WENO_Z) && ((icomp+ncomp)==NVAR) ) { - ncomp_end -= 2; - if (moist_use_WENO && spatial_order_WENO==3) { - AdvectionSrcForScalarsWrapper_N(bx, 2, moist_off, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (moist_use_WENO && spatial_order_WENO==5) { - AdvectionSrcForScalarsWrapper_N(bx, 2, moist_off, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (moist_use_WENO_Z && spatial_order_WENO==3) { - AdvectionSrcForScalarsWrapper_N(bx, 2, moist_off, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (moist_use_WENO_Z && spatial_order_WENO==5) { - AdvectionSrcForScalarsWrapper_N(bx, 2, moist_off, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else { - AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); - } - } - // Inline with 2nd order for efficiency - bool not_WENO = !(all_use_WENO || all_use_WENO_Z); - if (std::max(horiz_spatial_order,vert_spatial_order) == 2 && not_WENO) { - amrex::ParallelFor(bx, ncomp_end, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + 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 { Real invdetJ = (use_terrain) ? 1. / detJ(i,j,k) : 1.; @@ -332,48 +288,53 @@ AdvectionSrcForScalars (const Box& bx, const int &icomp, const int &ncomp, }); // Template higher order methods (horizontal first) } else { - if (horiz_spatial_order == 2 && not_WENO) { - AdvectionSrcForScalarsVert_N(bx, ncomp_end, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_spatial_order); - } else if (horiz_spatial_order == 3 && not_WENO) { - AdvectionSrcForScalarsVert_N(bx, ncomp_end, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_spatial_order); - } else if (horiz_spatial_order == 4 && not_WENO) { - AdvectionSrcForScalarsVert_N(bx, ncomp_end, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m, vert_spatial_order); - } else if (horiz_spatial_order == 5 && not_WENO) { - AdvectionSrcForScalarsVert_N(bx, ncomp_end, icomp, + 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_spatial_order); - } else if (horiz_spatial_order == 6 && not_WENO) { - AdvectionSrcForScalarsVert_N(bx, ncomp_end, icomp, + 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_spatial_order); - } else if (all_use_WENO && spatial_order_WENO==3) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, + 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 (all_use_WENO && spatial_order_WENO==5) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, + } 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 (all_use_WENO_Z && spatial_order_WENO==3) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, + } 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 (all_use_WENO_Z && spatial_order_WENO==5) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, + } 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); diff --git a/Source/Advection/AdvectionSrcForState_N.H b/Source/Advection/AdvectionSrcForState_N.H index 7ebedd681..8a04e74c4 100644 --- a/Source/Advection/AdvectionSrcForState_N.H +++ b/Source/Advection/AdvectionSrcForState_N.H @@ -91,33 +91,33 @@ AdvectionSrcForRhoThetaVert_N(const amrex::Box& bx, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int& vert_spatial_order) + const int vert_adv_type) { - if (vert_spatial_order == 2) { - 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_spatial_order == 3) { + if (vert_adv_type == 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) { 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_spatial_order == 4) { - 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_spatial_order == 5) { + } else if (vert_adv_type == 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) { 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_spatial_order == 6) { - 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) { + 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 { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } @@ -129,9 +129,7 @@ AdvectionSrcForRhoThetaVert_N(const amrex::Box& bx, template void AdvectionSrcForScalarsWrapper_N(const amrex::Box& bx, - const int& ncomp_end, - const int& icomp, - const int& use_terrain, + 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, @@ -146,7 +144,7 @@ AdvectionSrcForScalarsWrapper_N(const amrex::Box& bx, InterpType_V interp_prim_v(cell_prim); auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2]; - amrex::ParallelFor(bx, ncomp_end, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + 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); @@ -180,9 +178,7 @@ AdvectionSrcForScalarsWrapper_N(const amrex::Box& bx, template void AdvectionSrcForScalarsVert_N(const amrex::Box& bx, - const int& ncomp_end, - const int& icomp, - const int& use_terrain, + 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, @@ -191,33 +187,33 @@ AdvectionSrcForScalarsVert_N(const amrex::Box& bx, const amrex::Array4& detJ, const amrex::GpuArray& cellSizeInv, const amrex::Array4& mf_m, - const int& vert_spatial_order) + const int vert_adv_type) { - if (vert_spatial_order == 2) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_spatial_order == 3) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_spatial_order == 4) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, - use_terrain, advectionSrc, cell_prim, - avg_xmom, avg_ymom, avg_zmom, detJ, - cellSizeInv, mf_m); - } else if (vert_spatial_order == 5) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, + 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_spatial_order == 6) { - AdvectionSrcForScalarsWrapper_N(bx, ncomp_end, icomp, + } 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 b6d7d7004..0b1384f3d 100644 --- a/Source/Advection/AdvectionSrcForState_T.H +++ b/Source/Advection/AdvectionSrcForState_T.H @@ -108,38 +108,38 @@ AdvectionSrcForRhoThetaVert_T(const amrex::Box& bx, const amrex::Array4& mf_m, const amrex::Array4& mf_u, const amrex::Array4& mf_v, - const int& vert_spatial_order) + const int vert_adv_type) { - if (vert_spatial_order == 2) { - AdvectionSrcForRhoThetaWrapper_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); - } else if (vert_spatial_order == 3) { + if (vert_adv_type == AdvType::Centered_2nd) { + AdvectionSrcForRhoThetaWrapper_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); + } else if (vert_adv_type == AdvType::Upwind_3rd) { AdvectionSrcForRhoThetaWrapper_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); - } else if (vert_spatial_order == 4) { - AdvectionSrcForRhoThetaWrapper_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); - } else if (vert_spatial_order == 5) { + } else if (vert_adv_type == AdvType::Centered_4th) { + AdvectionSrcForRhoThetaWrapper_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); + } else if (vert_adv_type == AdvType::Upwind_5th) { AdvectionSrcForRhoThetaWrapper_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); - } else if (vert_spatial_order == 6) { - AdvectionSrcForRhoThetaWrapper_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); + } else if (vert_adv_type == AdvType::Centered_6th) { + AdvectionSrcForRhoThetaWrapper_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); } else { AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!"); } diff --git a/Source/DataStruct.H b/Source/DataStruct.H index 527e9337d..5ff9876a1 100644 --- a/Source/DataStruct.H +++ b/Source/DataStruct.H @@ -10,6 +10,7 @@ #include #include +#include enum class ABLDriverType { None, PressureGradient, GeostrophicWind @@ -201,9 +202,109 @@ struct SolverChoice { amrex::Abort("If you specify incompressible, you must specific no_substepping"); } - // Order of spatial discretization - pp.query("horiz_spatial_order", horiz_spatial_order); - pp.query("vert_spatial_order", vert_spatial_order); + // Order and type of spatial discretizations used in advection + std::string dycore_horiz_adv_string = "" ; std::string dycore_vert_adv_string = ""; + std::string dryscal_horiz_adv_string = "" ; std::string dryscal_vert_adv_string = ""; + pp.query("dycore_horiz_adv_type" , dycore_horiz_adv_string); + pp.query("dycore_vert_adv_type" , dycore_vert_adv_string); + pp.query("dryscal_horiz_adv_type" , dryscal_horiz_adv_string); + pp.query("dryscal_vert_adv_type" , dryscal_vert_adv_string); + +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + std::string moistscal_horiz_adv_string = ""; std::string moistscal_vert_adv_string = ""; + pp.query("moistscal_horiz_adv_type", moistscal_horiz_adv_string); + pp.query("moistscal_vert_adv_type" , moistscal_vert_adv_string); +#endif + + if ( (dycore_horiz_adv_string == "Centered_2nd") || + (dycore_horiz_adv_string == "Upwind_3rd" ) || + (dycore_horiz_adv_string == "Centered_4th") || + (dycore_horiz_adv_string == "Upwind_5th" ) || + (dycore_horiz_adv_string == "Centered_6th") ) + { + dycore_horiz_adv_type = adv_type_convert_string_to_int(dycore_horiz_adv_string); + amrex::Print() << "Using dycore_horiz_adv_type: " << dycore_horiz_adv_string << std::endl; + } else { + amrex::Warning("Using default dycore_horiz_adv_type"); + } + + if ( (dycore_vert_adv_string == "Centered_2nd") || + (dycore_vert_adv_string == "Upwind_3rd" ) || + (dycore_vert_adv_string == "Centered_4th") || + (dycore_vert_adv_string == "Upwind_5th" ) || + (dycore_vert_adv_string == "Centered_6th") ) + { + dycore_vert_adv_type = adv_type_convert_string_to_int(dycore_vert_adv_string); + amrex::Print() << "Using dycore_vert_adv_type: " << dycore_vert_adv_string << std::endl; + } else { + amrex::Warning("Using default dycore_vert_adv_type"); + } + + if ( (dryscal_horiz_adv_string == "Centered_2nd") || + (dryscal_horiz_adv_string == "Upwind_3rd" ) || + (dryscal_horiz_adv_string == "Centered_4th") || + (dryscal_horiz_adv_string == "Upwind_5th" ) || + (dryscal_horiz_adv_string == "Centered_6th") || + (dryscal_horiz_adv_string == "Weno_3" ) || + (dryscal_horiz_adv_string == "Weno_3Z" ) || + (dryscal_horiz_adv_string == "Weno_3MZQ" ) || + (dryscal_horiz_adv_string == "Weno_5Z" ) ) + { + dryscal_horiz_adv_type = adv_type_convert_string_to_int(dryscal_horiz_adv_string); + amrex::Print() << "Using dryscal_horiz_adv_type: " << dryscal_horiz_adv_string << std::endl; + } else { + amrex::Warning("Using default dryscal_horiz_adv_type"); + } + + if ( (dryscal_vert_adv_string == "Centered_2nd") || + (dryscal_vert_adv_string == "Upwind_3rd" ) || + (dryscal_vert_adv_string == "Centered_4th") || + (dryscal_vert_adv_string == "Upwind_5th" ) || + (dryscal_vert_adv_string == "Centered_6th") || + (dryscal_vert_adv_string == "Weno_3" ) || + (dryscal_vert_adv_string == "Weno_3Z" ) || + (dryscal_vert_adv_string == "Weno_3MZQ" ) || + (dryscal_vert_adv_string == "Weno_5Z" ) ) + { + dryscal_vert_adv_type = adv_type_convert_string_to_int(dryscal_vert_adv_string); + amrex::Print() << "Using dryscal_vert_adv_type: " << dryscal_vert_adv_string << std::endl; + } else { + amrex::Warning("Using default dryscal_vert_adv_type"); + } + +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + if ( (moistscal_horiz_adv_string == "Centered_2nd") || + (moistscal_horiz_adv_string == "Upwind_3rd" ) || + (moistscal_horiz_adv_string == "Centered_4th") || + (moistscal_horiz_adv_string == "Upwind_5th" ) || + (moistscal_horiz_adv_string == "Centered_6th") || + (moistscal_horiz_adv_string == "Weno_3" ) || + (moistscal_horiz_adv_string == "Weno_3Z" ) || + (moistscal_horiz_adv_string == "Weno_3MZQ" ) || + (moistscal_horiz_adv_string == "Weno_5Z" ) ) + { + moistscal_horiz_adv_type = adv_type_convert_string_to_int(moistscal_horiz_adv_string); + amrex::Print() << "Using moistscal_horiz_adv_type: " << moistscal_horiz_adv_string << std::endl; + } else { + amrex::Warning("Using default moistscal_horiz_adv_type"); + } + + if ( (moistscal_vert_adv_string == "Centered_2nd") || + (moistscal_vert_adv_string == "Upwind_3rd" ) || + (moistscal_vert_adv_string == "Centered_4th") || + (moistscal_vert_adv_string == "Upwind_5th" ) || + (moistscal_vert_adv_string == "Centered_6th") || + (moistscal_vert_adv_string == "Weno_3" ) || + (moistscal_vert_adv_string == "Weno_3Z" ) || + (moistscal_vert_adv_string == "Weno_3MZQ" ) || + (moistscal_vert_adv_string == "Weno_5Z" ) ) + { + moistscal_vert_adv_type = adv_type_convert_string_to_int(moistscal_vert_adv_string); + amrex::Print() << "Using moistscal_vert_adv_type: " << moistscal_vert_adv_string << std::endl; + } else { + amrex::Warning("Using default moistscal_vert_adv_type"); + } +#endif // Include Coriolis forcing? pp.query("use_coriolis", use_coriolis); @@ -263,37 +364,13 @@ struct SolverChoice { build_coriolis_forcings(); } - static std::string ic_bc_type_string = "Ideal"; - pp.query("ic_bc_type", ic_bc_type_string); - pp.query("Ave_Plane", ave_plane); #ifdef ERF_USE_MOISTURE pp.query("mp_clouds", do_cloud); pp.query("mp_precip", do_precip); - pp.query("moist_use_WENO_Z", moist_use_WENO_Z); - if(!moist_use_WENO_Z) moist_use_WENO = true; #endif - // Use WENO for advection? - pp.query("all_use_WENO", all_use_WENO); - if (all_use_WENO || moist_use_WENO) { - pp.query("spatial_order_WENO", spatial_order_WENO); - AMREX_ASSERT_WITH_MESSAGE(( (spatial_order_WENO==3) || (spatial_order_WENO==5) ), - "WENO advection only supports orders 3 & 5."); - } - - // Use WENO-Z for advection? - pp.query("all_use_WENO_Z", all_use_WENO_Z); - if (all_use_WENO_Z || moist_use_WENO_Z) { - pp.query("spatial_order_WENO", spatial_order_WENO); - AMREX_ASSERT_WITH_MESSAGE(( (spatial_order_WENO==3) || (spatial_order_WENO==5) ), - "WENO advection only supports orders 3 & 5."); - AMREX_ASSERT_WITH_MESSAGE(( (!all_use_WENO) || (!moist_use_WENO) ), - "WENO and WENO-Z cannot be specified together."); - } - - // Use numerical diffusion? pp.query("use_NumDiff",use_NumDiff); if(use_NumDiff) { @@ -327,8 +404,14 @@ struct SolverChoice { amrex::Print() << "sigma_k : " << sigma_k << std::endl; amrex::Print() << "Pr_t : " << Pr_t << std::endl; amrex::Print() << "Sc_t : " << Sc_t << std::endl; - amrex::Print() << "horiz spatial_order : " << horiz_spatial_order << std::endl; - amrex::Print() << "vert spatial_order : " << vert_spatial_order << std::endl; + amrex::Print() << "dycore_horiz_adv_type : " << adv_type_convert_int_to_string(dycore_horiz_adv_type) << std::endl; + amrex::Print() << "dycore_vert_adv_type : " << adv_type_convert_int_to_string(dycore_vert_adv_type) << std::endl; + amrex::Print() << "dryscal_horiz_adv_type : " << adv_type_convert_int_to_string(dryscal_horiz_adv_type) << std::endl; + amrex::Print() << "dryscal_vert_adv_type : " << adv_type_convert_int_to_string(dryscal_vert_adv_type) << std::endl; +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + amrex::Print() << "moistscal_horiz_adv_type : " << adv_type_convert_int_to_string(moistscal_horiz_adv_type) << std::endl; + amrex::Print() << "moistscal_vert_adv_type : " << adv_type_convert_int_to_string(moistscal_vert_adv_type) << std::endl; +#endif if (abl_driver_type == ABLDriverType::None) { amrex::Print() << "ABL Driver Type: " << "None" << std::endl; @@ -397,6 +480,60 @@ struct SolverChoice { } } + std::string + adv_type_convert_int_to_string (int adv_int) + { + if (adv_int == AdvType::Centered_2nd) { + return "Centered_2nd"; + } else if (adv_int == AdvType::Upwind_3rd) { + return "Upwind_3rd"; + } else if (adv_int == AdvType::Centered_4th) { + return "Centered_4th"; + } else if (adv_int == AdvType::Upwind_5th) { + return "Upwind_5th"; + } else if (adv_int == AdvType::Centered_6th) { + return "Centered_6th"; + } else if (adv_int == AdvType::Weno_3) { + return "Weno_3"; + } else if (adv_int == AdvType::Weno_3Z) { + return "Weno_3Z"; + } else if (adv_int == AdvType::Weno_5) { + return "Weno_5"; + } else if (adv_int == AdvType::Weno_5Z) { + return "Weno_5Z"; + } else if (adv_int == AdvType::Weno_3MZQ) { + return "Weno_3MZQ"; + } + } + + + int adv_type_convert_string_to_int (std::string adv_string) + { + if (adv_string == "Centered_2nd") { + return AdvType::Centered_2nd; + } else if (adv_string == "Upwind_3rd") { + return AdvType::Upwind_3rd; + } else if (adv_string == "Centered_4th") { + return AdvType::Centered_4th; + } else if (adv_string == "Upwind_5th") { + return AdvType::Upwind_5th; + } else if (adv_string == "Centered_6th") { + return AdvType::Centered_6th; + } else if (adv_string == "Weno_3") { + return AdvType::Weno_3; + } else if (adv_string == "Weno_3Z") { + return AdvType::Weno_3Z; + } else if (adv_string == "Weno_5") { + return AdvType::Weno_5; + } else if (adv_string == "Weno_5Z") { + return AdvType::Weno_5Z; + } else if (adv_string == "Weno_3MZQ") { + return AdvType::Weno_3MZQ; + } else { + return -1; + } + } + // Default prefix std::string pp_prefix {"erf"}; @@ -509,15 +646,15 @@ struct SolverChoice { amrex::Real sinphi = 0.0; // Spatial discretization - int horiz_spatial_order = 2; - int vert_spatial_order = 2; - - // Positive definite advection scheme (3/5 order) - bool all_use_WENO{false}; - bool moist_use_WENO{false}; - bool all_use_WENO_Z{false}; - bool moist_use_WENO_Z{false}; - int spatial_order_WENO{3}; + // Order and type of spatial discretizations used in advection + int dycore_horiz_adv_type = AdvType::Centered_2nd; + int dycore_vert_adv_type = AdvType::Centered_2nd; + int dryscal_horiz_adv_type = AdvType::Centered_2nd; + int dryscal_vert_adv_type = AdvType::Centered_2nd; +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + int moistscal_horiz_adv_type = AdvType::Weno_3; + int moistscal_vert_adv_type = AdvType::Weno_3; +#endif // Numerical diffusion bool use_NumDiff{false}; diff --git a/Source/ERF.H b/Source/ERF.H index 22b6d4140..5dfd1a0cc 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -742,29 +742,47 @@ private: static AMREX_FORCE_INLINE int - ComputeGhostCells (const int& horiz_spatial_order, - const int& vert_spatial_order, - const bool& use_numerical_diff) + ComputeGhostCells (const SolverChoice& solverChoice) { - int spatial_order = std::max(horiz_spatial_order, vert_spatial_order); - - int nGhostCells = 0; - if (use_numerical_diff) { - nGhostCells = 3; - } else { - switch (spatial_order) { - case 2: case 3: case 4: - nGhostCells = 2; // We need this many to compute the eddy viscosity in the ghost cells - break; - case 5: case 6: - nGhostCells = 3; - break; - default: - amrex::Error("Must specify spatial order to be 2,3,4,5 or 6"); - } - } - - return nGhostCells; + if (solverChoice.use_NumDiff) + { + return 3; + } else { + if ( + (solverChoice.dycore_horiz_adv_type == AdvType::Centered_6th) + || (solverChoice.dycore_vert_adv_type == AdvType::Centered_6th) +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + || (solverChoice.moistscal_horiz_adv_type == AdvType::Centered_6th) + || (solverChoice.moistscal_horiz_adv_type == AdvType::Centered_6th) +#endif + || (solverChoice.dryscal_horiz_adv_type == AdvType::Centered_6th) + || (solverChoice.dryscal_vert_adv_type == AdvType::Centered_6th) ) + { return 3; } + else if ( + (solverChoice.dycore_horiz_adv_type == AdvType::Upwind_5th) + || (solverChoice.dycore_vert_adv_type == AdvType::Upwind_5th) +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + || (solverChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th) + || (solverChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th) +#endif + || (solverChoice.dryscal_horiz_adv_type == AdvType::Upwind_5th) + || (solverChoice.dryscal_vert_adv_type == AdvType::Upwind_5th) ) + { return 3; } + else if ( + (solverChoice.dryscal_horiz_adv_type == AdvType::Weno_5) + || (solverChoice.dryscal_vert_adv_type == AdvType::Weno_5) +#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) + || (solverChoice.moistscal_horiz_adv_type == AdvType::Weno_5) + || (solverChoice.moistscal_vert_adv_type == AdvType::Weno_5) + || (solverChoice.moistscal_horiz_adv_type == AdvType::Weno_5Z) + || (solverChoice.moistscal_vert_adv_type == AdvType::Weno_5Z) +#endif + || (solverChoice.dryscal_horiz_adv_type == AdvType::Weno_5Z) + || (solverChoice.dryscal_vert_adv_type == AdvType::Weno_5Z) ) + { return 3; } + else + { return 2; } + } } AMREX_FORCE_INLINE diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index dd41c539d..667e81134 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -36,12 +36,8 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, // The number of ghost cells for density must be 1 greater than that for velocity // so that we can go back in forth betwen velocity and momentum on all faces - int ngrow_state = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff)+1; - int ngrow_vels = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff); + int ngrow_state = ComputeGhostCells(solverChoice) + 1; + int ngrow_vels = ComputeGhostCells(solverChoice); auto& lev_new = vars_new[lev]; auto& lev_old = vars_old[lev]; @@ -148,9 +144,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, ba_nd.surroundingNodes(); // We need this to be one greater than the ghost cells to handle levels > 0 - int ngrow = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff)+2; + int ngrow = ComputeGhostCells(solverChoice) + 2; z_phys_nd[lev] = std::make_unique(ba_nd,dm,1,IntVect(ngrow,ngrow,1)); if (solverChoice.terrain_type > 0) { z_phys_nd_new[lev] = std::make_unique(ba_nd,dm,1,IntVect(ngrow,ngrow,1)); @@ -282,12 +276,8 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp Vector temp_lev_new(Vars::NumTypes); Vector temp_lev_old(Vars::NumTypes); - int ngrow_state = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff)+1; - int ngrow_vels = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff); + int ngrow_state = ComputeGhostCells(solverChoice) + 1; + int ngrow_vels = ComputeGhostCells(solverChoice); temp_lev_new[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state); temp_lev_old[Vars::cons].define(ba, dm, Cons::NumVars, ngrow_state); diff --git a/Source/IO/NCCheckpoint.cpp b/Source/IO/NCCheckpoint.cpp index e403d40a9..7dc281567 100644 --- a/Source/IO/NCCheckpoint.cpp +++ b/Source/IO/NCCheckpoint.cpp @@ -153,12 +153,9 @@ ERF::ReadNCCheckpointFile () ncf.var("dt") .get(dt.data(), {0}, {static_cast(ndt)}); ncf.var("t_new").get(t_new.data(), {0}, {static_cast(ntime)}); - int ngrow_state = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff)+1; - int ngrow_vels = ComputeGhostCells(solverChoice.horiz_spatial_order, - solverChoice.vert_spatial_order, - solverChoice.use_NumDiff); + int ngrow_state = ComputeGhostCells(solverChoice) + 1; + int ngrow_vels = ComputeGhostCells(solverChoice); + for (int lev = 0; lev <= finest_level; ++lev) { int num_box = static_cast(ncf.dim("Nbox_"+std::to_string(lev)).len()); diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index b6e660034..3fc90bfdd 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -164,7 +164,8 @@ enum struct ERF_BC { // in amrex/Src/Base/AMReX_BC_TYPES.H. We had extras at // the end to use locally namespace ERFBCType { -enum mathematicalBndryTypes : int { bogus = -666, +enum mathematicalBndryTypes : int { + bogus = -666, reflect_odd = -1, int_dir = 0, reflect_even = 1, @@ -185,4 +186,17 @@ namespace IntVar { NumVars }; } + +enum AdvType : int { + Centered_2nd = 101, + Upwind_3rd = 102, + Centered_4th = 103, + Upwind_5th = 104, + Centered_6th = 105, + Weno_3 = 106, + Weno_3Z = 107, + Weno_5 = 108, + Weno_5Z = 109, + Weno_3MZQ = 110 +}; #endif diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 375c5a0eb..20621350e 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -116,8 +116,8 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, int num_comp = 2; int end_comp = start_comp + num_comp - 1; - const int l_horiz_spatial_order = solverChoice.horiz_spatial_order; - const int l_vert_spatial_order = solverChoice.vert_spatial_order; + const int l_horiz_adv_type = solverChoice.dycore_horiz_adv_type; + const int l_vert_adv_type = solverChoice.dycore_vert_adv_type; const bool l_use_terrain = solverChoice.use_terrain; AMREX_ALWAYS_ASSERT (!l_use_terrain); @@ -595,7 +595,7 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, avg_xmom, avg_ymom, avg_zmom, // these are being defined from the rho fluxes cell_prim, z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, - l_horiz_spatial_order, l_vert_spatial_order, l_use_terrain); + horiz_adv_type, vert_adv_type, l_use_terrain); if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); @@ -665,7 +665,7 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, rho_u_rhs, rho_v_rhs, rho_w_rhs, u, v, w, rho_u , rho_v , omega_arr, z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, - l_horiz_spatial_order, l_vert_spatial_order, l_use_terrain, domhi_z); + horiz_adv_type, vert_adv_type, l_use_terrain, domhi_z); if (l_use_diff) { DiffusionSrcForMom_N(tbx, tby, tbz, diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index ba7ff426a..60db0b0de 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -89,8 +89,6 @@ 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 int l_horiz_spatial_order = solverChoice.horiz_spatial_order; - const int l_vert_spatial_order = solverChoice.vert_spatial_order; const bool l_use_terrain = solverChoice.use_terrain; const bool l_moving_terrain = (solverChoice.terrain_type == 1); if (l_moving_terrain) AMREX_ALWAYS_ASSERT(l_use_terrain); @@ -104,11 +102,6 @@ void erf_slow_rhs_post (int /*level*/, const bool l_use_turb = ( solverChoice.les_type == LESType::Smagorinsky || solverChoice.les_type == LESType::Deardorff || solverChoice.pbl_type == PBLType::MYNN25 ); - const bool l_all_WENO = solverChoice.all_use_WENO; - const bool l_moist_WENO = solverChoice.moist_use_WENO; - const bool l_all_WENO_Z = solverChoice.all_use_WENO_Z; - const bool l_moist_WENO_Z = solverChoice.moist_use_WENO_Z; - const int l_spatial_order_WENO = solverChoice.spatial_order_WENO; const amrex::BCRec* bc_ptr = domain_bcs_type_d.data(); @@ -241,27 +234,47 @@ void erf_slow_rhs_post (int /*level*/, 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, l_all_WENO, l_moist_WENO, l_all_WENO_Z, l_moist_WENO_Z, - l_spatial_order_WENO,l_horiz_spatial_order, l_vert_spatial_order, + cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + solverChoice.dryscal_horiz_adv_type, + solverChoice.dryscal_vert_adv_type, l_use_terrain); } 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, l_all_WENO, l_moist_WENO, l_all_WENO_Z, l_moist_WENO_Z, - l_spatial_order_WENO,l_horiz_spatial_order, l_vert_spatial_order, + cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + solverChoice.dryscal_horiz_adv_type, + solverChoice.dryscal_vert_adv_type, l_use_terrain); } + + // This is simply an advected scalar for convenience start_comp = RhoScalar_comp; - num_comp = S_data[IntVar::cons].nComp() - start_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, + solverChoice.dryscal_horiz_adv_type, + solverChoice.dryscal_vert_adv_type, + l_use_terrain); + +#ifdef ERF_USE_MOISTURE + start_comp = RhoQt_comp; + num_comp = 2; AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, - cur_prim, cell_rhs, detJ_arr, - dxInv, mf_m, l_all_WENO, l_moist_WENO, l_all_WENO_Z, l_moist_WENO_Z, - l_spatial_order_WENO,l_horiz_spatial_order, l_vert_spatial_order, + cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + solverChoice.moistscal_horiz_adv_type, + solverChoice.moistscal_vert_adv_type, l_use_terrain); +#elif defined(ERF_USE_WARM_NO_PRECIP) + start_comp = RhoQv_comp; + num_comp = 2; + AdvectionSrcForScalars(tbx, start_comp, num_comp, avg_xmom, avg_ymom, avg_zmom, + cur_prim, cell_rhs, detJ_arr, dxInv, mf_m, + solverChoice.moistscal_horiz_adv_type, + solverChoice.moistscal_vert_adv_type, + l_use_terrain); +#endif if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index dbcfeb907..5d58d3082 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -114,8 +114,8 @@ void erf_slow_rhs_pre (int /*level*/, int nrk, int num_comp = 2; int end_comp = start_comp + num_comp - 1; - const int l_horiz_spatial_order = solverChoice.horiz_spatial_order; - const int l_vert_spatial_order = solverChoice.vert_spatial_order; + const int l_horiz_adv_type = solverChoice.dycore_horiz_adv_type; + const int l_vert_adv_type = solverChoice.dycore_vert_adv_type; const bool l_use_terrain = solverChoice.use_terrain; const bool l_moving_terrain = (solverChoice.terrain_type == 1); if (l_moving_terrain) AMREX_ALWAYS_ASSERT (l_use_terrain); @@ -638,7 +638,7 @@ void erf_slow_rhs_pre (int /*level*/, int nrk, avg_xmom, avg_ymom, avg_zmom, // these are being defined from the rho fluxes cell_prim, z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, - l_horiz_spatial_order, l_vert_spatial_order, l_use_terrain); + l_horiz_adv_type, l_vert_adv_type, l_use_terrain); if (l_use_diff) { Array4 diffflux_x = dflux_x->array(mfi); @@ -720,7 +720,7 @@ void erf_slow_rhs_pre (int /*level*/, int nrk, rho_u_rhs, rho_v_rhs, rho_w_rhs, u, v, w, rho_u , rho_v , omega_arr, z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v, - l_horiz_spatial_order, l_vert_spatial_order, l_use_terrain, domhi_z); + l_horiz_adv_type, l_vert_adv_type, l_use_terrain, domhi_z); if (l_use_diff) { if (l_use_terrain) { diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index c281c51b0..34c2147e1 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -104,8 +104,8 @@ // Remember this does NOT maintain dp_0 / dz = -rho_0 g, so we can no longer // discretize the vertical pressure gradient in perturbational form. - AMREX_ALWAYS_ASSERT(solverChoice.horiz_spatial_order == 2); - AMREX_ALWAYS_ASSERT(solverChoice.vert_spatial_order == 2); + AMREX_ALWAYS_ASSERT(solverChoice.dycore_horiz_adv_type == AdvType::Centered_2nd); + AMREX_ALWAYS_ASSERT(solverChoice.dycore_vert_adv_type == AdvType::Centered_2nd); Real dt_base = (new_stage_time - old_step_time); @@ -247,7 +247,7 @@ const Real old_step_time, const Real old_stage_time, const Real new_stage_time, - const int nrk) + const int /*nrk*/) { if (verbose) Print() << "Making slow rhs at time " << old_stage_time << " for slow variables advancing from " << diff --git a/Source/Utils/Interpolation.H b/Source/Utils/Interpolation.H index 0c6c6b06b..1d602e681 100644 --- a/Source/Utils/Interpolation.H +++ b/Source/Utils/Interpolation.H @@ -15,26 +15,20 @@ AMREX_FORCE_INLINE amrex::Real interpolatedVal (amrex::Real avg1, amrex::Real avg2, amrex::Real avg3, amrex::Real diff1, amrex::Real diff2, amrex::Real diff3, - amrex::Real scaled_upw, int spatial_order) + amrex::Real scaled_upw, const int adv_type) { - amrex::Real myInterpolatedVal; - switch (spatial_order) { - case 2: - myInterpolatedVal = 0.5 * avg1; - break; - case 3: - myInterpolatedVal = (7.0/12.0)*avg1 -(1.0/12.0)*avg2 + (scaled_upw/12.0)*(diff2 - 3.0*diff1); - break; - case 4: - myInterpolatedVal = (7.0/12.0)*avg1 -(1.0/12.0)*avg2; - break; - case 5: - myInterpolatedVal = (37.0/60.0)*avg1 -(2.0/15.0)*avg2 +(1.0/60.0)*avg3 - -(scaled_upw/60.0)*(diff3 - 5.0*diff2 + 10.0*diff1); - break; - case 6: - myInterpolatedVal = (37.0/60.0)*avg1 -(2.0/15.0)*avg2 +(1.0/60.0)*avg3; - break; + amrex::Real myInterpolatedVal(0.); + if (adv_type == Centered_2nd) { + myInterpolatedVal = 0.5 * avg1; + } else if (adv_type == AdvType::Upwind_3rd) { + myInterpolatedVal = (7.0/12.0)*avg1 -(1.0/12.0)*avg2 + (scaled_upw/12.0)*(diff2 - 3.0*diff1); + } else if (adv_type == AdvType::Centered_4th) { + myInterpolatedVal = (7.0/12.0)*avg1 -(1.0/12.0)*avg2; + } else if (adv_type == AdvType::Upwind_5th) { + myInterpolatedVal = (37.0/60.0)*avg1 -(2.0/15.0)*avg2 +(1.0/60.0)*avg3 + -(scaled_upw/60.0)*(diff3 - 5.0*diff2 + 10.0*diff1); + } else if (adv_type == AdvType::Centered_6th) { + myInterpolatedVal = (37.0/60.0)*avg1 -(2.0/15.0)*avg2 +(1.0/60.0)*avg3; } return myInterpolatedVal; } @@ -43,10 +37,10 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateInX (int i, int j, int k, const amrex::Array4& qty, - int qty_index, amrex::Real upw, int spatial_order) + int qty_index, amrex::Real upw, const int adv_type) { - if (spatial_order == 2) { + if (adv_type == Centered_2nd) { return 0.5 * (qty(i,j,k,qty_index) + qty(i-1,j,k,qty_index)); } else { @@ -62,12 +56,12 @@ InterpolateInX (int i, int j, int k, const amrex::Array4& qty diff1 = (qty(i, j, k, qty_index) - qty(i-1, j, k, qty_index)); avg2 = (qty(i+1, j, k, qty_index) + qty(i-2, j, k, qty_index)); diff2 = (qty(i+1, j, k, qty_index) - qty(i-2, j, k, qty_index)); - if (spatial_order > 4) + if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th) { avg3 = (qty(i+2, j, k, qty_index) + qty(i-3, j, k, qty_index)); diff3 = (qty(i+2, j, k, qty_index) - qty(i-3, j, k, qty_index)); } - return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,spatial_order); + return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type); } } @@ -75,9 +69,9 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateInY (int i, int j, int k, const amrex::Array4& qty, - int qty_index, amrex::Real upw, int spatial_order) + int qty_index, amrex::Real upw, const int adv_type) { - if (spatial_order == 2) { + if (adv_type == Centered_2nd) { return 0.5 * (qty(i,j,k,qty_index) + qty(i,j-1,k,qty_index)); } else { @@ -93,12 +87,12 @@ InterpolateInY (int i, int j, int k, const amrex::Array4& qty diff1 = (qty(i, j , k, qty_index) - qty(i, j-1, k, qty_index)); avg2 = (qty(i, j+1, k, qty_index) + qty(i, j-2, k, qty_index)); diff2 = (qty(i, j+1, k, qty_index) - qty(i, j-2, k, qty_index)); - if (spatial_order > 4) + if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th) { avg3 = (qty(i, j+2, k, qty_index) + qty(i, j-3, k, qty_index)); diff3 = (qty(i, j+2, k, qty_index) - qty(i, j-3, k, qty_index)); } - return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,spatial_order); + return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type); } } @@ -106,9 +100,9 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateInZ (int i, int j, int k, const amrex::Array4& qty, - int qty_index, amrex::Real upw, int spatial_order) + int qty_index, amrex::Real upw, const int adv_type) { - if (spatial_order == 2) { + if (adv_type == Centered_2nd) { return 0.5 * (qty(i,j,k,qty_index) + qty(i,j,k-1,qty_index)); } else { @@ -123,12 +117,12 @@ InterpolateInZ (int i, int j, int k, const amrex::Array4& qty diff1 = (qty(i, j, k , qty_index) - qty(i, j, k-1, qty_index)); avg2 = (qty(i, j, k+1, qty_index) + qty(i, j, k-2, qty_index)); diff2 = (qty(i, j, k+1, qty_index) - qty(i, j, k-2, qty_index)); - if (spatial_order > 4) + if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th) { avg3 = (qty(i, j, k+2, qty_index) + qty(i, j, k-3, qty_index)); diff3 = (qty(i, j, k+2, qty_index) - qty(i, j, k-3, qty_index)); } - return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,spatial_order); + return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type); } } @@ -137,8 +131,8 @@ AMREX_FORCE_INLINE amrex::Real InterpolatePertFromCell (int i, int j, int k, const amrex::Array4& qty, - int qty_index, amrex::Real upw, Coord coordDir, int spatial_order, - const amrex::Array4& r0_arr) + int qty_index, amrex::Real upw, Coord coordDir, + const int adv_type, const amrex::Array4& r0_arr) { amrex::Real avg1 = 0.; amrex::Real avg2 = 0.; amrex::Real avg3 = 0.; amrex::Real diff1 = 0.; amrex::Real diff2 = 0.; amrex::Real diff3 = 0.; @@ -152,13 +146,13 @@ InterpolatePertFromCell (int i, int j, int k, avg1 = (qty(i , j, k, qty_index) + qty(i-1, j, k, qty_index)); avg1 -= (r0_arr(i,j,k) + r0_arr(i-1,j,k)); diff1 = (qty(i , j, k, qty_index) - qty(i-1, j, k, qty_index)); - if (spatial_order > 2) + if (adv_type != Centered_2nd) { avg2 = (qty(i+1, j, k, qty_index) + qty(i-2, j, k, qty_index)); avg2 -= (r0_arr(i+1,j,k) + r0_arr(i-2,j,k)); diff2 = (qty(i+1, j, k, qty_index) - qty(i-2, j, k, qty_index)); } - if (spatial_order > 4) + if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th) { avg3 = (qty(i+2, j, k, qty_index) + qty(i-3, j, k, qty_index)); avg3 -= (r0_arr(i+2,j,k) + r0_arr(i-3,j,k)); @@ -168,13 +162,13 @@ InterpolatePertFromCell (int i, int j, int k, avg1 = (qty(i, j , k, qty_index) + qty(i, j-1, k, qty_index)); avg1 -= (r0_arr(i,j,k) + r0_arr(i,j-1,k)); diff1 = (qty(i, j , k, qty_index) - qty(i, j-1, k, qty_index)); - if (spatial_order > 2) + if (adv_type != Centered_2nd) { avg2 = (qty(i, j+1, k, qty_index) + qty(i, j-2, k, qty_index)); avg2 -= (r0_arr(i,j+1,k) + r0_arr(i,j-2,k)); diff2 = (qty(i, j+1, k, qty_index) - qty(i, j-2, k, qty_index)); } - if (spatial_order > 4) + if (adv_type == AdvType::Upwind_5th || adv_type == Centered_6th) { avg3 = (qty(i, j+2, k, qty_index) + qty(i, j-3, k, qty_index)); avg3 -= (r0_arr(i,j+2,k) + r0_arr(i,j-3,k)); @@ -185,14 +179,15 @@ InterpolatePertFromCell (int i, int j, int k, diff1 = (qty(i, j, k , qty_index) - qty(i, j, k-1, qty_index)); avg1 -= (r0_arr(i,j,k) + r0_arr(i,j,k-1)); diff1 -= (r0_arr(i,j,k) - r0_arr(i,j,k-1)); - if (spatial_order > 2) + + if (adv_type != Centered_2nd) { avg2 = (qty(i, j, k+1, qty_index) + qty(i, j, k-2, qty_index)); diff2 = (qty(i, j, k+1, qty_index) - qty(i, j, k-2, qty_index)); avg2 -= (r0_arr(i,j,k+1) + r0_arr(i,j,k-2)); diff2 -= (r0_arr(i,j,k+1) - r0_arr(i,j,k-2)); } - if (spatial_order > 4) + if (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th) { avg3 = (qty(i, j, k+2, qty_index) + qty(i, j, k-3, qty_index)); diff3 = (qty(i, j, k+2, qty_index) - qty(i, j, k-3, qty_index)); @@ -201,16 +196,16 @@ InterpolatePertFromCell (int i, int j, int k, } } - return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,spatial_order); + return interpolatedVal(avg1,avg2,avg3,diff1,diff2,diff3,scaled_upw,adv_type); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpolateDensityPertFromCellToFace (int i, int j, int k, const amrex::Array4& cons_in, - amrex::Real upw, Coord coordDir, int spatial_order, + amrex::Real upw, Coord coordDir, const int adv_type, const amrex::Array4& r0_arr) { - return InterpolatePertFromCell(i, j, k, cons_in, Rho_comp, upw, coordDir, spatial_order, r0_arr); + return InterpolatePertFromCell(i, j, k, cons_in, Rho_comp, upw, coordDir, adv_type, r0_arr); } #endif diff --git a/Source/Utils/Interpolation_UPW.H b/Source/Utils/Interpolation_UPW.H index 69c8982eb..6a1d20d32 100644 --- a/Source/Utils/Interpolation_UPW.H +++ b/Source/Utils/Interpolation_UPW.H @@ -4,11 +4,11 @@ #include "DataStruct.H" /** - * Interpolation operators used for 2nd order central scheme + * Interpolation operators used for 2nd order centered scheme */ -struct UPWIND2 +struct CENTERED2 { - UPWIND2(const amrex::Array4& phi) + CENTERED2(const amrex::Array4& phi) : m_phi(phi) {} AMREX_GPU_DEVICE @@ -254,11 +254,11 @@ private: /** - * Interpolation operators used for 4th order central scheme + * Interpolation operators used for 4th order centered scheme */ -struct UPWIND4 +struct CENTERED4 { - UPWIND4(const amrex::Array4& phi) + CENTERED4(const amrex::Array4& phi) : m_phi(phi) {} AMREX_GPU_DEVICE @@ -528,11 +528,11 @@ private: }; /** - * Interpolation operators used for 6th order central scheme + * Interpolation operators used for 6th order centered scheme */ -struct UPWIND6 +struct CENTERED6 { - UPWIND6(const amrex::Array4& phi) + CENTERED6(const amrex::Array4& phi) : m_phi(phi) {} AMREX_GPU_DEVICE @@ -678,25 +678,25 @@ struct UPWINDALL const int& qty_index, amrex::Real& val_lo, amrex::Real upw_lo, - const int& spatial_order) const + const int adv_type) const { - if (spatial_order == 2) { + if (adv_type == AdvType::Centered_2nd) { val_lo = 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 sp2 = (spatial_order > 4) ? m_phi(i , j , k+2, qty_index): 0.; + amrex::Real sp2 = (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th ) ? m_phi(i , j , k+2, qty_index): 0.; 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); - amrex::Real sm3 = (spatial_order > 4) ? m_phi(i , j , k-3, qty_index) : 0.; + amrex::Real sm3 = (adv_type == AdvType::Upwind_5th || adv_type == AdvType::Centered_6th ) ? m_phi(i , j , k-3, qty_index) : 0.; // Upwinding flags if (upw_lo != 0.) upw_lo = (upw_lo > 0) ? 1. : -1.; // Interpolate lo - val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3,upw_lo,spatial_order); + val_lo = Evaluate(sp2,sp1,s,sm1,sm2,sm3,upw_lo,adv_type); } } @@ -709,25 +709,25 @@ struct UPWINDALL const int& qty_index, amrex::Real& val_hi, amrex::Real upw_hi, - const int& spatial_order) const + const int adv_type) const { - if (spatial_order == 2) { + 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 = (spatial_order > 4) ? m_phi(i , j , k+3, qty_index): 0.; + 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 = (spatial_order > 4) ? m_phi(i , j , k-2, qty_index) : 0.; + 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,spatial_order); + val_hi = Evaluate(sp3,sp2,sp1,s,sm1,sm2,upw_hi,adv_type); } } @@ -741,7 +741,7 @@ struct UPWINDALL const amrex::Real& sm2, const amrex::Real& sm3, const amrex::Real& upw, - const int& spatial_order) const + const int adv_type) const { // Averages and diffs amrex::Real a1 = (s + sm1); @@ -752,23 +752,17 @@ struct UPWINDALL amrex::Real d3 = (sp2 - sm3); // Interpolated value - amrex::Real iv; - switch (spatial_order) { - case 2: + amrex::Real iv(0.0); + if (adv_type == AdvType::Centered_2nd) { iv = 0.5 * a1; - break; - case 3: + } else if (adv_type == AdvType::Upwind_3rd) { iv = g1_3_4*a1 - g2_3_4*a2 + upw * g2_3_4 * (d2 - 3.0*d1); - break; - case 4: + } else if (adv_type == AdvType::Centered_4th) { iv = g1_3_4*a1 - g2_3_4*a2; - break; - case 5: + } else if (adv_type == AdvType::Upwind_5th) { iv = g1_5_6*a1 - g2_5_6*a2 + g3_5_6*a3 - upw * g3_5_6 * (d3 - 5.0*d2 + 10.0*d1); - break; - case 6: + } else if (adv_type == AdvType::Centered_6th) { iv = g1_5_6*a1 - g2_5_6*a2 + g3_5_6*a3; - break; } return ( iv ); } diff --git a/Source/Utils/Interpolation_WENO_Z.H b/Source/Utils/Interpolation_WENO_Z.H index c3e9ff240..216437d76 100644 --- a/Source/Utils/Interpolation_WENO_Z.H +++ b/Source/Utils/Interpolation_WENO_Z.H @@ -168,6 +168,171 @@ private: static constexpr amrex::Real g2=(2.0/3.0); }; +/** + * Interpolation operators used for WENO_MZQ3 scheme + */ +struct WENO_MZQ3 +{ + WENO_MZQ3(const amrex::Array4& phi) + : m_phi(phi) {} + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void + InterpolateInX(const int& i, + 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) { + val_lo = Evaluate(sp1,s ,sm1); + } else { + val_lo = 0.5 * (s + sm1); + } + } + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void + InterpolateInY(const int& i, + 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) { + val_lo = Evaluate(sp1,s ,sm1); + } else { + val_lo = 0.5 * (s + sm1); + } + } + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void + InterpolateInZ_lo(const int& i, + const int& j, + const int& k, + const int& qty_index, + amrex::Real& val_lo, + amrex::Real upw_lo) 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); + amrex::Real sm1 = m_phi(i , j , k-1, qty_index); + amrex::Real sm2 = m_phi(i , j , k-2, qty_index); + + if (upw_lo > tol) { + val_lo = Evaluate(sm2,sm1,s ); + } else if (upw_lo < -tol) { + val_lo = Evaluate(sp1,s ,sm1); + } else { + val_lo = 0.5 * (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 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 + Evaluate(const amrex::Real& sm1, + const amrex::Real& s , + const amrex::Real& sp1) const + { + // Smoothing factors + amrex::Real b1 = (s - sm1) * (s - sm1); + amrex::Real b2 = (sp1 - s) * (sp1 - s); + + // Weight factors + amrex::Real t5 = std::abs(b2 - b1); + amrex::Real w1 = g1 * ( 1.0 + (t5*t5) / ((eps + b1) * (eps + b1)) ); + amrex::Real w2 = g2 * ( 1.0 + (t5*t5) / ((eps + b2) * (eps + b2)) ); + + // Weight factor norm + amrex::Real wsum = w1 + w2; + + // Taylor expansions + amrex::Real v1 = -sm1 + 3.0 * s; + amrex::Real v2 = s + sp1; + + // Interpolated value + return ( (w1 * v1 + w2 * v2) / (2.0 * wsum) ); + } + +private: + amrex::Array4 m_phi; // Quantity to interpolate + const amrex::Real eps=1.0e-6; + const amrex::Real tol=1.0e-12; + static constexpr amrex::Real g1=(1.0/3.0); + static constexpr amrex::Real g2=(2.0/3.0); +}; + /** * Interpolation operators used for WENO_Z-5 scheme */ diff --git a/Tests/test_files/Bubble_DensityCurrent/Bubble_DensityCurrent.i b/Tests/test_files/Bubble_DensityCurrent/Bubble_DensityCurrent.i index 492cefbf7..80f7d1b9d 100644 --- a/Tests/test_files/Bubble_DensityCurrent/Bubble_DensityCurrent.i +++ b/Tests/test_files/Bubble_DensityCurrent/Bubble_DensityCurrent.i @@ -46,8 +46,6 @@ erf.alpha_C = 0.0 erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" diff --git a/Tests/test_files/CouetteFlow/CouetteFlow.i b/Tests/test_files/CouetteFlow/CouetteFlow.i index a13fc1399..a77b035f4 100644 --- a/Tests/test_files/CouetteFlow/CouetteFlow.i +++ b/Tests/test_files/CouetteFlow/CouetteFlow.i @@ -48,9 +48,6 @@ erf.les_type = "None" erf.molec_diff_type = "Constant" erf.dynamicViscosity = 0.1 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 300.0 diff --git a/Tests/test_files/Deardorff_stationary/Deardorff_stationary.i b/Tests/test_files/Deardorff_stationary/Deardorff_stationary.i index 232baa5aa..2b49b6355 100644 --- a/Tests/test_files/Deardorff_stationary/Deardorff_stationary.i +++ b/Tests/test_files/Deardorff_stationary/Deardorff_stationary.i @@ -48,7 +48,6 @@ erf.plot_int_1 = 10 # number of timesteps between plotfiles (DEBUG) erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoKE #pres_hse dens_hse # SOLVER CHOICES -erf.spatial_order = 2 erf.use_gravity = false erf.use_coriolis = false erf.use_rayleigh_damping = false diff --git a/Tests/test_files/DensityCurrent/DensityCurrent.i b/Tests/test_files/DensityCurrent/DensityCurrent.i index b38b8694a..545903e3f 100644 --- a/Tests/test_files/DensityCurrent/DensityCurrent.i +++ b/Tests/test_files/DensityCurrent/DensityCurrent.i @@ -48,8 +48,6 @@ erf.alpha_C = 0.0 erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" diff --git a/Tests/test_files/DensityCurrent_detJ2/DensityCurrent_detJ2.i b/Tests/test_files/DensityCurrent_detJ2/DensityCurrent_detJ2.i index ef7280819..e3feceeba 100644 --- a/Tests/test_files/DensityCurrent_detJ2/DensityCurrent_detJ2.i +++ b/Tests/test_files/DensityCurrent_detJ2/DensityCurrent_detJ2.i @@ -48,8 +48,6 @@ erf.alpha_C = 0.0 erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" 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 483374dd3..1382e7ba3 100644 --- a/Tests/test_files/DensityCurrent_detJ2_MT/DensityCurrent_detJ2_MT.i +++ b/Tests/test_files/DensityCurrent_detJ2_MT/DensityCurrent_detJ2_MT.i @@ -48,8 +48,6 @@ erf.alpha_C = 0.0 erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" diff --git a/Tests/test_files/DensityCurrent_detJ2_nosub/DensityCurrent_detJ2_nosub.i b/Tests/test_files/DensityCurrent_detJ2_nosub/DensityCurrent_detJ2_nosub.i index 466fdb58a..0d5a6f719 100644 --- a/Tests/test_files/DensityCurrent_detJ2_nosub/DensityCurrent_detJ2_nosub.i +++ b/Tests/test_files/DensityCurrent_detJ2_nosub/DensityCurrent_detJ2_nosub.i @@ -48,9 +48,6 @@ erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - erf.les_type = "None" erf.molec_diff_type = "ConstantAlpha" # diffusion = 75 m^2/s, rho_0 = 1e5/(287*300) = 1.1614401858 diff --git a/Tests/test_files/EkmanSpiral/EkmanSpiral.i b/Tests/test_files/EkmanSpiral/EkmanSpiral.i index 9d5e18693..4bd46a90e 100644 --- a/Tests/test_files/EkmanSpiral/EkmanSpiral.i +++ b/Tests/test_files/EkmanSpiral/EkmanSpiral.i @@ -47,9 +47,6 @@ erf.les_type = "None" erf.molec_diff_type = "Constant" erf.dynamicViscosity = 5.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - erf.use_coriolis = true erf.abl_driver_type = "GeostrophicWind" erf.latitude = 90. diff --git a/Tests/test_files/IsentropicVortexAdvecting/IsentropicVortexAdvecting.i b/Tests/test_files/IsentropicVortexAdvecting/IsentropicVortexAdvecting.i index ddf7a1c2c..5bc7bc6b1 100644 --- a/Tests/test_files/IsentropicVortexAdvecting/IsentropicVortexAdvecting.i +++ b/Tests/test_files/IsentropicVortexAdvecting/IsentropicVortexAdvecting.i @@ -45,9 +45,6 @@ erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.p_inf = 1e5 # reference pressure [Pa] prob.T_inf = 300. # reference temperature [K] diff --git a/Tests/test_files/IsentropicVortexStationary/IsentropicVortexStationary.i b/Tests/test_files/IsentropicVortexStationary/IsentropicVortexStationary.i index 95edd2994..e9b5eac69 100644 --- a/Tests/test_files/IsentropicVortexStationary/IsentropicVortexStationary.i +++ b/Tests/test_files/IsentropicVortexStationary/IsentropicVortexStationary.i @@ -41,9 +41,6 @@ erf.alpha_T = 0.0 erf.alpha_C = 0.0 erf.use_gravity = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 diff --git a/Tests/test_files/MSF_NoSub_IsentropicVortexAdv/MSF_NoSub_IsentropicVortexAdv.i b/Tests/test_files/MSF_NoSub_IsentropicVortexAdv/MSF_NoSub_IsentropicVortexAdv.i index 625922dc1..e16035dfa 100644 --- a/Tests/test_files/MSF_NoSub_IsentropicVortexAdv/MSF_NoSub_IsentropicVortexAdv.i +++ b/Tests/test_files/MSF_NoSub_IsentropicVortexAdv/MSF_NoSub_IsentropicVortexAdv.i @@ -48,9 +48,6 @@ erf.les_type = "None" erf.molec_diff_type = "Constant" erf.dynamicViscosity = 1.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.p_inf = 1e5 # reference pressure [Pa] prob.T_inf = 300. # reference temperature [K] diff --git a/Tests/test_files/MSF_Sub_IsentropicVortexAdv/MSF_Sub_IsentropicVortexAdv.i b/Tests/test_files/MSF_Sub_IsentropicVortexAdv/MSF_Sub_IsentropicVortexAdv.i index 43d45f959..d8bcea2d7 100644 --- a/Tests/test_files/MSF_Sub_IsentropicVortexAdv/MSF_Sub_IsentropicVortexAdv.i +++ b/Tests/test_files/MSF_Sub_IsentropicVortexAdv/MSF_Sub_IsentropicVortexAdv.i @@ -48,9 +48,6 @@ erf.les_type = "None" erf.molec_diff_type = "Constant" erf.dynamicViscosity = 1.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.p_inf = 1e5 # reference pressure [Pa] prob.T_inf = 300. # reference temperature [K] diff --git a/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i b/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i index 989990f2b..be3d153a2 100644 --- a/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i +++ b/Tests/test_files/MovingTerrain_nosub/MovingTerrain_nosub.i @@ -46,9 +46,6 @@ erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - erf.buoyancy_type = 1 # MULTILEVEL diff --git a/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i b/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i index 7576830e1..ef688d2ce 100644 --- a/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i +++ b/Tests/test_files/MovingTerrain_sub/MovingTerrain_sub.i @@ -49,9 +49,6 @@ erf.use_gravity = true erf.use_coriolis = false erf.use_rayleigh_damping = false -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - erf.buoyancy_type = 1 # MULTILEVEL diff --git a/Tests/test_files/PoiseuilleFlow/PoiseuilleFlow.i b/Tests/test_files/PoiseuilleFlow/PoiseuilleFlow.i index 262218554..ffb9d33c8 100644 --- a/Tests/test_files/PoiseuilleFlow/PoiseuilleFlow.i +++ b/Tests/test_files/PoiseuilleFlow/PoiseuilleFlow.i @@ -52,9 +52,6 @@ erf.use_coriolis = false erf.abl_driver_type = "PressureGradient" erf.abl_pressure_grad = -0.2 0. 0. -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 300.0 diff --git a/Tests/test_files/RayleighDamping/RayleighDamping.i b/Tests/test_files/RayleighDamping/RayleighDamping.i index b01361a59..ddbc84636 100644 --- a/Tests/test_files/RayleighDamping/RayleighDamping.i +++ b/Tests/test_files/RayleighDamping/RayleighDamping.i @@ -52,9 +52,6 @@ erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvDiff_order2/ScalarAdvDiff_order2.i b/Tests/test_files/ScalarAdvDiff_order2/ScalarAdvDiff_order2.i index 0b8a14b16..38616dbcc 100644 --- a/Tests/test_files/ScalarAdvDiff_order2/ScalarAdvDiff_order2.i +++ b/Tests/test_files/ScalarAdvDiff_order2/ScalarAdvDiff_order2.i @@ -53,9 +53,6 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.A_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvDiff_order3/ScalarAdvDiff_order3.i b/Tests/test_files/ScalarAdvDiff_order3/ScalarAdvDiff_order3.i index bcdf4adfb..a3b626e8d 100644 --- a/Tests/test_files/ScalarAdvDiff_order3/ScalarAdvDiff_order3.i +++ b/Tests/test_files/ScalarAdvDiff_order3/ScalarAdvDiff_order3.i @@ -53,8 +53,10 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 3 -erf.vert_spatial_order = 3 +erf.dycore_horiz_adv_type = Upwind_3rd +erf.dycore_vert_adv_type = Upwind_3rd +erf.dryscal_horiz_adv_type = Upwind_3rd +erf.dryscal_vert_adv_type = Upwind_3rd # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvDiff_order4/ScalarAdvDiff_order4.i b/Tests/test_files/ScalarAdvDiff_order4/ScalarAdvDiff_order4.i index 4913d6345..d84995afd 100644 --- a/Tests/test_files/ScalarAdvDiff_order4/ScalarAdvDiff_order4.i +++ b/Tests/test_files/ScalarAdvDiff_order4/ScalarAdvDiff_order4.i @@ -53,8 +53,10 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 4 -erf.vert_spatial_order = 4 +erf.dycore_horiz_adv_type = Centered_4th +erf.dycore_vert_adv_type = Centered_4th +erf.dryscal_horiz_adv_type = Centered_4th +erf.dryscal_vert_adv_type = Centered_4th # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvDiff_order5/ScalarAdvDiff_order5.i b/Tests/test_files/ScalarAdvDiff_order5/ScalarAdvDiff_order5.i index eb59f9b2f..dd29b60dc 100644 --- a/Tests/test_files/ScalarAdvDiff_order5/ScalarAdvDiff_order5.i +++ b/Tests/test_files/ScalarAdvDiff_order5/ScalarAdvDiff_order5.i @@ -53,8 +53,10 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 5 -erf.vert_spatial_order = 5 +erf.dycore_horiz_adv_type = Upwind_5th +erf.dycore_vert_adv_type = Upwind_5th +erf.dryscal_horiz_adv_type = Upwind_5th +erf.dryscal_vert_adv_type = Upwind_5th # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvDiff_order6/ScalarAdvDiff_order6.i b/Tests/test_files/ScalarAdvDiff_order6/ScalarAdvDiff_order6.i index 37e7adc0b..46117bc14 100644 --- a/Tests/test_files/ScalarAdvDiff_order6/ScalarAdvDiff_order6.i +++ b/Tests/test_files/ScalarAdvDiff_order6/ScalarAdvDiff_order6.i @@ -53,8 +53,10 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 6 -erf.vert_spatial_order = 6 +erf.dycore_horiz_adv_type = Centered_6th +erf.dycore_vert_adv_type = Centered_6th +erf.dryscal_horiz_adv_type = Centered_6th +erf.dryscal_vert_adv_type = Centered_6th # PROBLEM PARAMETERS prob.rho_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvectionShearedU/ScalarAdvectionShearedU.i b/Tests/test_files/ScalarAdvectionShearedU/ScalarAdvectionShearedU.i index 59b1567a9..1ae466f7d 100644 --- a/Tests/test_files/ScalarAdvectionShearedU/ScalarAdvectionShearedU.i +++ b/Tests/test_files/ScalarAdvectionShearedU/ScalarAdvectionShearedU.i @@ -45,9 +45,6 @@ erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 diff --git a/Tests/test_files/ScalarAdvectionUniformU/ScalarAdvectionUniformU.i b/Tests/test_files/ScalarAdvectionUniformU/ScalarAdvectionUniformU.i index 8f99d0966..e28cc1064 100644 --- a/Tests/test_files/ScalarAdvectionUniformU/ScalarAdvectionUniformU.i +++ b/Tests/test_files/ScalarAdvectionUniformU/ScalarAdvectionUniformU.i @@ -41,13 +41,15 @@ erf.alpha_T = 0.0 erf.alpha_C = 0.0 erf.use_gravity = false +erf.dryscal_horiz_adv_type = "Centered_2nd" +erf.dryscal_vert_adv_type = "Centered_2nd" +erf.moistscal_horiz_adv_type = "Centered_2nd" +erf.moistscal_vert_adv_type = "Centered_2nd" + erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 diff --git a/Tests/test_files/ScalarDiffusionGaussian/ScalarDiffusionGaussian.i b/Tests/test_files/ScalarDiffusionGaussian/ScalarDiffusionGaussian.i index a1f0d2c52..d80b3aa9b 100644 --- a/Tests/test_files/ScalarDiffusionGaussian/ScalarDiffusionGaussian.i +++ b/Tests/test_files/ScalarDiffusionGaussian/ScalarDiffusionGaussian.i @@ -45,9 +45,6 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.A_0 = 1.0 diff --git a/Tests/test_files/ScalarDiffusionSine/ScalarDiffusionSine.i b/Tests/test_files/ScalarDiffusionSine/ScalarDiffusionSine.i index 018e553ca..ba7d729c5 100644 --- a/Tests/test_files/ScalarDiffusionSine/ScalarDiffusionSine.i +++ b/Tests/test_files/ScalarDiffusionSine/ScalarDiffusionSine.i @@ -45,9 +45,6 @@ erf.molec_diff_type = "Constant" erf.rho0_trans = 1.0 erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.A_0 = 0.0 diff --git a/Tests/test_files/TaylorGreenAdvecting/TaylorGreenAdvecting.i b/Tests/test_files/TaylorGreenAdvecting/TaylorGreenAdvecting.i index c79f83628..49ba19cff 100644 --- a/Tests/test_files/TaylorGreenAdvecting/TaylorGreenAdvecting.i +++ b/Tests/test_files/TaylorGreenAdvecting/TaylorGreenAdvecting.i @@ -44,9 +44,6 @@ erf.les_type = "None" erf.molec_diff_type = "None" erf.dynamicViscosity = 0.0 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.A_0 = 1.0 diff --git a/Tests/test_files/TaylorGreenAdvectingDiffusing/TaylorGreenAdvectingDiffusing.i b/Tests/test_files/TaylorGreenAdvectingDiffusing/TaylorGreenAdvectingDiffusing.i index 13bf474c3..fc44afc6b 100644 --- a/Tests/test_files/TaylorGreenAdvectingDiffusing/TaylorGreenAdvectingDiffusing.i +++ b/Tests/test_files/TaylorGreenAdvectingDiffusing/TaylorGreenAdvectingDiffusing.i @@ -43,9 +43,6 @@ erf.les_type = "None" erf.molec_diff_type = "Constant" erf.dynamicViscosity = 6.25e-4 # 1.5e-5 -erf.horiz_spatial_order = 2 -erf.vert_spatial_order = 2 - # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.A_0 = 1.0