Skip to content

Commit

Permalink
This seems to be working but dies with MOST due to 0 eddy viscosity i…
Browse files Browse the repository at this point in the history 
…n ghost cell.
  • Loading branch information
@AMLattanzi
AMLattanzi committed Nov 1, 2023
1 parent 7e9821f commit e1bbe47
Show file tree
Hide file tree
Showing 8 changed files with 188 additions and 96 deletions.
130 changes: 82 additions & 48 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,14 @@ public:
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd,
int ng_in)
: m_geom(geom), m_ma(geom,vars_old,Theta_prim,z_phys_nd)
amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>& sst_lev,
amrex::Vector<amrex::Vector<std::unique_ptr<amrex::iMultiFab>>>& lmask_lev,
amrex::Real& start_bdy_time,
amrex::Real& bdy_time_interval)
: m_start_bdy_time(start_bdy_time),
m_bdy_time_interval(bdy_time_interval),
m_geom(geom),
m_ma(geom,vars_old,Theta_prim,z_phys_nd)
{
amrex::ParmParse pp("erf");
pp.query("most.z0", z0_const);
Expand Down Expand Up @@ -78,95 +84,118 @@ public:
amrex::Abort("Undefined MOST roughness type!");
}

// Size the MOST params for all levels
int nlevs = m_geom.size();
z_0.resize(nlevs);
u_star.resize(nlevs);
t_star.resize(nlevs);
t_surf.resize(nlevs);
olen.resize(nlevs);

for (int lev = 0; lev < nlevs; lev++)
{
// Z0 heights
// Get pointers to SST and LANDMASK data
m_sst_lev.resize(nlevs);
m_lmask_lev.resize(nlevs);
for (int lev(0); lev<nlevs; ++lev) {
int nt_tot = sst_lev[lev].size();
m_sst_lev[lev].resize(nt_tot);
m_lmask_lev[lev].resize(nt_tot);
for (int nt(0); nt<nt_tot; ++nt) {
m_sst_lev[lev][nt] = sst_lev[lev][nt].get();
m_lmask_lev[lev][nt] = lmask_lev[lev][nt].get();
}
}

for (int lev = 0; lev < nlevs; lev++) {
// Attributes for MFs and FABs
//--------------------------------------------------------
amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng_in);
auto& mf = vars_old[lev][Vars::cons];
// Create a 2D ba, dm, & ghost cells
amrex::BoxArray ba = mf.boxArray();
amrex::BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
b.setRange(2,0);
}
amrex::BoxArray ba2d(std::move(bl2d));
const amrex::DistributionMapping& dm = mf.DistributionMap();
const int ncomp = 1;
amrex::IntVect ng = mf.nGrowVect(); ng[2]=0;

// Z0 heights FAB
//--------------------------------------------------------
amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng);
bx.setSmall(2,0);
bx.setBig(2,0);
z_0[lev].resize(bx,1);
z_0[lev].setVal<amrex::RunOn::Device>(z0_const);

// 2D MFs for U*, T*, T_surf
//--------------------------------------------------------
{ // CC vars
auto& mf = vars_old[lev][Vars::cons];
// Create a 2D ba, dm, & ghost cells
amrex::BoxArray ba = mf.boxArray();
amrex::BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
b.setRange(2,0);
}
amrex::BoxArray ba2d(std::move(bl2d));
const amrex::DistributionMapping& dm = mf.DistributionMap();
const int ncomp = 1;
amrex::IntVect ng = mf.nGrowVect(); ng[2]=0;

u_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
u_star[lev]->setVal(1.E34);

t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
t_star[lev]->setVal(1.E34);

olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
olen[lev]->setVal(1.E34);

t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
if (erf_st) {
t_surf[lev]->setVal(surf_temp);
} else {
t_surf[lev]->setVal(0.0);
}
u_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
u_star[lev]->setVal(1.E34);

t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
t_star[lev]->setVal(1.E34);

olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
olen[lev]->setVal(1.E34);

t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);

if (m_sst_lev[lev][0]) { // Valid SST data at t==0
theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
amrex::MultiFab::Copy(*(t_surf[lev]), *(m_sst_lev[lev][0]), 0, 0, 1, ng);
} else if (erf_st) { // Constant temp
t_surf[lev]->setVal(surf_temp);
} else {
t_surf[lev]->setVal(0.0);
}
}// lev
}

