Skip to content

Commit

Permalink
FastEddy microphysics clean up
Browse files Browse the repository at this point in the history
  • Loading branch information
Mahesh Natarajan committed Dec 11, 2023
1 parent bb7f3a8 commit 56dac28
Show file tree
Hide file tree
Showing 5 changed files with 40 additions and 418 deletions.
33 changes: 0 additions & 33 deletions Source/Microphysics/FastEddy/Diagnose_FE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,37 +5,4 @@
* from the existing Microphysics variables.
*/
void FastEddy::Diagnose () {

auto qt = mic_fab_vars[MicVar_FE::qt];
auto qp = mic_fab_vars[MicVar_FE::qp];
auto qv = mic_fab_vars[MicVar_FE::qv];
auto qn = mic_fab_vars[MicVar_FE::qn];
auto qcl = mic_fab_vars[MicVar_FE::qcl];
auto qci = mic_fab_vars[MicVar_FE::qci];
auto qpl = mic_fab_vars[MicVar_FE::qpl];
auto qpi = mic_fab_vars[MicVar_FE::qpi];
auto tabs = mic_fab_vars[MicVar_FE::tabs];

for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qt_array = qt->array(mfi);
auto qp_array = qp->array(mfi);
auto qv_array = qv->array(mfi);
auto qn_array = qn->array(mfi);
auto qcl_array = qcl->array(mfi);
auto qci_array = qci->array(mfi);
auto qpl_array = qpl->array(mfi);
auto qpi_array = qpi->array(mfi);

const auto& box3d = mfi.tilebox();

amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
qv_array(i,j,k) = std::max(0.0,qt_array(i,j,k) - qn_array(i,j,k));
amrex::Real omn = 1.0;
qcl_array(i,j,k) = qn_array(i,j,k)*omn;
qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn);
amrex::Real omp = 1.0;;
qpl_array(i,j,k) = qp_array(i,j,k)*omp;
qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp);
});
}
}
15 changes: 2 additions & 13 deletions Source/Microphysics/FastEddy/FastEddy.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,22 +20,11 @@
namespace MicVar_FE {
enum {
// independent variables
qt = 0,
qp,
qv = 0,
qc,
theta, // liquid/ice water potential temperature
tabs, // temperature
rho, // density
pres, // pressure
// derived variables
qr, // rain water
qv, // water vapor
qn, // cloud condensate (liquid+ice), initial to zero
qci, // cloud ice
qcl, // cloud water
qpl, // precip rain
qpi, // precip ice
// temporary variable
omega,
NumVars
};
}
Expand Down
201 changes: 29 additions & 172 deletions Source/Microphysics/FastEddy/FastEddy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,201 +13,58 @@ using namespace amrex;
*/
void FastEddy::AdvanceFE() {

auto qt = mic_fab_vars[MicVar_FE::qt];
auto qp = mic_fab_vars[MicVar_FE::qp];
auto qn = mic_fab_vars[MicVar_FE::qn];
auto tabs = mic_fab_vars[MicVar_FE::tabs];
auto tabs = mic_fab_vars[MicVar_FE::tabs];

auto qcl = mic_fab_vars[MicVar_FE::qcl];
auto theta = mic_fab_vars[MicVar_FE::theta];
auto qv = mic_fab_vars[MicVar_FE::qv];
auto rho = mic_fab_vars[MicVar_FE::rho];

auto dz = m_geom.CellSize(2);
auto domain = m_geom.Domain();
int k_lo = domain.smallEnd(2);
int k_hi = domain.bigEnd(2);

MultiFab fz;
auto ba = tabs->boxArray();
auto dm = tabs->DistributionMap();
fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells

for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
auto fz_array = fz.array(mfi);
const Box& tbz = mfi.tilebox();

ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {

Real rho_avg, qp_avg;

if (k==k_lo) {
rho_avg = rho_array(i,j,k);
qp_avg = qp_array(i,j,k);
} else if (k==k_hi+1) {
rho_avg = rho_array(i,j,k-1);
qp_avg = qp_array(i,j,k-1);
} else {
rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
}

qp_avg = std::max(0.0, qp_avg);

Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s

fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;

/*if(k==0){
fz_array(i,j,k) = 0;
}*/
});
}

