Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Input sounding initialization with moisture. #1333

Closed
wants to merge 2 commits into from
Closed

Input sounding initialization with moisture. #1333

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
52c1489
Input sounding initialization with moisture.
Dec 7, 2023
cb0525c
Fix without moisture.
Dec 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions Exec/RegTests/Bubble/input_sounding
View file
Open in desktop
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
Open in desktop
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
16 changes: 9 additions & 7 deletions Source/Initialization/ERF_init_from_input_sounding.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,9 +180,7 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx,
const Real* z_inp_sound = inputSoundingData.z_inp_sound_d.dataPtr();
const Real* rho_inp_sound = inputSoundingData.rho_inp_sound_d.dataPtr();
const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d.dataPtr();
#if defined(ERF_USE_MOISTURE)
const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d.dataPtr();
#elif defined(ERF_USE_WARM_NO_PRECIP)
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d.dataPtr();
#endif
const int inp_sound_size = inputSoundingData.size();
Expand All @@ -198,7 +196,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 +210,13 @@ 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 = 0.0;
#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP)
qv_k = interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
#endif
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
Open in desktop
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
Open in desktop
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
Loading