From bfcaf08c6ec9d848581d1605c16554e47e52e329 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Wed, 25 Oct 2023 10:54:11 -0700 Subject: [PATCH] Fix MOST with grid stretching in z. --- Source/BoundaryConditions/ABLMost.H | 5 +- Source/BoundaryConditions/ABLMost.cpp | 53 +++++++++++++-------- Source/BoundaryConditions/ERF_FillPatch.cpp | 3 +- Source/BoundaryConditions/MOSTStress.H | 18 +++---- Source/ERF.H | 2 +- Source/TimeIntegration/TI_utils.H | 4 +- 6 files changed, 51 insertions(+), 34 deletions(-) diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 77d8780a2..73d3f26ea 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -132,13 +132,16 @@ public: void impose_most_bcs (const int lev, const amrex::Vector& mfs, - amrex::MultiFab* eddyDiffs); + amrex::MultiFab* eddyDiffs, + amrex::MultiFab* z_phys); template void compute_most_bcs (const int lev, const amrex::Vector& mfs, amrex::MultiFab* eddyDiffs, + amrex::MultiFab* z_phys, + const amrex::Real& dz_no_terrain, const FluxCalc& flux_comp); void diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index bc270cb8c..96bbadc01 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -108,15 +108,16 @@ ABLMost::compute_fluxes (const int& lev, void ABLMost::impose_most_bcs (const int lev, const Vector& mfs, - MultiFab* eddyDiffs) + MultiFab* eddyDiffs, + MultiFab* z_phys) { const int zlo = 0; if (flux_type == FluxCalcType::DONELAN) { - donelan_flux flux_comp(zlo,m_geom[lev].CellSize(2)); - compute_most_bcs(lev,mfs,eddyDiffs,flux_comp); + donelan_flux flux_comp(zlo); + compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp); } else { - moeng_flux flux_comp(zlo,m_geom[lev].CellSize(2)); - compute_most_bcs(lev,mfs,eddyDiffs,flux_comp); + moeng_flux flux_comp(zlo); + compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp); } } @@ -134,6 +135,8 @@ void ABLMost::compute_most_bcs (const int lev, const Vector& mfs, MultiFab* eddyDiffs, + MultiFab* z_phys, + const Real& dz_no_terrain, const FluxCalc& flux_comp) { const int zlo = 0; @@ -141,10 +144,11 @@ ABLMost::compute_most_bcs (const int lev, for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi) { // Get field arrays - const auto cons_arr = mfs[Vars::cons]->array(mfi); - const auto velx_arr = mfs[Vars::xvel]->array(mfi); - const auto vely_arr = mfs[Vars::yvel]->array(mfi); - const auto eta_arr = eddyDiffs->array(mfi); + const auto cons_arr = mfs[Vars::cons]->array(mfi); + const auto velx_arr = mfs[Vars::xvel]->array(mfi); + const auto vely_arr = mfs[Vars::yvel]->array(mfi); + const auto eta_arr = eddyDiffs->array(mfi); + const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4{}; // Get average arrays const auto *const u_mean = m_ma.get_average(lev,0); @@ -168,36 +172,45 @@ ABLMost::compute_most_bcs (const int lev, auto dest_arr = (*mfs[var_idx])[mfi].array(); if (var_idx == Vars::cons) { - amrex::Box b2d = bx; // Copy constructor + Box b2d = bx; b2d.setBig(2,zlo-1); int n = Cons::RhoTheta; ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - flux_comp.compute_t_flux(i,j,k,n,icomp,cons_arr,velx_arr,vely_arr,eta_arr, - umm_arr,tm_arr,u_star_arr,t_star_arr,t_surf_arr,dest_arr); + Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain; + flux_comp.compute_t_flux(i, j, k, n, icomp, dz, + cons_arr, velx_arr, vely_arr, eta_arr, + umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr, + dest_arr); }); - } else if (var_idx == Vars::xvel || var_idx == Vars::xmom) { //for velx + } else if (var_idx == Vars::xvel || var_idx == Vars::xmom) { - amrex::Box xb2d = surroundingNodes(bx,0); // Copy constructor + Box xb2d = surroundingNodes(bx,0); xb2d.setBig(2,zlo-1); ParallelFor(xb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - flux_comp.compute_u_flux(i,j,k,icomp,var_idx,cons_arr,velx_arr,vely_arr,eta_arr, - umm_arr,um_arr,u_star_arr,dest_arr); + Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain; + flux_comp.compute_u_flux(i, j, k, icomp, var_idx, dz, + cons_arr, velx_arr, vely_arr, eta_arr, + umm_arr, um_arr, u_star_arr, + dest_arr); }); - } else if (var_idx == Vars::yvel || var_idx == Vars::ymom) { //for vely + } else if (var_idx == Vars::yvel || var_idx == Vars::ymom) { - amrex::Box yb2d = surroundingNodes(bx,1); // Copy constructor + Box yb2d = surroundingNodes(bx,1); yb2d.setBig(2,zlo-1); ParallelFor(yb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - flux_comp.compute_v_flux(i,j,k,icomp,var_idx,cons_arr,velx_arr,vely_arr,eta_arr, - umm_arr,vm_arr,u_star_arr,dest_arr); + Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain; + flux_comp.compute_v_flux(i, j, k, icomp, var_idx, dz, + cons_arr, velx_arr, vely_arr, eta_arr, + umm_arr, vm_arr, u_star_arr, + dest_arr); }); } } // var_idx diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index ba9cc3c7e..7d2a4e4ed 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -158,7 +158,6 @@ ERF::FillIntermediatePatch (int lev, Real time, const Vector& mfs, int ng_cons, int ng_vel, bool cons_only, int icomp_cons, int ncomp_cons, - MultiFab* eddyDiffs, bool allow_most_bcs) { BL_PROFILE_VAR("FillIntermediatePatch()",FillIntermediatePatch); @@ -260,7 +259,7 @@ ERF::FillIntermediatePatch (int lev, Real time, // MOST boundary conditions if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs) - m_most->impose_most_bcs(lev,mfs,eddyDiffs); + m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get()); } // diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 7c021d7d9..d6661705c 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -587,9 +587,8 @@ private: */ struct moeng_flux { - moeng_flux (int l_zlo, - amrex::Real l_dz) - : zlo(l_zlo), dz(l_dz) {} + moeng_flux (int l_zlo) + : zlo(l_zlo) {} AMREX_GPU_DEVICE @@ -600,6 +599,7 @@ struct moeng_flux const int& k, const int& n, const int& icomp, + const amrex::Real& dz, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, @@ -663,6 +663,7 @@ struct moeng_flux const int& k, const int& icomp, const int& var_idx, + const amrex::Real& dz, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, @@ -725,6 +726,7 @@ struct moeng_flux const int& k, const int& icomp, const int& var_idx, + const amrex::Real& dz, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, @@ -781,7 +783,6 @@ struct moeng_flux private: int zlo; - amrex::Real dz; }; @@ -790,9 +791,8 @@ private: */ struct donelan_flux { - donelan_flux (int l_zlo, - amrex::Real l_dz) - : zlo(l_zlo), dz(l_dz) {} + donelan_flux (int l_zlo) + : zlo(l_zlo) {} AMREX_GPU_DEVICE @@ -803,6 +803,7 @@ struct donelan_flux const int& k, const int& n, const int& icomp, + const amrex::Real& dz, const amrex::Array4& cons_arr, const amrex::Array4& /*velx_arr*/, const amrex::Array4& /*vely_arr*/, @@ -849,6 +850,7 @@ struct donelan_flux const int& k, const int& icomp, const int& var_idx, + const amrex::Real& dz, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, @@ -915,6 +917,7 @@ struct donelan_flux const int& k, const int& icomp, const int& var_idx, + const amrex::Real& dz, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, @@ -975,6 +978,5 @@ struct donelan_flux private: int zlo; - amrex::Real dz; }; #endif diff --git a/Source/ERF.H b/Source/ERF.H index 02beaddb7..ca240c279 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -369,7 +369,7 @@ private: void FillIntermediatePatch (int lev, amrex::Real time, const amrex::Vector& mfs, int ng_cons, int ng_vel, bool cons_only, int icomp_cons, int ncomp_cons, - amrex::MultiFab* eddyDiffs, bool allow_most_bcs = true); + bool allow_most_bcs = true); // Fill all multifabs (and all components) in a vector of multifabs corresponding to the // grid variables defined in vars_old and vars_new just as FillCoarsePatch. diff --git a/Source/TimeIntegration/TI_utils.H b/Source/TimeIntegration/TI_utils.H index 664357b82..88b0b587a 100644 --- a/Source/TimeIntegration/TI_utils.H +++ b/Source/TimeIntegration/TI_utils.H @@ -80,7 +80,7 @@ FillIntermediatePatch(level, time_for_fp, {&S_data[IntVar::cons], &xvel_new, &yvel_new, &zvel_new}, - ng_cons_to_use, 0, cons_only, scomp_cons, ncomp_cons, eddyDiffs); + ng_cons_to_use, 0, cons_only, scomp_cons, ncomp_cons); // Here we don't use include any of the ghost region because we have only updated // momentum on valid faces @@ -132,7 +132,7 @@ FillIntermediatePatch(level, time_for_fp, {&S_data[IntVar::cons], &xvel_new, &yvel_new, &zvel_new}, ng_cons_to_use, ng_vel, cons_only, scomp_cons, ncomp_cons, - eddyDiffs, allow_most_bcs); + allow_most_bcs); // Now we can convert back to momentum on valid+ghost since we have // filled the ghost regions for both velocity and density