void
impose_most_bcs (const int lev,
update_fluxes (const int& lev,
const amrex::Real& time,
int max_iters = 25);

template <typename FluxIter>
void
compute_fluxes (const int& lev,
const int& max_iters,
const FluxIter& most_flux);

void
impose_most_bcs (const int& lev,
const amrex::Vector<amrex::MultiFab*>& mfs,
amrex::MultiFab* eddyDiffs);

template<typename FluxCalc>
void
compute_most_bcs (const int lev,
compute_most_bcs (const int& lev,
const amrex::Vector<amrex::MultiFab*>& mfs,
amrex::MultiFab* eddyDiffs,
const FluxCalc& flux_comp);

void
update_fluxes (int lev, int max_iters = 25);
time_interp_tsurf(const int& lev,
const amrex::Real& time);

template <typename FluxIter>
void
compute_fluxes (const int& lev,
const int& max_iters,
const FluxIter& most_flux);

void
update_mac_ptrs (int lev,
update_mac_ptrs (const int& lev,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim)
{ m_ma.update_field_ptrs(lev,vars_old,Theta_prim); }

const amrex::MultiFab*
get_u_star (int lev) { return u_star[lev].get(); }
get_u_star (const int& lev) { return u_star[lev].get(); }

const amrex::MultiFab*
get_t_star (int lev) { return t_star[lev].get(); }
get_t_star (const int& lev) { return t_star[lev].get(); }

const amrex::MultiFab*
get_olen (int lev) { return olen[lev].get(); }
get_olen (const int& lev) { return olen[lev].get(); }

const amrex::MultiFab*
get_mac_avg (int lev, int comp) { return m_ma.get_average(lev,comp); }
get_mac_avg (const int& lev, int comp) { return m_ma.get_average(lev,comp); }

enum struct FluxCalcType {
MOENG = 0, ///< Moeng functional form
Expand Down Expand Up @@ -195,6 +224,8 @@ private:
amrex::Real surf_temp_flux{0};
amrex::Real cnk_a{0.0185};
amrex::Real depth{30.0};
amrex::Real m_start_bdy_time;
amrex::Real m_bdy_time_interval;
amrex::Vector<amrex::Geometry> m_geom;
amrex::Vector<amrex::FArrayBox> z_0;

Expand All @@ -203,6 +234,9 @@ private:
amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_star;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> olen;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_surf;

amrex::Vector<amrex::Vector<amrex::MultiFab*>> m_sst_lev;
amrex::Vector<amrex::Vector<amrex::iMultiFab*>> m_lmask_lev;
};

#endif /* ABLMOST_H */
55 changes: 43 additions & 12 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,13 @@ using namespace amrex;
* @param[in] max_iters maximum iterations to use
*/
void
ABLMost::update_fluxes (int lev,
ABLMost::update_fluxes (const int& lev,
const Real& time,
int max_iters)
{
// Update SST data if we have a valid pointer
if (m_sst_lev[lev][0]) time_interp_tsurf(lev, time);

// Compute plane averages for all vars (regardless of flux type)
m_ma.compute_averages(lev);

Expand Down Expand Up @@ -73,12 +77,9 @@ ABLMost::compute_fluxes (const int& lev,
const auto *const tm_ptr = m_ma.get_average(lev,2);
const auto *const umm_ptr = m_ma.get_average(lev,3);

// Ghost cells
amrex::IntVect ng = u_star[lev]->nGrowVect(); ng[2]=0;

for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
{
amrex::Box gtbx = mfi.growntilebox(ng);
Box gtbx = mfi.growntilebox();

auto u_star_arr = u_star[lev]->array(mfi);
auto t_star_arr = t_star[lev]->array(mfi);
Expand Down Expand Up @@ -106,7 +107,7 @@ ABLMost::compute_fluxes (const int& lev,
* @param[in] eddyDiffs Diffusion coefficients from turbulence model
*/
void
ABLMost::impose_most_bcs (const int lev,
ABLMost::impose_most_bcs (const int& lev,
const Vector<MultiFab*>& mfs,
MultiFab* eddyDiffs)
{
Expand All @@ -131,7 +132,7 @@ ABLMost::impose_most_bcs (const int lev,
*/
template<typename FluxCalc>
void
ABLMost::compute_most_bcs (const int lev,
ABLMost::compute_most_bcs (const int& lev,
const Vector<MultiFab*>& mfs,
MultiFab* eddyDiffs,
const FluxCalc& flux_comp)
Expand Down Expand Up @@ -168,7 +169,7 @@ 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;

Expand All @@ -178,9 +179,9 @@ ABLMost::compute_most_bcs (const int lev,
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)
Expand All @@ -189,9 +190,9 @@ ABLMost::compute_most_bcs (const int lev,
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)
Expand All @@ -203,3 +204,33 @@ ABLMost::compute_most_bcs (const int lev,
} // var_idx
} // mfiter
}

void
ABLMost::time_interp_tsurf(const int& lev,
const Real& time)
{
// Time interpolation
Real dT = m_bdy_time_interval;
Real time_since_start = time - m_start_bdy_time;
int n_time = static_cast<int>( time_since_start / dT);
amrex::Real alpha = (time_since_start - n_time * dT) / dT;
AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
amrex::Real oma = 1.0 - alpha;
AMREX_ALWAYS_ASSERT( (n_time >= 0) && (n_time < (m_sst_lev[lev].size()-1)));

// Populate t_surf
for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi)
{
Box gtbx = mfi.growntilebox();

auto t_surf_arr = t_surf[lev]->array(mfi);
const auto sst_hi_arr = m_sst_lev[lev][n_time+1]->const_array(mfi);
const auto sst_lo_arr = m_sst_lev[lev][n_time ]->const_array(mfi);

ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k)
+ alpha * sst_hi_arr(i,j,k);
});
}
}
11 changes: 5 additions & 6 deletions Source/BoundaryConditions/BoundaryConditions_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,8 @@ ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs,

// Time interpolation
Real dT = bdy_time_interval;
// amrex::Print() << " BoundaryConditions_metgrid.cpp ERF::fill_from_metgrid \ttime \t" << time << std::endl;
Real time_since_start = time;
int n_time = static_cast<int>( time_since_start / dT);
Real time_since_start = time - start_bdy_time;
int n_time = static_cast<int>( time_since_start / dT);
amrex::Real alpha = (time_since_start - n_time * dT) / dT;
AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0);
amrex::Real oma = 1.0 - alpha;
Expand All @@ -34,15 +33,15 @@ ERF::fill_from_metgrid (const Vector<MultiFab*>& mfs,
#if defined(ERF_USE_MOISTURE)
// Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQt RhoQp NumVars]
Vector<int> cons_read = {1, 1, 0, 0, 0, 1, 0};
Vector<int> cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0};
Vector<int> cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0};
#elif defined(ERF_USE_WARM_NO_PRECIP)
// Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQv RhoQc NumVars]
Vector<int> cons_read = {1, 1, 0, 0, 0, 1, 0};
Vector<int> cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0};
Vector<int> cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0};
# else
// Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar NumVars]
Vector<int> cons_read = {1, 1, 0, 0, 0};
Vector<int> cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0};
Vector<int> cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0};
#endif

