Skip to content

Commit

Permalink
WIP debugging metgrid initialization and forcing
Browse files Browse the repository at this point in the history
  • Loading branch information
@wiersema1
wiersema1 committed Jan 23, 2024
1 parent f9f1bca commit a51636c
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 27 deletions.
8 changes: 5 additions & 3 deletions Exec/DevTests/MetGrid/inputs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 1
max_step = 10

amrex.fpe_trap_invalid = 1
amrex.fpe_trap_zero = 1
Expand Down Expand Up @@ -41,6 +41,7 @@ erf.check_int = 100 # number of timesteps between checkpoints
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 1 # number of timesteps between plotfiles
erf.plot_vars_1 = qt qv qc density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres pert_dens KE QKE
#erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres pert_dens KE QKE

# SOLVER CHOICE
erf.alpha_T = 1.0
Expand All @@ -53,14 +54,15 @@ erf.les_type = "Deardorff"
erf.Cs = 0.1

erf.moisture_model = "Kessler"
#erf.moisture_model = "Null"
#erf.do_cloud = false
#erf.do_precip = false

#erf.terrain_z_levels = 0 130 354 583 816 1054 1549 2068 2615 3193 3803 4450 5142 5892 6709 7603 8591 9702 10967 12442 14230 16610 18711 20752 22133 23960 26579 28493 31236 33699 36068 40000

# INITIALIZATION WITH ATM DATA
erf.metgrid_bdy_width = 5
erf.metgrid_bdy_set_width = 5
erf.metgrid_bdy_width = 10
erf.metgrid_bdy_set_width = 1
erf.init_type = "metgrid"
erf.nc_init_file_0 = "met_em.d01.2022-06-18_00:00:00.nc" "met_em.d01.2022-06-18_06:00:00.nc" "met_em.d01.2022-06-18_12:00:00.nc" "met_em.d01.2022-06-18_18:00:00.nc"

Expand Down
50 changes: 30 additions & 20 deletions Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -479,7 +479,10 @@ ERF::init_from_metgrid (int lev)
auto xhi_arr = bdy_data_xhi[it][ivar].array();
auto ylo_arr = bdy_data_ylo[it][ivar].array();
auto yhi_arr = bdy_data_yhi[it][ivar].array();
const Array4<Real const>& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array();
// const Array4<Real const>& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array();

// TODO: stationary boundary conditions, remove or comment when not needed for debugging.
const Array4<Real const>& fabs_for_bcs_arr = fabs_for_bcs[0][ivar].const_array();

if (ivar == MetGridBdyVars::U) {
xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag;
Expand Down Expand Up @@ -675,12 +678,12 @@ init_state_from_metgrid (const bool use_moisture,

// ---------------------------------------------------------------------------------------
// DJW: remove this after debugging is complete. This enforces quiescent conditions..
if (it == 0) {
x_vel_fab.template setVal<RunOn::Device>(0.0);
y_vel_fab.template setVal<RunOn::Device>(0.0);
}
fabs_for_bcs[it][MetGridBdyVars::U].template setVal<RunOn::Device>(0.0);
fabs_for_bcs[it][MetGridBdyVars::V].template setVal<RunOn::Device>(0.0);
// if (it == 0) {
// x_vel_fab.template setVal<RunOn::Device>(0.0);
// y_vel_fab.template setVal<RunOn::Device>(0.0);
// }
// fabs_for_bcs[it][MetGridBdyVars::U].template setVal<RunOn::Device>(0.0);
// fabs_for_bcs[it][MetGridBdyVars::V].template setVal<RunOn::Device>(0.0);
// ---------------------------------------------------------------------------------------

// ********************************************************
Expand Down Expand Up @@ -782,12 +785,12 @@ init_state_from_metgrid (const bool use_moisture,
}

// TODO: TEMPORARY CODE TO RUN QUIESCENT, REMOVE WHEN NOT NEEDED.
// if (it == 0) {
// x_vel_fab.template setVal<RunOn::Device>(0.0); // TODO: temporary code to initialize with quiescent atmosphere.
// y_vel_fab.template setVal<RunOn::Device>(0.0); // TODO: temporary code to initialize with quiescent atmosphere.
// }
// fabs_for_bcs[it][MetGridBdyVars::U].template setVal<RunOn::Device>(0.0); // TODO: temporary code to force with quiescent atmosphere.
// fabs_for_bcs[it][MetGridBdyVars::V].template setVal<RunOn::Device>(0.0); // TODO: temporary code to force with quiescent atmosphere.
if (it == 0) {
x_vel_fab.template setVal<RunOn::Device>(0.0); // TODO: temporary code to initialize with quiescent atmosphere.
y_vel_fab.template setVal<RunOn::Device>(0.0); // TODO: temporary code to initialize with quiescent atmosphere.
}
fabs_for_bcs[it][MetGridBdyVars::U].template setVal<RunOn::Device>(0.0); // TODO: temporary code to force with quiescent atmosphere.
fabs_for_bcs[it][MetGridBdyVars::V].template setVal<RunOn::Device>(0.0); // TODO: temporary code to force with quiescent atmosphere.

} // it
}
Expand Down Expand Up @@ -863,29 +866,35 @@ init_base_state_from_metgrid (const bool use_moisture,
Q_vec[k] = (use_moisture) ? new_data(i,j,k,RhoQ_comp) : 0.0;
}
z_vec[kmax+1] = new_z(i,j,kmax+1);
if ((i == 40) && (j == 40)) {
for (int k=0; k<=kmax; k++) {
Print() << " z_vec[" << k << "] = " << z_vec[k] << std::endl;
}
}

