Skip to content

Commit

Permalink
Turbulent viscosity corner cells (#1281)
Browse files Browse the repository at this point in the history 
* Corner cells outside of the domain were not copied correctly.

* Fix tolerance on GPU.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Nov 2, 2023
1 parent 885a8c4 commit 17ba397
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 19 deletions.
27 changes: 17 additions & 10 deletions Source/BoundaryConditions/MOSTStress.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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;
};


Expand Down Expand Up @@ -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);
Expand All @@ -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);

Expand Down Expand Up @@ -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) );
Expand Down Expand Up @@ -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) );
Expand All @@ -784,7 +787,7 @@ struct moeng_flux

private:
int zlo;
static constexpr amrex::Real EPS = 1e-16;
const amrex::Real eps = 1e-16;
};


Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -980,5 +986,6 @@ struct donelan_flux

private:
int zlo;
const amrex::Real eps = 1e-16;
};
#endif
20 changes: 11 additions & 9 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -202,19 +202,21 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu

const Array4<Real>& 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)) {
Expand Down Expand Up @@ -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;
}
}
}
}
}
}
}

/**
Expand Down

0 comments on commit 17ba397

Please sign in to comment.