Real dtn = dt;

/*for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi);
auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi);
const auto& box3d = mfi.tilebox();
ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));
if(qt_array(i,j,k) == 0.0){
qv_array(i,j,k) = 0.0;
qn_array(i,j,k) = 0.0;
}
});
}
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);
auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
const auto& box3d = mfi.tilebox();
auto fz_array = fz.array(mfi);
ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn;
qp_array(i,j,k) = qp_array(i,j,k) + dq_sed;
});
}*/




// get the temperature, dentisy, theta, qt and qp from input
// get the temperature, dentisy, theta, qt and qc from input
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
auto qn_array = mic_fab_vars[MicVar_FE::qn]->array(mfi);
auto qt_array = mic_fab_vars[MicVar_FE::qt]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_FE::qp]->array(mfi);

auto theta_array = theta->array(mfi);
auto qv_array = qv->array(mfi);
auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto qv_array = mic_fab_vars[MicVar_FE::qv]->array(mfi);
auto qc_array = mic_fab_vars[MicVar_FE::qc]->array(mfi);
auto tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi);
auto rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);

const auto& box3d = mfi.tilebox();

auto fz_array = fz.array(mfi);

// Expose for GPU
Real d_fac_cond = m_fac_cond;

ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {

qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));

if(qt_array(i,j,k) == 0.0){
qv_array(i,j,k) = 0.0;
qn_array(i,j,k) = 0.0;
}
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));

//------- Autoconversion/accretion
Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat;
Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat;

Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0;
erf_qsatw(tabs_array(i,j,k), pressure, qsat);

// If there is precipitating water (i.e. rain), and the cell is not saturated
// then the rain water can evaporate leading to extraction of latent heat, hence
// reducing temperature and creating negative buoyancy

dq_clwater_to_rain = 0.0;
dq_rain_to_vapor = 0.0;
dq_vapor_to_clwater = 0.0;
dq_clwater_to_vapor = 0.0;

// If there is precipitating water (i.e. rain), and the cell is not saturated
// then the rain water can evaporate leading to extraction of latent heat, hence
// reducing temperature and creating negative buoyancy

//Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
dq_vapor_to_clwater = 0.0;
dq_clwater_to_vapor = 0.0;

// If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
if(qv_array(i,j,k) > qsat){
dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
}
//Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));


// If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and
// reducing temperature
if(qv_array(i,j,k) < qsat and qn_array(i,j,k) > 0.0){
dq_clwater_to_vapor = std::min(qn_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
}

if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) {
Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046);
dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/
(5.4e5 + 2.55e6/(pressure*qsat))*dtn;
// The negative sign is to make this variable (vapor formed from evaporation)
// a poistive quantity (as qv/qs < 1)
dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor});

// Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature
}

// If there is cloud water present then do accretion and autoconversion to rain

if (qn_array(i,j,k) > 0.0) {
qcc = qn_array(i,j,k);

autor = 0.0;
if (qcc > qcw0) {
autor = 0.001;
}

accrr = 0.0;
accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875);
dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - 0.001));

// If the amount of change is more than the amount of qc present, then dq = qc
dq_clwater_to_rain = std::min(dq_clwater_to_rain, qn_array(i,j,k));
// If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
if(qv_array(i,j,k) > qsat){
dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
}
// If water vapor is less than the satruated value, then the cloud water can evaporate, leading to evaporative cooling and
// reducing temperature
if(qv_array(i,j,k) < qsat and qc_array(i,j,k) > 0.0){
dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
}

qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor;
qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor;

Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn;

qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain;
qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor;
qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;

theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);
theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond*(dq_vapor_to_clwater - dq_clwater_to_vapor);

qt_array(i,j,k) = std::max(0.0, qt_array(i,j,k));
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
qn_array(i,j,k) = std::max(0.0, qn_array(i,j,k));
qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));

});
}
Expand Down
Loading

0 comments on commit 56dac28

Please sign in to comment.