Skip to content

Commit

Permalink
Cloud SAM (erf-model#1632)
Browse files Browse the repository at this point in the history 
* Extend changes in SAM cloud.

* Make inline function static to avoid implicit capture issue.
  • Loading branch information
@AMLattanzi
AMLattanzi authored May 28, 2024
1 parent cef8354 commit d80cad7
Show file tree
Hide file tree
Showing 4 changed files with 171 additions and 122 deletions.
9 changes: 3 additions & 6 deletions Source/BoundaryConditions/MOSTStress.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1174,7 +1174,7 @@ struct custom_flux
const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
const amrex::Array4<amrex::Real>& dest_arr) const
{
amrex::Real rho, eta, qv;
amrex::Real rho, eta;

int ic, jc;
ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
Expand All @@ -1183,7 +1183,6 @@ struct custom_flux
jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;

rho = cons_arr(ic,jc,zlo,Rho_comp);
qv = cons_arr(ic,jc,zlo,RhoQ1_comp) / rho;

amrex::Real qstar = q_star_arr(ic,jc,zlo);
amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0;
Expand Down Expand Up @@ -1220,12 +1219,12 @@ struct custom_flux
const amrex::Array4<const amrex::Real>& /*vely_arr*/,
const amrex::Array4<const amrex::Real>& /*umm_arr*/,
const amrex::Array4<const amrex::Real>& /*tm_arr*/,
const amrex::Array4<const amrex::Real>& u_star_arr,
const amrex::Array4<const amrex::Real>& /*u_star_arr*/,
const amrex::Array4<const amrex::Real>& t_star_arr,
const amrex::Array4<const amrex::Real>& /*t_surf_arr*/,
const amrex::Array4<amrex::Real>& dest_arr) const
{
amrex::Real rho, eta, theta;
amrex::Real rho, eta;

int ic, jc;
ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i;
Expand All @@ -1234,9 +1233,7 @@ struct custom_flux
jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc;

rho = cons_arr(ic,jc,zlo,Rho_comp);
theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho;

amrex::Real ustar = u_star_arr(ic,jc,zlo);
amrex::Real tstar = t_star_arr(ic,jc,zlo);
amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar : 0.0;
amrex::Real deltaz = dz * (zlo - k);
Expand Down
162 changes: 52 additions & 110 deletions Source/Microphysics/SAM/Cloud_SAM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,6 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
Real fac_fus = m_fac_fus;
Real rdOcp = m_rdOcp;

Real tol = 1.0e-4;

SolverChoice sc = solverChoice;

int SAM_moisture_type = 1;
Expand Down Expand Up @@ -47,16 +45,12 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
// Saturation moisture fractions
Real omn, domn;
Real qsat, dqsat;
Real qsatw, dqsatw;
Real qsati, dqsati;
Real omn;
Real qsat;
Real qsatw;
Real qsati;

// Newton iteration vars
int niter;
Real fff, dfff, dtabs;
Real lstar, dlstar;
Real lstarw, lstari;
Real delta_qv, delta_qc, delta_qi;

// NOTE: Conversion before iterations is necessary to
Expand Down Expand Up @@ -118,122 +112,40 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
pres_array(i,j,k) *= 0.01;
}

// Initial guess for temperature & pressure
Real tabs = tabs_array(i,j,k);
Real pres = pres_array(i,j,k);

// Saturation moisture fractions
erf_qsatw(tabs, pres, qsatw);
erf_qsati(tabs, pres, qsati);
erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
qsat = omn * qsatw + (1.0-omn) * qsati;

