ERF_ImplicitPre.H

Go to the documentation of this file.

        // *****************************************************************************
        // Do semi-implicit solve for diffusion: theta, u, and v
        // *****************************************************************************
        {
            const Real l_vert_implicit_fac = solverChoice.vert_implicit_fac[level][nrk];
 
            if (l_vert_implicit_fac > zero) {
 
                const bool l_do_implicit_theta = solverChoice.implicit_thermal_diffusion;
                const bool l_do_implicit_mom   = solverChoice.implicit_momentum_diffusion;
 
                // If we're doing an implicit solve for momenta (u and v only),
                // then we use the explicit solution at this point to derive
                // the surface velocity gradient -- to be consistent with the
                // hoextrap used with the surface layer BC. We then subtract
                // out the contribution from the gradient of tau_corr times
                // (1 - implicit_fac) and add in the implicit solution scaled
                // by implicit_fac.
                MultiFab* Tau13corr = (l_do_implicit_mom) ? Tau_corr[level][0].get() : nullptr;
                MultiFab* Tau23corr = (l_do_implicit_mom) ? Tau_corr[level][1].get() : nullptr;
 
                // BCs
                const BCRec* bc_ptr_h = domain_bcs_type.data();
                GpuArray<Real, AMREX_SPACEDIM*2> l_bc_neumann_vals_d;
                for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) {
                    l_bc_neumann_vals_d[ori] = m_bc_neumann_vals[RhoTheta_comp][ori];
                }
                const bool l_use_SurfLayer = (m_SurfaceLayer[Orientation(Direction::z, Orientation::low)] != nullptr);
 
                const bool l_use_turb       = solverChoice.turbChoice[level].use_kturb;
 
                const bool l_use_stretched_dz = (solverChoice.mesh_type == MeshType::StretchedDz);
 
                for ( MFIter mfi(S_old[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
                {
                    Box bx  = mfi.tilebox();
                    Box tbx = mfi.nodaltilebox(0);
                    Box tby = mfi.nodaltilebox(1);
 
                    const Array4<      Real>& cell_data  = scratch.array(mfi);
 
                    const Array4<      Real>& rho_u      = (l_do_implicit_mom) ? scratch_xmom.array(mfi) : Array4<Real>{};
                    const Array4<      Real>& rho_v      = (l_do_implicit_mom) ? scratch_ymom.array(mfi) : Array4<Real>{};
                    const Array4<const Real>& tau13_corr = (l_do_implicit_mom) ? Tau13corr->array(mfi) : Array4<Real>{};
                    const Array4<const Real>& tau23_corr = (l_do_implicit_mom) ? Tau23corr->array(mfi) : Array4<Real>{};
 
 
                    const Array4<const Real>& z_nd_arr = z_phys_nd[level]->const_array(mfi);
                    const Array4<const Real>& detJ_arr = detJ_cc[level]->const_array(mfi);
 
                    // erf.implicit_terrain_metric (terrain-fitted meshes): the map factors the
                    // terrain-metric coefficient needs at the cells and at the u and v faces
                    const bool l_implicit_metric = solverChoice.implicit_terrain_metric &&
                                                   l_use_terrain_fitted_coords;
                    const Array4<const Real> mf_mx = mapfac[level][MapFacType::m_x]->const_array(mfi);
                    const Array4<const Real> mf_my = mapfac[level][MapFacType::m_y]->const_array(mfi);
                    const Array4<const Real> mf_ux = mapfac[level][MapFacType::u_x]->const_array(mfi);
                    const Array4<const Real> mf_uy = mapfac[level][MapFacType::u_y]->const_array(mfi);
                    const Array4<const Real> mf_vx = mapfac[level][MapFacType::v_x]->const_array(mfi);
                    const Array4<const Real> mf_vy = mapfac[level][MapFacType::v_y]->const_array(mfi);
 
                    const Array4<const Real>& mu_turb = l_use_turb ? eddyDiffs->const_array(mfi) : Array4<const Real>{};
 
                    // Vertical extent of the grid column over each (i,j), with a one-cell halo
                    // so a box can also see the column on the other side of each of its faces
                    const Array4<const int>& col_kext = column_kextent[level]->const_array(mfi);
 
                    const Array4<const Real> tau13  = Tau[level][TauType::tau13]->array(mfi);
                    const Array4<const Real> tau23  = Tau[level][TauType::tau23]->array(mfi);
                    [[maybe_unused]] const Array4<const Real> tau33  = Tau[level][TauType::tau33]->array(mfi);
                    const Array4<const Real>& hfx_z = Hfx3->const_array(mfi);
                    const bool l_use_mrf_cg = solverChoice.turbChoice[level].enable_mrf_countergradient
                                            && (solverChoice.turbChoice[level].pbl_type == PBLType::MRF);
                    // NOTE: Only ComputeDiffusivityYSUNew fills EddyDiff::HGAMU_v/HGAMV_v;
                    //       the legacy PBLType::YSU never writes them.  Gating on YSU
                    //       therefore both disabled the term where it is produced and read
                    //       an unwritten component where it was not.
                    const bool l_use_ysu_mom_cg = solverChoice.turbChoice[level].enable_ysu_countergradient
                                               && (solverChoice.turbChoice[level].pbl_type == PBLType::YSUNew);
                    if (l_use_stretched_dz) {
                        if (l_do_implicit_theta) {
 
                            ImplicitDiffForStateLU_S(bx, fine_geom.Domain(), level, RhoTheta_comp,
                                                     stage_dt, l_bc_neumann_vals_d, cell_data,
                                                     stretched_dz_d[level], hfx_z,
                                                     mu_turb, solverChoice,
                                                     bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                     l_use_mrf_cg);
                        }
                        if (l_do_implicit_mom) {
                            ImplicitDiffForMomLU_S<0>(tbx, fine_geom.Domain(), level, stage_dt,
                                                      col_kext,
                                                      cell_data, rho_u, tau13, tau13_corr,
                                                      stretched_dz_d[level],
                                                      mu_turb, solverChoice,
                                                      bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                      l_use_ysu_mom_cg);
 
                            ImplicitDiffForMomLU_S<1>(tby, fine_geom.Domain(), level, stage_dt,
                                                      col_kext,
                                                      cell_data, rho_v, tau23, tau23_corr,
                                                      stretched_dz_d[level],
                                                      mu_turb, solverChoice,
                                                      bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                      l_use_ysu_mom_cg);
                        }
                    } else if (l_use_terrain_fitted_coords) {
                        if (l_do_implicit_theta) {
                            ImplicitDiffForStateLU_T(bx, fine_geom.Domain(), level, RhoTheta_comp,
                                                     stage_dt, l_bc_neumann_vals_d, cell_data,
                                                     z_nd_arr, detJ_arr, dxInv, hfx_z,
                                                     mu_turb, solverChoice,
                                                     bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                     l_use_mrf_cg,
                                                     mf_mx, mf_my, l_implicit_metric);
                        }
                        if (l_do_implicit_mom) {
                            ImplicitDiffForMomLU_T<0>(tbx, fine_geom.Domain(), level, stage_dt,
                                                      col_kext,
                                                      cell_data, rho_u, tau13, tau13_corr,
                                                      z_nd_arr, detJ_arr, dxInv,
                                                      mu_turb, solverChoice,
                                                      bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                      l_use_ysu_mom_cg,
                                                      mf_ux, mf_uy, l_implicit_metric);
 
                            ImplicitDiffForMomLU_T<1>(tby, fine_geom.Domain(), level, stage_dt,
                                                      col_kext,
                                                      cell_data, rho_v, tau23, tau23_corr,
                                                      z_nd_arr, detJ_arr, dxInv,
                                                      mu_turb, solverChoice,
                                                      bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                      l_use_ysu_mom_cg,
                                                      mf_vx, mf_vy, l_implicit_metric);
                        }
                    } else { // no terrain
                        if (l_do_implicit_theta) {
                            ImplicitDiffForStateLU_N(bx, fine_geom.Domain(), level, RhoTheta_comp,
                                                     stage_dt, l_bc_neumann_vals_d, cell_data,
                                                     dxInv, hfx_z, mu_turb, solverChoice,
                                                     bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                     l_use_mrf_cg,
                                                     (ib_wall_face33[level]) ? ib_wall_face33[level]->const_array(mfi)
                                                                             : Array4<const Real>{});
                        }
                        if (l_do_implicit_mom) {
                            // Fraction-stress immersed wall law: no coupling through the wall face
                            const Array4<const Real> wall13 = (ib_wall_face13[level])
                                ? ib_wall_face13[level]->const_array(mfi) : Array4<const Real>{};
                            const Array4<const Real> wall23 = (ib_wall_face23[level])
                                ? ib_wall_face23[level]->const_array(mfi) : Array4<const Real>{};
                            ImplicitDiffForMomLU_N<0>(tbx, fine_geom.Domain(), level, stage_dt,
                                                      col_kext,
                                                      cell_data, rho_u, tau13, tau13_corr,
                                                      dxInv,
                                                      mu_turb, solverChoice,
                                                      bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                      l_use_ysu_mom_cg, wall13);
 
                            ImplicitDiffForMomLU_N<1>(tby, fine_geom.Domain(), level, stage_dt,
                                                      col_kext,
                                                      cell_data, rho_v, tau23, tau23_corr,
                                                      dxInv,
                                                      mu_turb, solverChoice,
                                                      bc_ptr_h, l_use_SurfLayer, l_vert_implicit_fac,
                                                      l_use_ysu_mom_cg, wall23);
                        }
                    }
                } // mfi
 
            } // if do implicit solve for diffusive contribution to (rho theta) update
        } // wrapper