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

Buoyancy Fix & WRFInput Rebalance #2162

Open
wants to merge 1 commit into
base: development
Choose a base branch
from
+176 −60
Open

Buoyancy Fix & WRFInput Rebalance #2162

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
8 changes: 8 additions & 0 deletions Source/DataStructs/ERF_DataStruct.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,6 +228,11 @@ struct SolverChoice {
}
}

// Check for rebalancing with wrfinput
if (init_type == InitType::WRFInput) {
pp.query("rebalance_wrfinput",rebalance_wrfinput);
}


if (terrain_type == TerrainType::StaticFittedMesh ||
terrain_type == TerrainType::MovingFittedMesh) {
Expand Down Expand Up @@ -745,6 +750,9 @@ struct SolverChoice {
// Monotonic advection limiter
bool use_mono_adv{false};

// Rebalance wrfinput
bool rebalance_wrfinput{false};

CouplingType coupling_type;
MoistureType moisture_type;
WindFarmType windfarm_type;
Expand Down
123 changes: 117 additions & 6 deletions Source/Initialization/ERF_InitFromWRFInput.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,8 @@ init_base_state_from_wrfinput (const Box& domain,
MultiFab& qv_hse,
MultiFab& r_hse,
const MultiFab& mf_PB,
const MultiFab& mf_P);
const MultiFab& mf_P,
const bool& use_P_eos);

/**
* ERF function that initializes data from a WRF dataset
Expand Down Expand Up @@ -516,6 +517,112 @@ ERF::init_from_wrfinput (int lev)
make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
}

// **************************************************************************
// Rebalance the base state if needed
// **************************************************************************
if (solverChoice.rebalance_wrfinput) {
int ncomp = lev_new[Vars::cons].nComp();
int k_dom_lo = geom[lev].Domain().smallEnd(2);
int k_dom_hi = geom[lev].Domain().bigEnd(2);
Real tol = 1.0e-10;
Real grav = CONST_GRAV;
for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
Box bx = mfi.tilebox();
int klo = bx.smallEnd(2);
int khi = bx.bigEnd(2);
AMREX_ALWAYS_ASSERT((klo == k_dom_lo) && (khi == k_dom_hi));
bx.makeSlab(2,klo);

const Array4< Real>& con_arr = lev_new[Vars::cons].array(mfi);
const Array4<const Real>& z_arr = z_phys_nd[lev]->const_array(mfi);

ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
{
// integrate from surface to domain top
Real Factor;
Real dz, F, C;
Real rho_tot_hi, rho_tot_lo;
Real z_lo, z_hi;
Real R_lo, R_hi;
Real qv_lo, qv_hi;
Real Th_lo, Th_hi;
Real P_lo, P_hi;

// First integrate from sea level to the height at klo
{
// Vertical grid spacing
z_lo = 0.0;
z_hi = 0.125 * (z_arr(i,j,klo ) + z_arr(i+1,j,klo ) + z_arr(1,j+1,klo ) + z_arr(i+1,j+1,klo )
+z_arr(i,j,klo+1) + z_arr(i+1,j,klo+1) + z_arr(1,j+1,klo+1) + z_arr(i+1,j+1,klo+1));
dz = z_hi - z_lo;

// Establish known constant
qv_lo = con_arr(i,j,klo,RhoQ1_comp) / con_arr(i,j,klo,Rho_comp);
Th_lo = con_arr(i,j,klo,RhoTheta_comp) / con_arr(i,j,klo,Rho_comp);
P_lo = p_0;
R_lo = getRhogivenThetaPress(Th_lo, P_lo, R_d/Cp_d, qv_lo);
rho_tot_lo = R_lo * (1. + qv_lo);
C = -P_lo + 0.5*rho_tot_lo*grav*dz;

// Initial guess and residual
qv_hi = con_arr(i,j,klo,RhoQ1_comp) / con_arr(i,j,klo,Rho_comp);
Th_hi = con_arr(i,j,klo,RhoTheta_comp) / con_arr(i,j,klo,Rho_comp);
P_hi = p_0;
R_hi = getRhogivenThetaPress(Th_hi, P_hi, R_d/Cp_d, qv_hi);
rho_tot_hi = R_hi * (1. + qv_hi);
F = P_hi + 0.5*rho_tot_hi*grav*dz + C;

// Do iterations
HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz,
grav, C, Th_hi,
qv_hi, qv_hi,
P_hi, R_hi, F);

// Assign data
Factor = R_hi / con_arr(i,j,klo,Rho_comp);
con_arr(i,j,klo,Rho_comp) = R_hi;
for (int n(1); n<ncomp; ++n) { con_arr(i,j,klo,n) *= Factor; }
P_lo = P_hi;
z_lo = z_hi;
}

for (int k(klo+1); k<=khi; ++k) {
// Vertical grid spacing
z_hi = 0.125 * (z_arr(i,j,k ) + z_arr(i+1,j,k ) + z_arr(1,j+1,k ) + z_arr(i+1,j+1,k )
+z_arr(i,j,k+1) + z_arr(i+1,j,k+1) + z_arr(1,j+1,k+1) + z_arr(i+1,j+1,k+1));
dz = z_hi - z_lo;

// Establish known constant
qv_lo = con_arr(i,j,k,RhoQ1_comp) / con_arr(i,j,k,Rho_comp);
Th_lo = con_arr(i,j,k,RhoTheta_comp) / con_arr(i,j,k,Rho_comp);
R_lo = getRhogivenThetaPress(Th_lo, P_lo, R_d/Cp_d, qv_lo);
rho_tot_lo = R_lo * (1. + qv_lo);
C = -P_lo + 0.5*rho_tot_lo*grav*dz;

// Initial guess and residual
qv_hi = con_arr(i,j,k,RhoQ1_comp) / con_arr(i,j,k,Rho_comp);
Th_hi = con_arr(i,j,k,RhoTheta_comp) / con_arr(i,j,k,Rho_comp);
R_hi = getRhogivenThetaPress(Th_hi, P_hi, R_d/Cp_d, qv_hi);
rho_tot_hi = R_hi * (1. + qv_hi);
F = P_hi + 0.5*rho_tot_hi*grav*dz + C;

// Do iterations
HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz,
grav, C, Th_hi,
qv_hi, qv_hi,
P_hi, R_hi, F);

// Assign data
Factor = R_hi / con_arr(i,j,k,Rho_comp);
con_arr(i,j,k,Rho_comp) = R_hi;
for (int n(1); n<ncomp; ++n) { con_arr(i,j,k,n) *= Factor; }
P_lo = P_hi;
z_lo = z_hi;
}
});
} // mfi
} // rebalance_wrfinput

// **************************************************************************
// Initialize the base state
// **************************************************************************
Expand All @@ -530,9 +637,11 @@ ERF::init_from_wrfinput (int lev)

int n_qstate = micro->Get_Qstate_Size();

bool use_P_eos = (solverChoice.rebalance_wrfinput);

init_base_state_from_wrfinput(domain, l_rdOcp, solverChoice.moisture_type, n_qstate,
lev_new[Vars::cons], p_hse, pi_hse, th_hse, qv_hse, r_hse,
mf_PB, mf_P);
mf_PB, mf_P, use_P_eos);

// FillBoundary to populate the internal ghost cells (no averaging in above call)
r_hse.FillBoundary(geom[lev].periodicity());
Expand Down Expand Up @@ -623,7 +732,8 @@ init_base_state_from_wrfinput (const Box& domain,
MultiFab& qv_hse,
MultiFab& r_hse,
const MultiFab& mf_PB,
const MultiFab& mf_P)
const MultiFab& mf_P,
const bool& use_P_eos)
{
const auto& dom_lo = lbound(domain);
const auto& dom_hi = ubound(domain);
Expand Down Expand Up @@ -660,6 +770,7 @@ init_base_state_from_wrfinput (const Box& domain,
cons_arr(ii,jj,kk,RhoQ1_comp) / cons_arr(ii,jj,kk,Rho_comp) : 0.0;
Real RT = cons_arr(ii,jj,kk,RhoTheta_comp);
Real P_eos = getPgivenRTh(RT, Qv);
if (use_P_eos) { Ptot = P_eos; }
Real DelP = std::fabs(Ptot - P_eos);

// NOTE: Ghost cells don't contain valid data
Expand Down Expand Up @@ -755,7 +866,7 @@ compute_terrain_top_and_bottom (Real& terrain_bottom_min,
//
// This loop computes the min and max values of the bottom surface
//
ParallelFor(Fab2dBox_lo, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(Fab2dBox_lo, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
{
int ii = std::max(std::min(i,ihi-1),ilo+1);
int jj = std::max(std::min(j,jhi-1),jlo+1);
Expand All @@ -770,7 +881,7 @@ compute_terrain_top_and_bottom (Real& terrain_bottom_min,
//
// This loop computes the max value of the top surface
//
ParallelFor(Fab2dBox_hi, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(Fab2dBox_hi, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
{
int ii = std::max(std::min(i,ihi-1),ilo+1);
int jj = std::max(std::min(j,jhi-1),jlo+1);
Expand All @@ -784,7 +895,7 @@ compute_terrain_top_and_bottom (Real& terrain_bottom_min,
//
// This loop computes the max value of the layer just below the top surface
//
ParallelFor(Fab2dBox_hi_m1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(Fab2dBox_hi_m1, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/) noexcept
{
int ii = std::max(std::min(i,ihi-1),ilo+1);
int jj = std::max(std::min(j,jhi-1),jlo+1);
Expand Down
95 changes: 45 additions & 50 deletions Source/SourceTerms/ERF_BuoyancyUtils.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,8 @@ buoyancy_dry_anelastic (int& i, int& j, int& k,
const amrex::Array4<const amrex::Real>& cell_data)
{
// Note: this is the same term as the moist anelastic buoyancy when qv = qc = qt = 0
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp);
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);

amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi);
amrex::Real theta_d0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1));
Expand All @@ -32,29 +32,35 @@ buoyancy_moist_anelastic (int& i, int& j, int& k,
amrex::Real const& rv_over_rd,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& qv0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp) /cell_data(i,j,k-1,Rho_comp);
amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp) /cell_data(i,j,k-1,Rho_comp);
amrex::Real qt_lo = qv_lo + qc_lo;
amrex::Real theta_v_lo = theta_d_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo);
amrex::Real Fact = (1.0 - rv_over_rd); // = -0.61

amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/cell_data(i,j,k,Rho_comp);
amrex::Real qv_hi = cell_data(i,j,k,RhoQ1_comp) /cell_data(i,j,k,Rho_comp);
amrex::Real qc_hi = cell_data(i,j,k,RhoQ2_comp) /cell_data(i,j,k,Rho_comp);
amrex::Real qt_hi = qv_hi + qc_hi;
amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi);

// NOTE: The following has three errors that depend upon storing a qv0 base state
// 1. theta_v0_wface needs to incorporate qv0
// 2. rho0_wface should be rhod0_wface (rho0 includes (1 + qv0) from balance w/ moisture)
// 3. another term with qv' = (qv - qv0) should be included
amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi);
amrex::Real theta_v0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1));
amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));
amrex::Real theta_v_hi = theta_d_hi * (1.0 - Fact*qt_hi - rv_over_rd*qc_hi);

amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp) /cell_data(i,j,k-1,Rho_comp);
amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp) /cell_data(i,j,k-1,Rho_comp);
amrex::Real qt_lo = qv_lo + qc_lo;
amrex::Real theta_v_lo = theta_d_lo * (1.0 - Fact*qt_lo - rv_over_rd*qc_lo);

amrex::Real theta_v0_hi = th0_arr(i,j,k ) * (1.0 - Fact*qv0_arr(i,j,k ));
amrex::Real theta_v0_lo = th0_arr(i,j,k-1) * (1.0 - Fact*qv0_arr(i,j,k-1));
amrex::Real qv0_hi = qv0_arr(i,j,k );
amrex::Real qv0_lo = qv0_arr(i,j,k-1);

amrex::Real q_hi = (qv_hi-qv0_hi) - qc_hi + (theta_v_hi - theta_v0_hi) / theta_v0_hi;
amrex::Real q_lo = (qv_lo-qv0_lo) - qc_lo + (theta_v_lo - theta_v0_lo) / theta_v0_lo;
amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo);

amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));

return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v0_wface) / theta_v0_wface);
return (-r0avg * grav_gpu * qavg);
}

AMREX_GPU_DEVICE
Expand Down Expand Up @@ -110,10 +116,6 @@ buoyancy_dry_Thpert (int& i, int& j, int& k,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_prim)
{
//
// TODO: we are currently using Theta_prime/Theta_0 to replace T_prime/T_0 - P_prime/P_0
// but I don't think that is quite right.
//
amrex::Real thetaprime_hi = (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,k );
amrex::Real thetaprime_lo = (cell_prim(i,j,k-1,PrimTheta_comp) - th0_arr(i,j,k-1)) / th0_arr(i,j,k-1);

Expand All @@ -132,33 +134,32 @@ buoyancy_moist_Tpert (int& i, int& j, int& k,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& qv0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<const amrex::Real>& cell_data)
{
//
// Note: this currently assumes the base state qv0 is identically zero
// TODO: generalize this to allow for moist base state
//
amrex::Real qv_hi = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;

amrex::Real qc_hi = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;

amrex::Real qp_hi = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;
amrex::Real qp_lo = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;

amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp);
amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp);
amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;
amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;
amrex::Real qp_lo = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;

amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp), qv_hi);
amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp), qv_lo);

amrex::Real q_hi = 0.61 * qv_hi - (qc_hi + qp_hi) + (t_hi-t0_hi)/t0_hi;
amrex::Real q_lo = 0.61 * qv_lo - (qc_lo + qp_lo) + (t_lo-t0_lo)/t0_lo;
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k ), th0_arr(i,j,k ), rd_over_cp);
amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp);

amrex::Real qv0_hi = qv0_arr(i,j,k );
amrex::Real qv0_lo = qv0_arr(i,j,k-1);

amrex::Real q_hi = amrex::Real(0.61) * (qv_hi-qv0_hi) - (qc_hi + qp_hi) + (t_hi-t0_hi)/t0_hi;
amrex::Real q_lo = amrex::Real(0.61) * (qv_lo-qv0_lo) - (qc_lo + qp_lo) + (t_lo-t0_lo)/t0_lo;
amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo);

amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo);
amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));

return ( -r0avg * grav_gpu * qavg);
Expand All @@ -172,32 +173,26 @@ buoyancy_moist_Thpert (int& i, int& j, int& k,
amrex::Real const& grav_gpu,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& qv0_arr,
const amrex::Array4<const amrex::Real>& cell_prim)
{
//
// Note: this currently assumes the base state qv0 is identically zero
// TODO: generalize this to allow for moist base state
//
amrex::Real qv_hi = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;

amrex::Real qc_hi = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;

amrex::Real qp_hi = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;

amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;
amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;
amrex::Real qp_lo = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;

//
// TODO: we are currently using Theta_prime/Theta_0 to replace T_prime/T_0 - P_prime/P_0
// but I don't think that is quite right.
//
amrex::Real q_hi = amrex::Real(0.61) * qv_hi - (qc_hi + qp_hi)
+ (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,k );
amrex::Real qv0_hi = qv0_arr(i,j,k );
amrex::Real qv0_lo = qv0_arr(i,j,k-1);

amrex::Real q_lo = amrex::Real(0.61) * qv_lo - (qc_lo + qp_lo)
amrex::Real q_hi = amrex::Real(0.61) * (qv_hi-qv0_hi) - (qc_hi + qp_hi)
+ (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,k );
amrex::Real q_lo = amrex::Real(0.61) * (qv_lo-qv0_lo) - (qc_lo + qp_lo)
+ (cell_prim(i,j,k-1,PrimTheta_comp) - th0_arr(i,j,k-1)) / th0_arr(i,j,k-1);

amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo);

amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));

return ( -r0avg * grav_gpu * qavg);
Expand Down
Loading
Loading