Skip to content

Commit

Permalink
Modifications for BOMEX (erf-model#1460)
Browse files Browse the repository at this point in the history 
* This compiled in debug mode.

* Verified with debug run.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Feb 28, 2024
1 parent 5da3c98 commit 2e78b23
Show file tree
Hide file tree
Showing 10 changed files with 511 additions and 39 deletions.
70 changes: 70 additions & 0 deletions Exec/ABL/inputs_most_bomex
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 4000

amrex.fpe_trap_invalid = 1

fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 1024 1024 1024
amr.n_cell = 64 64 64

geometry.is_periodic = 1 1 0

# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA)
zlo.type = "Most"
erf.most.flux_type = "custom"
erf.most.surf_mom_flux = 0.28 # ustar
erf.most.surf_temp_flux = 8.0e-3 # theta flux
erf.most.surf_moist_flux = 5.2e-5 # qv flux
erf.most.z0 = 0.1
erf.most.zref = 8.0

zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.fixed_dt = 0.1 # fixed time step depending on grid resolution

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
erf.v = 1 # verbosity in ERF.cpp
amr.v = 1 # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = 100 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 10 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta qv qc qp

# SOLVER CHOICE
erf.alpha_T = 0.0
erf.alpha_C = 1.0
erf.use_gravity = false

erf.molec_diff_type = "None"
erf.les_type = "Smagorinsky"
erf.Cs = 0.1

erf.init_type = "uniform"

erf.moisture_model = "Kessler"

# PROBLEM PARAMETERS
prob.rho_0 = 1.0
prob.A_0 = 1.0
prob.T_0 = 300.0
prob.U_0 = 10.0
prob.V_0 = 0.0
prob.W_0 = 0.0

# Higher values of perturbations lead to instability
# Instability seems to be coming from BC
prob.U_0_Pert_Mag = 0.08
prob.V_0_Pert_Mag = 0.08 #
prob.W_0_Pert_Mag = 0.0
67 changes: 48 additions & 19 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ public:
explicit ABLMost (const amrex::Vector<amrex::Geometry>& geom,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& 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,
Expand All @@ -42,8 +43,12 @@ public:
: 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)
m_ma(geom,vars_old,Theta_prim,Qv_prim,z_phys_nd)
{
// We have a moisture model if Qv_prim is a valid pointer
use_moisture = (Qv_prim[0].get());

// Get roughness
amrex::ParmParse pp("erf");
pp.query("most.z0", z0_const);

Expand All @@ -54,28 +59,42 @@ public:
flux_type = FluxCalcType::DONELAN;
} else if (flux_string == "moeng") {
flux_type = FluxCalcType::MOENG;
} else if (flux_string == "custom") {
flux_type = FluxCalcType::CUSTOM;
} else {
amrex::Abort("Undefined MOST flux type!");
}

// Specify surface temperature or surface flux
// Get surface temperature
auto erf_st = pp.query("most.surf_temp", surf_temp);
if (erf_st) {
theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
surf_heating_rate = surf_heating_rate / 3600.0; // [K/s]
if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
}

// Custom type user must specify the fluxes
if (flux_type == FluxCalcType::CUSTOM) {
AMREX_ASSERT_WITH_MESSAGE(use_moisture, "Moisture model must be used with MOST CUSTOM!");
theta_type = ThetaCalcType::HEAT_FLUX;
pp.query("most.surf_mom_flux" , surf_mom_flux );
pp.query("most.surf_temp_flux" , surf_temp_flux );
pp.query("most.surf_moist_flux", surf_moist_flux);

// Specify surface temperature or surface flux
} else {
pp.query("most.surf_temp_flux", surf_temp_flux);
if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
}
if (std::abs(surf_temp_flux) > std::numeric_limits<amrex::Real>::epsilon()) {
theta_type = ThetaCalcType::HEAT_FLUX;
if (erf_st) {
theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h]
surf_heating_rate = surf_heating_rate / 3600.0; // [K/s]
if (pp.query("most.surf_temp_flux", surf_temp_flux)) {
amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
}
} else {
theta_type = ThetaCalcType::ADIABATIC;
pp.query("most.surf_temp_flux", surf_temp_flux);
if (pp.query("most.surf_heating_rate", surf_heating_rate)) {
amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate");
}
if (std::abs(surf_temp_flux) > std::numeric_limits<amrex::Real>::epsilon()) {
theta_type = ThetaCalcType::HEAT_FLUX;
} else {
theta_type = ThetaCalcType::ADIABATIC;
}
}
}

Expand All @@ -99,6 +118,7 @@ public:
z_0.resize(nlevs);
u_star.resize(nlevs);
t_star.resize(nlevs);
q_star.resize(nlevs);
t_surf.resize(nlevs);
olen.resize(nlevs);