// We have enough total moisture to relax to equilibrium
if (qt_array(i,j,k) > qsat) {
niter = 0;
dtabs = 1;
//==================================================
// Newton iteration to qv=qsat (cloud phase only)
//==================================================
do {
// Latent heats and their derivatives wrt to T
lstarw = fac_cond;
lstari = fac_fus;
domn = 0.0;

// Saturation moisture fractions
erf_qsatw(tabs, pres, qsatw);
erf_qsati(tabs, pres, qsati);
erf_dtqsatw(tabs, pres, dqsatw);
erf_dtqsati(tabs, pres, dqsati);

if (SAM_moisture_type == 1) {
// Cloud ice not permitted (condensation & fusion)
if(tabs >= tbgmax) {
omn = 1.0;
}
// Cloud water not permitted (sublimation & fusion)
else if(tabs <= tbgmin) {
omn = 0.0;
lstarw = fac_sub;
}
// Mixed cloud phase (condensation & fusion)
else {
omn = an*tabs-bn;
domn = an;
}
} else if (SAM_moisture_type == 2) {
omn = 1.0;
domn = 0.0;
}

// Linear combination of each component
qsat = omn * qsatw + (1.0-omn) * qsati;
dqsat = omn * dqsatw + (1.0-omn) * dqsati
+ domn * qsatw - domn * qsati;
lstar = omn * lstarw + (1.0-omn) * lstari;
dlstar = domn * lstarw - domn * lstari;

// Function for root finding:
// 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat);

// Derivative of function (T_new iterated on)
dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat;

// Update the temperature
dtabs = -fff/dfff;
tabs = tabs+dtabs;

// For now at least we perform this iteration at constant pressure
// For the moist bubble case, the results are indistinguisable
// between running with this used vs commented out
// pres = rho_array(i,j,k) * R_d * tabs
// * (1.0 + R_v/R_d * qsat) * 0.01;

// Update iteration
niter = niter+1;
} while(std::abs(dtabs) > tol && niter < 20);

// Update qsat from last iteration (dq = dq/dt * dt)
qsat += dqsat*dtabs;

// Changes in each component
delta_qv = qv_array(i,j,k) - qsat;
delta_qc = std::max(-qcl_array(i,j,k), delta_qv * omn);
delta_qi = std::max(-qci_array(i,j,k), delta_qv * (1.0-omn));

// Partition the change in non-precipitating q
qv_array(i,j,k) = qsat;
qcl_array(i,j,k) += delta_qc;
qci_array(i,j,k) += delta_qi;
qn_array(i,j,k) = qcl_array(i,j,k) + qci_array(i,j,k);
qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);

// Update temperature
tabs_array(i,j,k) = tabs;
tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type ,
fac_cond , fac_fus , fac_sub ,
an , bn ,
tabs_array, pres_array,
qv_array , qcl_array , qci_array,
qn_array , qt_array);

// Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
// This would be inconsistent with updating the pressure as part of the iteration above.
// Empirically based on the moist bubble rise case, getting the correct theta here
// depends on using the old (unchanged by the phase changes) pressure.

theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);

// Update pressure to be consistent with updated theta_array
pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));

// This was used in the earlier implmentation when we updated theta using this new pressure
// pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
// * (1.0 + R_v/R_d * qv_array(i,j,k));

// Pressure unit conversion
pres_array(i,j,k) *= 0.01;

//
// We cannot blindly relax to qsat, but we can convert qc/qi -> qv
// The concept here is that if we assume the temperature change due to this conversion
// doesn't change qsat, we can safely put all the moisture into qv without reaching qv > qsat
// because we are in the "else" part of the test on (qt_array(i,j,k) > qsat)
// We cannot blindly relax to qsat, but we can convert qc/qi -> qv.
// The concept here is that if we put all the moisture into qv and modify
// the temperature, we can then check if qv > qsat occurs (for final T/P/qv).
// If the reduction in T/qsat and increase in qv does trigger the
// aforementioned condition, we can do Newton iteration to drive qv = qsat.
//
} else {
// Changes in each component
Expand All @@ -251,15 +163,45 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
// Update temperature (endothermic since we evap/sublime)
tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi;

// Update pressure
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
// Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
// This would be inconsistent with updating the pressure as part of the iteration above.
// Empirically based on the moist bubble rise case, getting the correct theta here
// depends on using the old (unchanged by the phase changes) pressure.
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);

// Update theta from temperature
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
// Update pressure to be consistent with updated theta_array
pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));

// Pressure unit conversion
pres_array(i,j,k) *= 0.01;

// Verify assumption that qv > qsat does not occur
erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
qsat = omn * qsatw + (1.0-omn) * qsati;
if (qt_array(i,j,k) > qsat) {

// Update temperature
tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type ,
fac_cond , fac_fus , fac_sub ,
an , bn ,
tabs_array, pres_array,
qv_array , qcl_array , qci_array,
qn_array , qt_array);

// Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
// This would be inconsistent with updating the pressure as part of the iteration above.
// Empirically based on the moist bubble rise case, getting the correct theta here
// depends on using the old (unchanged by the phase changes) pressure.
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);

// Update pressure to be consistent with updated theta_array
pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));

