Skip to content

Commit

Permalink
Input sounding initialization with moisture.
Browse files Browse the repository at this point in the history
  • Loading branch information
Aaron Lattanzi committed Dec 7, 2023
1 parent c0a1c9f commit 52c1489
Show file tree
Hide file tree
Showing 5 changed files with 118 additions and 73 deletions.
4 changes: 4 additions & 0 deletions Exec/RegTests/Bubble/input_sounding
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
1000. 300.0 1.0
0.0 300.0 1.0 0.0 0.0
10000.0 300.0 1.0 0.0 0.0

78 changes: 78 additions & 0 deletions Exec/RegTests/Bubble/inputs_kessler
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 2000
stop_time = 3600.0

amrex.fpe_trap_invalid = 1

fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
geometry.prob_extent = 20000.0 312.5 10000.0
amr.n_cell = 256 4 128
geometry.is_periodic = 1 1 0
zlo.type = "SlipWall"
zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.fixed_dt = 0.5
erf.fixed_mri_dt_ratio = 4

# 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 = 100 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhotheta rhoQt rhoQp rhoadv_0 x_velocity y_velocity z_velocity pressure theta scalar temp pres_hse dens_hse pert_pres pert_dens qt qp qv qc qi

# SOLVER CHOICES
erf.use_gravity = true
erf.use_coriolis = false
erf.use_rayleigh_damping = false

erf.buoyancy_type = 1

erf.dycore_horiz_adv_type = "Upwind_3rd"
erf.dycore_vert_adv_type = "Upwind_3rd"
erf.dryscal_horiz_adv_type = "Upwind_3rd"
erf.dryscal_vert_adv_type = "Upwind_3rd"
erf.moistscal_horiz_adv_type = "Upwind_3rd"
erf.moistscal_vert_adv_type = "Upwind_3rd"

# PHYSICS OPTIONS
erf.les_type = "None"
erf.pbl_type = "None"
erf.moisture_model = "Kessler"
erf.use_moist_background = true

erf.molec_diff_type = "ConstantAlpha"
erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
erf.dynamicViscosity = 1.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
erf.alpha_T = 1.0 # [m^2/s]
erf.alpha_C = 1.0

# INITIAL CONDITIONS
erf.init_type = "input_sounding"
erf.init_sounding_ideal = true

#erf.init_type = ""

# PROBLEM PARAMETERS (optional)
# warm bubble input
prob.T_pert = 0.0 # theta pert magnitude
prob.x_c = 10000.0
prob.z_c = 2000.0
prob.x_r = 2000.0
prob.z_r = 2000.0
prob.T_0 = 300.0
prob.Qt_0 = 0.0
prob.T_pert_is_airtemp = false # Perturb theta
9 changes: 5 additions & 4 deletions Source/Initialization/ERF_init_from_input_sounding.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,7 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx,
const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
int ktop = bx.bigEnd(2);

Real rho_k, rhoTh_k;
Real rho_k, qv_k, rhoTh_k;

// Set the density
rho_k = interpolate_1d(z_inp_sound, rho_inp_sound, z, inp_sound_size);
Expand All @@ -212,9 +212,10 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx,
state(i, j, k, RhoScalar_comp) = 0;

// Update hse quantities with values calculated from InputSoundingData.calc_rho_p()
r_hse_arr (i, j, k) = rho_k;
p_hse_arr (i, j, k) = getPgivenRTh(rhoTh_k); // NOTE WE ONLY USE THE DRY AIR PRESSURE
pi_hse_arr(i, j, k) = getExnergivenRTh(rhoTh_k, l_rdOcp);
qv_k = interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
r_hse_arr (i, j, k) = rho_k * (1.0 + qv_k);
p_hse_arr (i, j, k) = getPgivenRTh(rhoTh_k, qv_k); // NOTE WE ONLY USE THE DRY AIR PRESSURE
pi_hse_arr(i, j, k) = getExnergivenRTh(rhoTh_k, l_rdOcp, qv_k);

