diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 838c64a6a..e6ce285d1 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -456,76 +456,86 @@ which means it is determined by the fluid speed rather than the sound speed and List of Parameters ------------------ -+----------------------------+----------------------+----------------+---------------------+ -| Parameter | Definition | Acceptable | Default | -| | | Values | | -+============================+======================+================+=====================+ -| **erf.substepping_type** | Should we substep in | "Implicit" or | "Implicit" if | -| | each RK stage? | "None" | compressible, | -| | | | "None" if anelastic | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.vert_implicit** | Do vertical implicit | Boolean | false | -| | solve for diffusion | | | -| | of u, v, and theta | | | -| | with default | | | -| | time-centering in | | | -| | each stage | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.vert_implicit_fac** | How much implicit | Real >= 0 | 0.0, 0.0, 0.0 | -| | vertical diffusion | (explicit) and | (fully explicit) | -| | to include in each | <= 1 (implicit)| | -| | RK stage? Currently, | | | -| | only applies to | | | -| | rho*theta component. | | | -| | | | | -| | Specify either one | | | -| | (the same for all | | | -| | stages) or three | | | -| | values, one per stage| | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.cfl** | CFL number used to | Real > 0 and | 0.8 | -| | compute level 0 dt | <= 1 | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.substepping_cfl** | CFL number used to | Real > 0 and | 1.0 | -| | compute the number | <= 1 | | -| | of substeps | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.fixed_dt** | set level 0 dt | Real > 0 | unused if not | -| | as this value | | set | -| | regardless of | | | -| | cfl or other | | | -| | settings | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.fixed_fast_dt** | set fast dt | Real > 0 | | -| | as this value | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.fixed_mri_dt_ratio** | set fast dt | even int > 0 | only relevant if | -| | as slow dt / | | substepping_type | -| | this ratio | | is not None | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.init_shrink** | factor by which | Real > 0 and | 1.0 | -| | to shrink the | <= 1 | | -| | initial dt | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.change_max** | factor by which | Real >= 1 | 1.1 | -| | dt can grow | | | -| | in subsequent | | | -| | steps | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.dt_max** | maximum adaptive | Real > 0 | 1e9 | -| | timestep | | | -| | allowed by time | | | -| | stepping | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.dt_max_initial** | maximum initial | Real > 0 | 1.0 | -| | timestep | | | -+----------------------------+----------------------+----------------+---------------------+ -| **erf.dt_ref_ratio** | ratio of coarse | Integer >= 1 | same as | -| | to fine grid | (one per level)| maximum over | -| | time steps between | | directions of | -| | subsequent | | ref_ratio | -| | levels | | | -+----------------------------+----------------------+----------------+---------------------+ ++-------------------------------------+----------------------+----------------+---------------------+ +| Parameter | Definition | Acceptable | Default | +| | | Values | | ++=====================================+======================+================+=====================+ +| **erf.substepping_type** | Should we substep in | "Implicit" or | "Implicit" if | +| | each RK stage? | "None" | compressible, | +| | | | "None" if anelastic | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.vert_implicit** | Do vertical implicit | Boolean | false | +| | solve for diffusion | | | +| | of u, v, theta, KE, | | | +| | and qv with default | | | +| | time-centering in | | | +| | each stage | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.vert_implicit_fac** | How much implicit | Real >= 0 | 0.0, 0.0, 0.0 | +| | vertical diffusion | (explicit) and | (fully explicit) | +| | to include in each | <= 1 (implicit)| | +| | RK stage? | | | +| | | | | +| | Specify either one | | | +| | (the same for all | | | +| | stages) or three | | | +| | values, one per stage| | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.implicit_thermal_diffusion** | Do vertical implicit | Boolean | true | +| | for theta? | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.implicit_moisture_diffusion** | Do vertical implicit | Boolean | true | +| | for Qv? | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.implicit_ke_diffusion** | Do vertical implicit | Boolean | true | +| | for KE? | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.implicit_momentum_diffusion** | Do vertical implicit | Boolean | true | +| | for U & V? | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.cfl** | CFL number used to | Real > 0 and | 0.8 | +| | compute level 0 dt | <= 1 | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.substepping_cfl** | CFL number used to | Real > 0 and | 1.0 | +| | compute the number | <= 1 | | +| | of substeps | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.fixed_dt** | set level 0 dt | Real > 0 | unused if not | +| | as this value | | set | +| | regardless of | | | +| | cfl or other | | | +| | settings | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.fixed_fast_dt** | set fast dt | Real > 0 | | +| | as this value | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.fixed_mri_dt_ratio** | set fast dt | even int > 0 | only relevant if | +| | as slow dt / | | substepping_type | +| | this ratio | | is not None | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.init_shrink** | factor by which | Real > 0 and | 1.0 | +| | to shrink the | <= 1 | | +| | initial dt | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.change_max** | factor by which | Real >= 1 | 1.1 | +| | dt can grow | | | +| | in subsequent | | | +| | steps | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.dt_max** | maximum adaptive | Real > 0 | 1e9 | +| | timestep | | | +| | allowed by time | | | +| | stepping | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.dt_max_initial** | maximum initial | Real > 0 | 1.0 | +| | timestep | | | ++-------------------------------------+----------------------+----------------+---------------------+ +| **erf.dt_ref_ratio** | ratio of coarse | Integer >= 1 | same as | +| | to fine grid | (one per level)| maximum over | +| | time steps between | | directions of | +| | subsequent | | ref_ratio | +| | levels | | | ++-------------------------------------+----------------------+----------------+---------------------+ Notes ----------------- diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 45a7ebfb9..a24d5c861 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -682,14 +682,19 @@ struct SolverChoice { // moisture diffusion pp.query("implicit_moisture_diffusion", implicit_moisture_diffusion); + // If true (default), include implicit contributions to vertical + // KE diffusion + pp.query("implicit_ke_diffusion", implicit_ke_diffusion); + // If true (default), include implicit contributions in tau13, tau23, // (and if ERF_IMPLICIT_W is set, tau33) to correct u, v, (and w). pp.query("implicit_momentum_diffusion", implicit_momentum_diffusion); - if (!implicit_thermal_diffusion && + if (!implicit_thermal_diffusion && !implicit_moisture_diffusion && + !implicit_ke_diffusion && !implicit_momentum_diffusion) { - amrex::Print() << "Thermal and momentum diffusion are both turned off -- turning off vertical implicit solve" << std::endl; + amrex::Print() << "Thermal, moisture, KE, and momentum diffusion are both turned off -- turning off vertical implicit solve" << std::endl; vert_implicit_fac[0] = zero; vert_implicit_fac[1] = zero; vert_implicit_fac[2] = zero; @@ -1170,6 +1175,7 @@ struct SolverChoice { // if any vert_implicit_fac > 0, then the following apply: bool implicit_thermal_diffusion = true; bool implicit_moisture_diffusion = true; + bool implicit_ke_diffusion = true; bool implicit_momentum_diffusion = true; bool implicit_before_substep = true; diff --git a/Source/Diffusion/ERF_ImplicitDiff_N.cpp b/Source/Diffusion/ERF_ImplicitDiff_N.cpp index af3421ab6..a55defb06 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_N.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_N.cpp @@ -36,7 +36,7 @@ ImplicitDiffForStateLU_N (const Box& bx, const GpuArray& cellSizeInv, const Array4& scalar_zflux, const Array4& mu_turb, - const SolverChoice &solverChoice, + const SolverChoice& solverChoice, const BCRec* bc_ptr, const bool use_SurfLayer, const Real implicit_fac) @@ -129,7 +129,7 @@ ImplicitDiffForStateLU_N (const Box& bx, a_tmp = -Fact * rhoAlpha_lo * dz_inv; c_tmp = -Fact * rhoAlpha_hi * dz_inv; b_tmp = cell_data(i,j,k,Rho_comp) - a_tmp - c_tmp; - inv_b2_tmp = one/ (b_tmp - a_tmp * coeffG_a(i,j,k-1)); + inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,k-1)); RHS_a(i,j,k) = cell_data(i,j,k,n); // NOTE: this is rho*phi; solution is phi @@ -147,7 +147,7 @@ ImplicitDiffForStateLU_N (const Box& bx, a_tmp = -Fact * rhoAlpha_lo * dz_inv; c_tmp = zero; b_tmp = cell_data(i,j,khi,Rho_comp) - a_tmp - c_tmp; - inv_b2_tmp = one/ (b_tmp - a_tmp * coeffG_a(i,j,khi-1)); + inv_b2_tmp = one / (b_tmp - a_tmp * coeffG_a(i,j,khi-1)); RHS_a(i,j,khi) = cell_data(i,j,khi,n); // NOTE: this is rho*phi; solution is phi if (neumann_on_zhi) { diff --git a/Source/Diffusion/ERF_ImplicitDiff_S.cpp b/Source/Diffusion/ERF_ImplicitDiff_S.cpp index 6fbad2f2c..1c62ad638 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_S.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_S.cpp @@ -34,7 +34,7 @@ ImplicitDiffForStateLU_S (const Box& bx, const Gpu::DeviceVector& stretched_dz_d, const Array4& scalar_zflux, const Array4& mu_turb, - const SolverChoice &solverChoice, + const SolverChoice& solverChoice, const BCRec* bc_ptr, const bool use_SurfLayer, const Real implicit_fac) diff --git a/Source/Diffusion/ERF_ImplicitDiff_T.cpp b/Source/Diffusion/ERF_ImplicitDiff_T.cpp index 16c937ad3..876db775b 100644 --- a/Source/Diffusion/ERF_ImplicitDiff_T.cpp +++ b/Source/Diffusion/ERF_ImplicitDiff_T.cpp @@ -38,7 +38,7 @@ ImplicitDiffForStateLU_T (const Box& bx, const GpuArray& cellSizeInv, const Array4& scalar_zflux, const Array4& mu_turb, - const SolverChoice &solverChoice, + const SolverChoice& solverChoice, const BCRec* bc_ptr, const bool use_SurfLayer, const Real implicit_fac) diff --git a/Source/IO/ERF_ReadFromWRFBdy.cpp b/Source/IO/ERF_ReadFromWRFBdy.cpp index 84c64ba11..848a8fc1d 100644 --- a/Source/IO/ERF_ReadFromWRFBdy.cpp +++ b/Source/IO/ERF_ReadFromWRFBdy.cpp @@ -163,7 +163,6 @@ convert_wrfbdy_data (const int itime, int ihi = domain.bigEnd()[0]; int jlo = domain.smallEnd()[1]; int jhi = domain.bigEnd()[1]; - int klo = domain.smallEnd()[2]; int khi = domain.bigEnd()[2]; // PH bounds limiting diff --git a/Source/TimeIntegration/ERF_ImplicitPost.H b/Source/TimeIntegration/ERF_ImplicitPost.H index d700060fe..43a35bd1c 100644 --- a/Source/TimeIntegration/ERF_ImplicitPost.H +++ b/Source/TimeIntegration/ERF_ImplicitPost.H @@ -1,55 +1,102 @@ // ***************************************************************************** // Do semi-implicit solve for diffusion: q1 // ***************************************************************************** + const bool l_do_implicit_ke = solverChoice.implicit_ke_diffusion; const bool l_do_implicit_moist = solverChoice.implicit_moisture_diffusion; const Real l_vert_implicit_fac = solverChoice.vert_implicit_fac[nrk]; - if ( (l_do_implicit_moist) && - (l_vert_implicit_fac > zero) && - (solverChoice.moisture_type != MoistureType::None) ) { + if ( l_vert_implicit_fac > zero ) { - // BCs - const BCRec* bc_ptr_h = domain_bcs_type.data(); - GpuArray l_bc_neumann_vals_d; - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { - l_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoQ1_comp][ori]; - } - const bool l_use_SurfLayer = (m_SurfaceLayer != nullptr); + const bool l_use_turb = solverChoice.turbChoice[level].use_kturb; + const bool l_use_stretched_dz = (solverChoice.mesh_type == MeshType::StretchedDz); - const bool l_use_turb = solverChoice.turbChoice[level].use_kturb; + if ( (l_do_implicit_ke) && + (solverChoice.turbChoice[level].use_tke)) { - const bool l_use_stretched_dz = (solverChoice.mesh_type == MeshType::StretchedDz); + // BCs for q1 + const BCRec* bc_ptr_h = domain_bcs_type.data(); + GpuArray l_bc_neumann_vals_d; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoKE_comp][ori]; + } + const bool l_use_SurfLayer = false; - for ( MFIter mfi(S_new[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) - { - Box bx = mfi.tilebox(); + for ( MFIter mfi(S_new[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) + { + Box bx = mfi.tilebox(); - // NOTE: We have updated the state, use that directly - const Array4< Real>& cell_data = S_new[IntVars::cons].array(mfi); + // NOTE: We have updated the state, use that directly + const Array4< Real>& cell_data = S_new[IntVars::cons].array(mfi); - const Array4& z_nd_arr = z_phys_nd[level]->const_array(mfi); - const Array4& detJ_arr = detJ_cc[level]->const_array(mfi); + const Array4& z_nd_arr = z_phys_nd[level]->const_array(mfi); + const Array4& detJ_arr = detJ_cc[level]->const_array(mfi); - const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : + const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4{}; - const Array4& q1fx_z = Q1fx3->const_array(mfi); - - if (l_use_stretched_dz) { - ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoQ1_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, - stretched_dz_d[level], q1fx_z, - mu_turb, solverChoice, - bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); - } else if (l_use_terrain_fitted_coords) { - ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoQ1_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, - z_nd_arr, detJ_arr, dxInv, q1fx_z, - mu_turb, solverChoice, - bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); - } else { // no terrain - ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoQ1_comp, - slow_dt, l_bc_neumann_vals_d, cell_data, - dxInv, q1fx_z, mu_turb, solverChoice, - bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + + if (l_use_stretched_dz) { + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoKE_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + stretched_dz_d[level], Array4{}, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } else if (l_use_terrain_fitted_coords) { + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoKE_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + z_nd_arr, detJ_arr, dxInv, Array4{}, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } else { // no terrain + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoKE_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + dxInv, Array4{}, mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } + } // mfi + } // do implicit tke and evolve tke + + if ( (l_do_implicit_moist) && + (solverChoice.moisture_type != MoistureType::None) ) { + + // BCs for q1 + const BCRec* bc_ptr_h = domain_bcs_type.data(); + GpuArray l_bc_neumann_vals_d; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoQ1_comp][ori]; } - } // mfi - } // if moist model and do implicit diff + const bool l_use_SurfLayer = (m_SurfaceLayer != nullptr); + + for ( MFIter mfi(S_new[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi) + { + Box bx = mfi.tilebox(); + + // NOTE: We have updated the state, use that directly + const Array4< Real>& cell_data = S_new[IntVars::cons].array(mfi); + + const Array4& z_nd_arr = z_phys_nd[level]->const_array(mfi); + const Array4& detJ_arr = detJ_cc[level]->const_array(mfi); + + const Array4& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : + Array4{}; + const Array4& q1fx_z = Q1fx3->const_array(mfi); + + if (l_use_stretched_dz) { + ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoQ1_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + stretched_dz_d[level], q1fx_z, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } else if (l_use_terrain_fitted_coords) { + ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoQ1_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + z_nd_arr, detJ_arr, dxInv, q1fx_z, + mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } else { // no terrain + ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoQ1_comp, + slow_dt, l_bc_neumann_vals_d, cell_data, + dxInv, q1fx_z, mu_turb, solverChoice, + bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac); + } + } // mfi + } // if moist model and do implicit diff + } // implicit vertical factor > 0 diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index f9a3bc4b4..6aa167e8a 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -422,8 +422,11 @@ void erf_slow_rhs_post (int level, int finest_level, if (l_use_diff) { // Allow for implicit moisture diffusion - const Real l_vert_implicit_fac = (solverChoice.implicit_moisture_diffusion) ? - solverChoice.vert_implicit_fac[nrk] : zero; + Real l_vert_implicit_fac = zero; + if ( (ivar == RhoKE_comp && solverChoice.implicit_ke_diffusion ) || + (ivar == RhoQ1_comp && solverChoice.implicit_moisture_diffusion) ) { + l_vert_implicit_fac = solverChoice.vert_implicit_fac[nrk]; + } const Array4 tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4{};