calc_rho_p(l_rdOcp,kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Theta_vec,
Q_vec,z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec);

for (int k=0; k<=kmax; k++) {
p_hse_arr(i,j,k) = Pd_vec[k];
r_hse_arr(i,j,k) = Rhom_vec[k];
r_hse_arr(i,j,k) = Rhod_vec[k];
}
});

ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k);
new_data(i,j,k,RhoScalar_comp) = 0.0;
// RhoTheta and RhoQ1 or RhoQv currently hold Theta and Q1 or Qv. Multiply by Rho.
new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k);
Real RhoTheta = r_hse_arr(i,j,k)*new_data(i,j,k,RhoTheta_comp);
new_data(i,j,k,RhoTheta_comp) = RhoTheta;
if (use_moisture){
Real RhoQ = r_hse_arr(i,j,k)*new_data(i,j,k,RhoQ_comp);
new_data(i,j,k,RhoQ_comp) = RhoQ;
}
pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
//pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp);
// p_hse_arr(i,j,k) = getPgivenRTh(RhoTheta);
//pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp);
pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp);
});
}

Expand Down Expand Up @@ -926,10 +935,11 @@ init_base_state_from_metgrid (const bool use_moisture,

for (int k=0; k<=kmax; k++) {
p_hse_arr(i,j,k) = Pd_vec[k];
// p_hse_arr(i,j,k) = getPgivenRTh(Rhod_vec[k]*Theta_vec[k]);
if (mask_c_arr(i,j,k)) {
r_hse_arr(i,j,k) = Rhom_vec[k];
if (use_moisture) Q_arr(i,j,k) = Rhom_vec[k]*Q_vec[k];
Theta_arr(i,j,k) = Rhom_vec[k]*Theta_vec[k];
r_hse_arr(i,j,k) = Rhod_vec[k];
if (use_moisture) Q_arr(i,j,k) = Rhod_vec[k]*Q_vec[k];
Theta_arr(i,j,k) = Rhod_vec[k]*Theta_vec[k];
}
} // k
});
Expand Down
14 changes: 10 additions & 4 deletions Source/Initialization/Metgrid_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,30 +122,36 @@ calc_rho_p (const amrex::Real& l_rdOcp,
}

// calculate moist density at the surface.
Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], l_rdOcp, Q_vec[0])*(1.0+Q_vec[0]);
// Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], l_rdOcp, Q_vec[0])*(1.0+Q_vec[0]);
Rhom_vec[0] = 1.0/(R_d/p_0*Theta_vec[0]*(1.0+(R_v/R_d)*Q_vec[0])*std::pow(Pm_vec[0]/p_0, -iGamma));

// integrate from the surface to the top boundary.
for (int k=1; k<=kmax; k++) {
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
amrex::Real dz = z_vec[k]-z_vec[k-1];
Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess.
for (int it=0; it<maxiter; it++) {
Pm_vec[k] = Pm_vec[k-1]-0.5*dz*(Rhom_vec[k]+Rhom_vec[k-1])*CONST_GRAV;
if (Pm_vec[k] < 0.0) Pm_vec[k] = 0.0;
Rhom_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pm_vec[k], l_rdOcp, Q_vec[k])*(1.0+Q_vec[k]);
// Rhom_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pm_vec[k], l_rdOcp, Q_vec[k])*(1.0+Q_vec[k]);
Rhom_vec[k] = 1.0/(R_d/p_0*Theta_vec[k]*qvf*std::pow(Pm_vec[k]/p_0, -iGamma));
} // it
} // k

// integrate from the top back down to get dry pressure and density.
Pd_vec[kmax] = Pm_vec[kmax];
Rhod_vec[kmax] = getRhogivenThetaPress(Theta_vec[kmax], Pd_vec[kmax], l_rdOcp, Q_vec[kmax]);
// Rhod_vec[kmax] = getRhogivenThetaPress(Theta_vec[kmax], Pd_vec[kmax], l_rdOcp, Q_vec[kmax]);
Rhod_vec[kmax] = 1.0/(R_d/p_0*Theta_vec[kmax]*std::pow(Pd_vec[kmax]/p_0, -iGamma));
for (int k=kmax-1; k>=0; k--) {
amrex::Real dz = z_vec[k+1]-z_vec[k];
Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
for (int it=0; it < maxiter; it++) {
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0;
Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], l_rdOcp, Q_vec[k]);
// Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], l_rdOcp, Q_vec[k]);
Rhod_vec[k] = 1.0/(R_d/p_0*Theta_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma));
} // it
Pd_vec[k] = getPgivenRTh(Rhod_vec[k]*Theta_vec[k]);
} // k
}

Expand Down

0 comments on commit a51636c

Please sign in to comment.