// Pressure unit conversion
pres_array(i,j,k) *= 0.01;

}
}
});
} // mfi
Expand Down
116 changes: 116 additions & 0 deletions Source/Microphysics/SAM/SAM.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,122 @@ public:
int
Qstate_Size () override { return SAM::m_qstate_size; }

AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
static amrex::Real
NewtonIterSat (int& i, int& j, int& k,
const int& SAM_moisture_type,
const amrex::Real& fac_cond,
const amrex::Real& fac_fus,
const amrex::Real& fac_sub,
const amrex::Real& an,
const amrex::Real& bn,
const amrex::Array4<amrex::Real>& tabs_array,
const amrex::Array4<amrex::Real>& pres_array,
const amrex::Array4<amrex::Real>& qv_array,
const amrex::Array4<amrex::Real>& qc_array,
const amrex::Array4<amrex::Real>& qi_array,
const amrex::Array4<amrex::Real>& qn_array,
const amrex::Array4<amrex::Real>& qt_array)
{
// Solution tolerance
amrex::Real tol = 1.0e-4;

// Saturation moisture fractions
amrex::Real omn, domn;
amrex::Real qsat, dqsat;
amrex::Real qsatw, dqsatw;
amrex::Real qsati, dqsati;

// Newton iteration vars
int niter;
amrex::Real fff, dfff, dtabs;
amrex::Real lstar, dlstar;
amrex::Real lstarw, lstari;
amrex::Real delta_qv, delta_qc, delta_qi;

// Initial guess for temperature & pressure
amrex::Real tabs = tabs_array(i,j,k);
amrex::Real pres = pres_array(i,j,k);

niter = 0;
dtabs = 1;
//==================================================
// Newton iteration to qv=qsat (cloud phase only)
//==================================================
do {
// Latent heats and their derivatives wrt to T
lstarw = fac_cond;
lstari = fac_fus;
domn = 0.0;

// Saturation moisture fractions
erf_qsatw(tabs, pres, qsatw);
erf_qsati(tabs, pres, qsati);
erf_dtqsatw(tabs, pres, dqsatw);
erf_dtqsati(tabs, pres, dqsati);

if (SAM_moisture_type == 1) {
// Cloud ice not permitted (condensation & fusion)
if(tabs >= tbgmax) {
omn = 1.0;
}
// Cloud water not permitted (sublimation & fusion)
else if(tabs <= tbgmin) {
omn = 0.0;
lstarw = fac_sub;
}
// Mixed cloud phase (condensation & fusion)
else {
omn = an*tabs-bn;
domn = an;
}
} else if (SAM_moisture_type == 2) {
omn = 1.0;
domn = 0.0;
}

// Linear combination of each component
qsat = omn * qsatw + (1.0-omn) * qsati;
dqsat = omn * dqsatw + (1.0-omn) * dqsati
+ domn * qsatw - domn * qsati;
lstar = omn * lstarw + (1.0-omn) * lstari;
dlstar = domn * lstarw - domn * lstari;

// Function for root finding:
// 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat);

// Derivative of function (T_new iterated on)
dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat;

// Update the temperature
dtabs = -fff/dfff;
tabs += dtabs;

// Update iteration
niter = niter+1;
} while(std::abs(dtabs) > tol && niter < 20);

// Update qsat from last iteration (dq = dq/dt * dt)
qsat += dqsat*dtabs;

// Changes in each component
delta_qv = qv_array(i,j,k) - qsat;
delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn);
delta_qi = std::max(-qi_array(i,j,k), delta_qv * (1.0-omn));

// Partition the change in non-precipitating q
qv_array(i,j,k) = qsat;
qc_array(i,j,k) += delta_qc;
qi_array(i,j,k) += delta_qi;
qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);

// Return to temperature
return tabs;
}

private:
// Number of qmoist variables (qt, qv, qcl, qci, qp, qpr, qps, qpg)
int m_qmoist_size = 11;
Expand Down
6 changes: 0 additions & 6 deletions Source/TimeIntegration/ERF_slow_rhs_pre.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,12 +202,6 @@ void erf_slow_rhs_pre (int level, int finest_level,
// Define updates and fluxes in the current RK stage
// *****************************************************************************

// Open bc will be imposed upon all vars (we only access cons here for simplicity)
const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open);
const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open);
const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open);
const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open);

#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand Down

0 comments on commit d80cad7

Please sign in to comment.