From 17ba397c2c8e1832436d6947de799ffc1ffd92a8 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 2 Nov 2023 15:38:26 -0700 Subject: [PATCH] Turbulent viscosity corner cells (#1281) * Corner cells outside of the domain were not copied correctly. * Fix tolerance on GPU. --- Source/BoundaryConditions/MOSTStress.H | 27 ++++++++++++------- .../Diffusion/ComputeTurbulentViscosity.cpp | 20 +++++++------- 2 files changed, 28 insertions(+), 19 deletions(-) diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 1e4204192..c654546f0 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -149,7 +149,7 @@ struct adiabatic_charnock private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -202,7 +202,7 @@ struct adiabatic_mod_charnock private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -260,7 +260,7 @@ struct surface_flux private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -323,7 +323,7 @@ struct surface_flux_charnock private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -387,7 +387,7 @@ struct surface_flux_mod_charnock private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -447,7 +447,7 @@ struct surface_temp private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -512,7 +512,7 @@ struct surface_temp_charnock private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -578,7 +578,7 @@ struct surface_temp_mod_charnock private: most_data mdata; similarity_funs sfuns; - amrex::Real tol = 1.0e-5; + const amrex::Real tol = 1.0e-5; }; @@ -639,6 +639,7 @@ struct moeng_flux rho = cons_arr(ic,jc,zlo,Rho_comp); theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] + eta = amrex::max(eta,eps); amrex::Real theta_mean = tm_arr(ic,jc,zlo); amrex::Real wsp_mean = umm_arr(ic,jc,zlo); @@ -649,7 +650,7 @@ struct moeng_flux amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real num1 = (theta-theta_mean)*wsp_mean; amrex::Real num2 = (theta_mean-theta_surf)*wsp; - amrex::Real moflux = (std::abs(tstar) > EPS) ? + amrex::Real moflux = (std::abs(tstar) > eps) ? tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0; amrex::Real deltaz = dz * (zlo - k); @@ -700,6 +701,7 @@ struct moeng_flux + cons_arr(ic ,jc,zlo,Rho_comp) ); eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eps); amrex::Real umean = um_arr(i,j,zlo); amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) ); @@ -763,6 +765,7 @@ struct moeng_flux + cons_arr(ic,jc ,zlo,Rho_comp) ); eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eps); amrex::Real vmean = vm_arr(i,j,zlo); amrex::Real wsp_mean = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) ); @@ -784,7 +787,7 @@ struct moeng_flux private: int zlo; - static constexpr amrex::Real EPS = 1e-16; + const amrex::Real eps = 1e-16; }; @@ -833,6 +836,7 @@ struct donelan_flux rho = cons_arr(ic,jc,zlo,Rho_comp); theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] + eta = amrex::max(eta,eps); amrex::Real Cd = 0.0012; amrex::Real wsp_mean = umm_arr(ic,jc,zlo); @@ -888,6 +892,7 @@ struct donelan_flux + cons_arr(ic ,jc,zlo,Rho_comp) ); eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eps); amrex::Real Cd = 0.001; const amrex::Real c = 7e-5; @@ -955,6 +960,7 @@ struct donelan_flux + cons_arr(ic,jc ,zlo,Rho_comp) ); eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eps); amrex::Real Cd = 0.001; const amrex::Real c = 7e-5; @@ -980,5 +986,6 @@ struct donelan_flux private: int zlo; + const amrex::Real eps = 1e-16; }; #endif diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 4827610ff..bf3709cf4 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -193,7 +193,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu for ( amrex::MFIter mfi(eddyViscosity,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box bxcc = mfi.tilebox(); - Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); + Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1); Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc); int i_lo = bxcc.smallEnd(0); int i_hi = bxcc.bigEnd(0); int j_lo = bxcc.smallEnd(1); int j_hi = bxcc.bigEnd(1); @@ -202,19 +202,21 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu const Array4& mu_turb = eddyViscosity.array(mfi); - // Extrapolate outside the domain in lateral directions + // Extrapolate outside the domain in lateral directions (planex owns corner cells) if (i_lo == domain.smallEnd(0)) { amrex::ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - mu_turb(i_lo-i, j, k, EddyDiff::Mom_h) = mu_turb(i_lo, j, k, EddyDiff::Mom_h); - mu_turb(i_lo-i, j, k, EddyDiff::Mom_v) = mu_turb(i_lo, j, k, EddyDiff::Mom_v); + int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1)); + mu_turb(i_lo-i, j, k, EddyDiff::Mom_h) = mu_turb(i_lo, lj, k, EddyDiff::Mom_h); + mu_turb(i_lo-i, j, k, EddyDiff::Mom_v) = mu_turb(i_lo, lj, k, EddyDiff::Mom_v); }); } if (i_hi == domain.bigEnd(0)) { amrex::ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - mu_turb(i_hi+i, j, k, EddyDiff::Mom_h) = mu_turb(i_hi, j, k, EddyDiff::Mom_h); - mu_turb(i_hi+i, j, k, EddyDiff::Mom_v) = mu_turb(i_hi, j, k, EddyDiff::Mom_v); + int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1)); + mu_turb(i_hi+i, j, k, EddyDiff::Mom_h) = mu_turb(i_hi, lj, k, EddyDiff::Mom_h); + mu_turb(i_hi+i, j, k, EddyDiff::Mom_v) = mu_turb(i_hi, lj, k, EddyDiff::Mom_v); }); } if (j_lo == domain.smallEnd(1)) { @@ -336,9 +338,9 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v); }); break; - } - } - } + } + } + } } /**