Expand Down Expand Up @@ -159,6 +179,9 @@ public:
t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
t_star[lev]->setVal(1.E34);

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

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

Expand Down Expand Up @@ -229,8 +252,9 @@ public:
void
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); }
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim)
{ m_ma.update_field_ptrs(lev,vars_old,Theta_prim,Qv_prim); }

const amrex::MultiFab*
get_u_star (const int& lev) { return u_star[lev].get(); }
Expand All @@ -246,7 +270,8 @@ public:

enum struct FluxCalcType {
MOENG = 0, ///< Moeng functional form
DONELAN ///< Donelan functional form
DONELAN, ///< Donelan functional form
CUSTOM ///< Custom constant flux functional form
};

enum struct ThetaCalcType {
Expand All @@ -266,10 +291,13 @@ public:
RoughCalcType rough_type;

private:
bool use_moisture;
amrex::Real z0_const;
amrex::Real surf_temp;
amrex::Real surf_heating_rate{0};
amrex::Real surf_mom_flux{0};
amrex::Real surf_temp_flux{0};
amrex::Real surf_moist_flux{0};
amrex::Real cnk_a{0.0185};
amrex::Real depth{30.0};
amrex::Real m_start_bdy_time;
Expand All @@ -280,6 +308,7 @@ private:
MOSTAverage m_ma;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> u_star;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_star;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> q_star;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> olen;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_surf;

Expand Down
29 changes: 25 additions & 4 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,11 @@ ABLMost::update_fluxes (const int& lev,
compute_fluxes(lev, max_iters, most_flux);
}
} // theta flux
} // Moeng flux
} else if (flux_type == FluxCalcType::CUSTOM) {
u_star[lev]->setVal(surf_mom_flux);
t_star[lev]->setVal(surf_temp_flux);
q_star[lev]->setVal(surf_moist_flux);
}
}


Expand All @@ -83,7 +87,7 @@ ABLMost::compute_fluxes (const int& lev,
{
// Pointers to the computed averages
const auto *const tm_ptr = m_ma.get_average(lev,2);
const auto *const umm_ptr = m_ma.get_average(lev,3);
const auto *const umm_ptr = m_ma.get_average(lev,4);

for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
{
Expand Down Expand Up @@ -121,11 +125,14 @@ ABLMost::impose_most_bcs (const int& lev,
MultiFab* z_phys)
{
const int zlo = 0;
if (flux_type == FluxCalcType::DONELAN) {
if (flux_type == FluxCalcType::MOENG) {
moeng_flux flux_comp(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
} else if (flux_type == FluxCalcType::DONELAN) {
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);
custom_flux flux_comp(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
}
}
Expand Down Expand Up @@ -177,6 +184,7 @@ ABLMost::compute_most_bcs (const int& lev,
// Get derived arrays
const auto u_star_arr = u_star[lev]->array(mfi);
const auto t_star_arr = t_star[lev]->array(mfi);
const auto q_star_arr = q_star[lev]->array(mfi);
const auto t_surf_arr = t_surf[lev]->array(mfi);

// Get LSM fluxes
Expand Down Expand Up @@ -209,6 +217,19 @@ ABLMost::compute_most_bcs (const int& lev,
}
});

// TODO: Generalize MOST q flux with MOENG & DONELAN flux types
if (flux_type == FluxCalcType::CUSTOM) {
n = RhoQ1_comp;
ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz,
cons_arr, velx_arr, vely_arr, eta_arr,
umm_arr, tm_arr, u_star_arr, q_star_arr, t_surf_arr,
dest_arr);
});
}

} else if (var_idx == Vars::xvel) {

Box xb2d = surroundingNodes(bx,0);
Expand Down
10 changes: 6 additions & 4 deletions Source/BoundaryConditions/MOSTAverage.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ public:
explicit MOSTAverage (amrex::Vector<amrex::Geometry> geom,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd);

// MOSTAverage() = default;
Expand All @@ -35,8 +36,9 @@ public:

// Reset the pointers to field MFs
void update_field_ptrs (int lev,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim);
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim);

// Compute ncells per plane
void set_plane_normalization ();
Expand Down Expand Up @@ -164,8 +166,8 @@ protected:

// General vars for multiple or all policies
//--------------------------------------------
int m_nvar{3}; // 3 fields for U/V/T
int m_navg{4}; // 4 averages for U/V/T/Umag
int m_nvar{4}; // 3 fields for U/V/T/Qv
int m_navg{5}; // 4 averages for U/V/T/Qv/Umag
int m_maxlev{0}; // Total number of levels
int m_policy{0}; // Policy for type of averaging
amrex::Real m_zref{10.0}; // Height above surface for MOST BC
Expand Down
Loading

0 comments on commit 2e78b23

Please sign in to comment.