Vector<Vector<int>> is_read;
Expand Down
14 changes: 9 additions & 5 deletions Source/BoundaryConditions/MOSTAverage.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,6 +262,9 @@ MOSTAverage::set_k_indices_N()
AMREX_ALWAYS_ASSERT(lk >= m_radius);

m_k_indx[lev]->setVal(lk);

amrex::Print() << "MA: " << m_zref << ' ' << m_zlo << ' ' << m_zhi << ' ' << m_dz << ' ' << lk << "\n";
exit(0);
}
// Specified k_indx & compute z_ref
} else if (read_k) {
Expand Down Expand Up @@ -304,12 +307,13 @@ MOSTAverage::set_k_indices_T()
auto k_arr = m_k_indx[lev]->array(mfi);
ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
k_arr(i,j,k) = 0;
Real z_target = d_zref + z_phys_arr(i,j,k);
for (int lk(0); lk<=kmax; ++lk) {
Real z_lo = 0.25 * ( z_phys_arr(i,j ,lk ) + z_phys_arr(i+1,j ,lk )
+ z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) );
+ z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) );
Real z_hi = 0.25 * ( z_phys_arr(i,j ,lk+1) + z_phys_arr(i+1,j ,lk+1)
+ z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) );
+ z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) );
if (z_target > z_lo && z_target < z_hi){
AMREX_ASSERT_WITH_MESSAGE(lk >= d_radius,
"K index must be larger than averaging radius!");
Expand Down Expand Up @@ -376,9 +380,9 @@ MOSTAverage::set_norm_indices_T()
Real z_target = delta_z + z_phys_arr(i,j,k);
for (int lk(0); lk<=kmax; ++lk) {
Real z_lo = 0.25 * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk )
+ z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
+ z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) );
Real z_hi = 0.25 * ( z_phys_arr(i_new,j_new ,lk+1) + z_phys_arr(i_new+1,j_new ,lk+1)
+ z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
+ z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) );
if (z_target > z_lo && z_target < z_hi){
AMREX_ASSERT_WITH_MESSAGE(lk >= d_radius,
"K index must be larger than averaging radius!");
Expand Down Expand Up @@ -532,7 +536,7 @@ MOSTAverage::compute_plane_averages(int lev)
// Peel back the level
auto& fields = m_fields[lev];
auto& averages = m_averages[lev];
const auto & geom = m_geom[lev];
const auto & geom = m_geom[lev];

auto& z_phys = m_z_phys_nd[lev];
auto& x_pos = m_x_pos[lev];
Expand Down
Loading

0 comments on commit e1bbe47

Please sign in to comment.