Skip to content

Commit

Permalink
Modify dry compressible buoyancy to use thetapert not Tpert
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Nov 20, 2024
1 parent d7d3c29 commit 23d6ef9
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 77 deletions.
155 changes: 88 additions & 67 deletions Source/SourceTerms/ERF_buoyancy_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,55 +7,34 @@
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_anelastic (int& i,
buoyancy_dry_thetapert (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>& p0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
//
// This returns - g rho0 (thetaprime / theta)
//
// 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)/r0_arr(i,j,k-1);
amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/r0_arr(i,j,k);

amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi);

amrex::Real theta_d_0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)) / r0_arr(i,j,k-1);
amrex::Real theta_d_0_hi = getRhoThetagivenP(p0_arr(i,j,k )) / r0_arr(i,j,k );
//
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_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi);

return (-grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1)));
}
// We don't do divide by 2 here because the "2" will cancel below
amrex::Real theta_d_zface = theta_d_lo + theta_d_hi;

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_anelastic_T (int& i,
int& j,
int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp);
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;
// We don't do divide by 2 here because the "2" will cancel below
amrex::Real theta_d_0_zface = th0_arr(i,j,k) + th0_arr(i,j,k-1);

amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;
// rho0 on z-face
amrex::Real r0_zface = 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1));

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
return -grav_gpu * r0_zface * (theta_d_zface - theta_d_0_zface) / theta_d_0_zface;
}


AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
Expand All @@ -68,6 +47,9 @@ buoyancy_moist_anelastic (int& i,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
//
// This returns - g rho0 (thetaprime / theta)
//
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1);
amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp) /r0_arr(i,j,k-1);
amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp) /r0_arr(i,j,k-1);
Expand All @@ -90,41 +72,17 @@ buoyancy_moist_anelastic (int& i,
return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1)));
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_default (int& i,
int& j,
int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp));
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;

amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type1 (int& i,
int& j,
int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
buoyancy_rhopert (int& i,
int& j,
int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real rhop_hi = cell_data(i,j,k ,Rho_comp) - r0_arr(i,j,k );
amrex::Real rhop_lo = cell_data(i,j,k-1,Rho_comp) - r0_arr(i,j,k-1);
Expand Down Expand Up @@ -272,4 +230,67 @@ buoyancy_type4 (int& i,

return (-qavg * r0avg * grav_gpu);
}

// ********************************************************************************
// The functions below this are not currently in use
// ********************************************************************************

#if 0
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_Tpert_comp (int& i,
int& j,
int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
//
// This returns - g rho0 (Tprime / T)
//
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp));
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;

amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_Tpert_anel (int& i,
int& j,
int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
//
// This returns - g rho0 (Tprime / T)
//
amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp);
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;

amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
}
#endif
#endif
19 changes: 9 additions & 10 deletions Source/SourceTerms/ERF_make_buoyancy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,15 +70,14 @@ void make_buoyancy (Vector<MultiFab>& S_data,

// Base state density and pressure
const Array4<const Real>& r0_arr = r0.const_array(mfi);
const Array4<const Real>& p0_arr = p0.const_array(mfi);
const Array4<const Real>& th0_arr = th0.const_array(mfi);

if (solverChoice.moisture_type == MoistureType::None) {
ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,
buoyancy_fab(i, j, k) = buoyancy_dry_thetapert(i,j,k,
grav_gpu[2],
r0_arr,p0_arr,cell_data);
r0_arr,th0_arr,cell_data);
});
} else {
// NOTE: For decomposition in the vertical direction, klo may not
Expand Down Expand Up @@ -121,7 +120,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,

ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data);
buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data);
});
} // mfi

Expand All @@ -147,17 +146,17 @@ void make_buoyancy (Vector<MultiFab>& S_data,

// Base state density and pressure
const Array4<const Real>& r0_arr = r0.const_array(mfi);
const Array4<const Real>& p0_arr = p0.const_array(mfi);
//const Array4<const Real>& p0_arr = p0.const_array(mfi);
const Array4<const Real>& th0_arr = th0.const_array(mfi);

const Array4<const Real> & cell_data = S_data[IntVars::cons].array(mfi);
const Array4< Real> & buoyancy_fab = buoyancy.array(mfi);

ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
buoyancy_fab(i, j, k) = buoyancy_dry_default(i,j,k,
grav_gpu[2],rd_over_cp,
r0_arr,p0_arr,th0_arr,cell_data);
buoyancy_fab(i, j, k) = buoyancy_dry_thetapert(i,j,k,
grav_gpu[2],
r0_arr,th0_arr,cell_data);
});
} // mfi
} // buoyancy_type
Expand Down Expand Up @@ -196,8 +195,8 @@ void make_buoyancy (Vector<MultiFab>& S_data,

ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_qstate,
grav_gpu[2],r0_arr,cell_data);
buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate,
grav_gpu[2],r0_arr,cell_data);
});
} // mfi

Expand Down

0 comments on commit 23d6ef9

Please sign in to comment.