// Boundary treatment
if (k==0)
Expand Down
98 changes: 30 additions & 68 deletions Source/Initialization/InputSoundingData.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ public:
iss_z >> z >> theta >> qv >> U >> V;
if (z == 0) {
AMREX_ALWAYS_ASSERT(theta == theta_inp_sound_tmp[0]);
AMREX_ALWAYS_ASSERT(qv == qv_inp_sound_tmp[0]);
AMREX_ALWAYS_ASSERT(qv*0.001 == qv_inp_sound_tmp[0]);
U_inp_sound_tmp[0] = U;
V_inp_sound_tmp[0] = V;
} else {
Expand Down Expand Up @@ -141,94 +141,56 @@ public:
rhod_integ.resize(Ninp);

// evaluate surface quantities (k=0)
amrex::Real thm_surf = theta_ref_inp_sound *
(1 + (R_v/R_d-1) * qv_ref_inp_sound); // _moist_ theta == virtual potential temperature
pm_integ[0] = press_ref_inp_sound; // _total_ pressure (incl moisture)
rhom_integ[0] = 1 / (
R_d/p_0 * thm_surf * std::pow(pm_integ[0]/p_0, -iGamma)); // density of _moist_ air
amrex::Real thm_surf = theta_ref_inp_sound * (1 + (R_v/R_d) * qv_ref_inp_sound); // moist theta
pm_integ[0] = press_ref_inp_sound; // _total_ pressure (incl moisture)
rhod_integ[0] = 1 / (R_d/p_0 * thm_surf * std::pow(pm_integ[0]/p_0, -iGamma)); // density of dry air

amrex::Print() << "ideal sounding init: surface density of moist air = "
<< rhom_integ[0] << " kg/m^3" << std::endl;
amrex::Print() << "ideal sounding init: surface density of dry air = "
<< rhod_integ[0] << " kg/m^3" << std::endl;

// Note:
// p_dry = rho_d R_d T
// p_tot = rho_m R_d T_v
// = rho_d(1 + q_v) R_d T_v
// p_tot = rho_d R_d T_v
// = rho_d R_d T ( 1 + Rv/Rd * qv)

// integrate from surface to domain top
amrex::Real qvf, dz;
#if 0 // Printing
#if 1 // Printing
// In this absence of moisture, this moist profile will match the
// following dry profile
amrex::Print() << "z p_m rho_m theta U V" << std::endl;
amrex::Print() << "z p_m rho_d theta qv U V" << std::endl;
amrex::Print() << z_inp_sound[0]
<< " " << pm_integ[0]
<< " " << rhom_integ[0]
<< " " << theta_inp_sound[0]
<< " " << U_inp_sound[0]
<< " " << V_inp_sound[0]
<< std::endl;
<< " " << pm_integ[0]
<< " " << rhod_integ[0]
<< " " << theta_inp_sound[0]
<< ' ' << qv_inp_sound[0]
<< " " << U_inp_sound[0]
<< " " << V_inp_sound[0]
<< std::endl;
#endif
for (int k=1; k < size(); ++k)
{
qvf = 1 + (R_v/R_d-1) * qv_inp_sound[k];
qvf = 1 + (R_v/R_d) * qv_inp_sound[k];
dz = z_inp_sound[k] - z_inp_sound[k-1];
rhom_integ[k] = rhom_integ[k-1]; // guess
rhod_integ[k] = rhod_integ[k-1]; // guess
for (int it=0; it < maxiter; ++it)
{
pm_integ[k] = pm_integ[k-1]
- 0.5*dz*(rhom_integ[k] + rhom_integ[k-1])*CONST_GRAV;
pm_integ[k] = pm_integ[k-1] - 0.5*dz*(rhod_integ[k] + rhod_integ[k-1])*(1.0 + qv_inp_sound[k])*CONST_GRAV;
AMREX_ALWAYS_ASSERT(pm_integ[k] > 0);
rhom_integ[k] = 1 / (
R_d/p_0 * theta_inp_sound[k]*qvf * std::pow(pm_integ[k]/p_0, -iGamma));
rhod_integ[k] = 1 / (R_d/p_0 * theta_inp_sound[k]*qvf * std::pow(pm_integ[k]/p_0, -iGamma));
}
#if 0 // Printing
#if 1 // Printing
amrex::Print() << z_inp_sound[k]
<< " " << pm_integ[k]
<< " " << rhom_integ[k]
<< " " << theta_inp_sound[k]
<< " " << U_inp_sound[k]
<< " " << V_inp_sound[k]
<< std::endl;
<< " " << pm_integ[k]
<< " " << rhod_integ[k]
<< " " << theta_inp_sound[k]
<< ' ' << qv_inp_sound[k]
<< " " << U_inp_sound[k]
<< " " << V_inp_sound[k]
<< std::endl;
#endif
}

// now, integrate from the top of the sounding (where it's assumed dry) back
// down to get the dry air column properties
pd_integ[Ninp-1] = pm_integ[Ninp-1];
rhod_integ[Ninp-1] = 1 / (
R_d/p_0 * theta_inp_sound[Ninp-1] * std::pow(pd_integ[Ninp-1]/p_0, -iGamma));
amrex::Print() << "z p_d rho_d theta U V" << std::endl;
amrex::Print() << z_inp_sound[Ninp-1]
<< " " << pd_integ[Ninp-1]
<< " " << rhod_integ[Ninp-1]
<< " " << theta_inp_sound[Ninp-1]
<< " " << U_inp_sound[Ninp-1]
<< " " << V_inp_sound[Ninp-1]
<< std::endl;
for (int k=Ninp-2; k >= 0; --k)
{
dz = z_inp_sound[k+1] - z_inp_sound[k];
rhod_integ[k] = rhod_integ[k+1]; // guess
for (int it=0; it < maxiter; ++it)
{
pd_integ[k] = pd_integ[k+1]
+ 0.5*dz*(rhod_integ[k] + rhod_integ[k+1])*CONST_GRAV;
AMREX_ALWAYS_ASSERT(pd_integ[k] > 0);
rhod_integ[k] = 1 / (
R_d/p_0 * theta_inp_sound[k] * std::pow(pd_integ[k]/p_0, -iGamma));
}
amrex::Print() << z_inp_sound[k]
<< " " << pd_integ[k]
<< " " << rhod_integ[k]
<< " " << theta_inp_sound[k]
<< " " << U_inp_sound[k]
<< " " << V_inp_sound[k]
<< std::endl;
}
// Note: at this point, the surface pressure, density of the dry air
// column is stored in pd_integ[0], rhod_integ[0]

// update
host_to_device();
}
Expand Down Expand Up @@ -266,7 +228,7 @@ public:
rhod_integ.begin(), rhod_integ.end(),
rho_inp_sound_d.begin());
amrex::Gpu::copy(amrex::Gpu::hostToDevice,
pd_integ.begin(), pd_integ.end(),
pm_integ.begin(), pm_integ.end(),
p_inp_sound_d.begin());
}
}
Expand Down
2 changes: 1 addition & 1 deletion Source/Microphysics/Kessler/Init_Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@ void Kessler::Init (const MultiFab& cons_in, MultiFab& qmoist,
});

// This fills qv
//Diagnose();
Diagnose();

#if 0
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int k, int j, int i) {
Expand Down

0 comments on commit 52c1489

Please sign in to comment.