Skip to content

Commit

Permalink
Merge develop and adhere to formatting.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Dec 11, 2023
1 parent 772dbc2 commit 2587dab
Show file tree
Hide file tree
Showing 3 changed files with 102 additions and 99 deletions.
83 changes: 42 additions & 41 deletions Source/Microphysics/FastEddy/FastEddy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,61 +11,62 @@ using namespace amrex;
/**
* Compute Precipitation-related Microphysics quantities.
*/
void FastEddy::AdvanceFE() {
void FastEddy::AdvanceFE () {

auto tabs = mic_fab_vars[MicVar_FE::tabs];
auto tabs = mic_fab_vars[MicVar_FE::tabs];

// get the temperature, dentisy, theta, qt and qc from input
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++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);
// get the temperature, dentisy, theta, qt and qc from input
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++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();
const auto& box3d = mfi.tilebox();

// Expose for GPU
Real d_fac_cond = m_fac_cond;
// Expose for GPU
Real d_fac_cond = m_fac_cond;

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

qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));

//------- Autoconversion/accretion
Real dq_clwater_to_vapor, dq_vapor_to_clwater, qsat;
//------- Autoconversion/accretion
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);
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
// 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_vapor_to_clwater = 0.0;
dq_clwater_to_vapor = 0.0;
dq_vapor_to_clwater = 0.0;
dq_clwater_to_vapor = 0.0;

//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));
//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 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));
}
// 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;
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;

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);
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);

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));
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));

});
}
});
}
}
75 changes: 38 additions & 37 deletions Source/Microphysics/FastEddy/Init_FE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,54 +23,55 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
const BoxArray& grids,
const Geometry& geom,
const Real& dt_advance)
{
m_geom = geom;
m_gtoe = grids;
{
m_geom = geom;
m_gtoe = grids;

dt = dt_advance;
dt = dt_advance;

// initialize microphysics variables
for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) {
mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect());
mic_fab_vars[ivar]->setVal(0.);
}
// initialize microphysics variables
for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) {
mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect());
mic_fab_vars[ivar]->setVal(0.);
}

// We must initialize to zero since we now need boundary values for the call to getP and we need all values filled
// The ghost cells of these arrays aren't filled in the boundary condition calls for the state
// We must initialize to zero since we now need boundary values for the call to getP and we need all values filled
// The ghost cells of these arrays aren't filled in the boundary condition calls for the state

qmoist.setVal(0.);
qmoist.setVal(0.);

for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {

const auto& box3d = mfi.tilebox();
const auto& box3d = mfi.tilebox();

const auto& lo = amrex::lbound(box3d);
const auto& hi = amrex::ubound(box3d);
const auto& lo = amrex::lbound(box3d);
const auto& hi = amrex::ubound(box3d);

nlev = box3d.length(2);
zlo = lo.z;
zhi = hi.z;
nlev = box3d.length(2);
zlo = lo.z;
zhi = hi.z;
}

// Get the temperature, density, theta, qt and qp from input
for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) {
auto states_array = cons_in.array(mfi);
// Get the temperature, density, theta, qt and qp from input
for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) {
auto states_array = cons_in.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 rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi);
auto temp_array = mic_fab_vars[MicVar_FE::tabs]->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 rho_array = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto theta_array = mic_fab_vars[MicVar_FE::theta]->array(mfi);
auto temp_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);

const auto& box3d = mfi.tilebox();
const auto& box3d = mfi.tilebox();

// Get pressure, theta, temperature, density, and qt, qp
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k));
});
}
// Get pressure, theta, temperature, density, and qt, qp
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k));
});
}
}
43 changes: 22 additions & 21 deletions Source/Microphysics/FastEddy/Update_FE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,27 +14,28 @@ void FastEddy::Update (amrex::MultiFab& cons,
amrex::MultiFab& qmoist)
{

// Get the temperature, density, theta, qt and qp from input
for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto states_arr = cons.array(mfi);

auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi);
auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi);
auto qc_arr = mic_fab_vars[MicVar_FE::qc]->array(mfi);
const auto& box3d = mfi.tilebox();

// get potential total density, temperature, qt, qp
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k);
states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k);
});
}

// Fill interior ghost cells and periodic boundaries
cons.FillBoundary(m_geom.periodicity());
qmoist.FillBoundary(m_geom.periodicity());
// Get the temperature, density, theta, qt and qp from input
for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto states_arr = cons.array(mfi);

auto rho_arr = mic_fab_vars[MicVar_FE::rho]->array(mfi);
auto theta_arr = mic_fab_vars[MicVar_FE::theta]->array(mfi);
auto qv_arr = mic_fab_vars[MicVar_FE::qv]->array(mfi);
auto qc_arr = mic_fab_vars[MicVar_FE::qc]->array(mfi);
const auto& box3d = mfi.tilebox();

// get potential total density, temperature, qt, qp
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k);
states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k);
});
}

// Fill interior ghost cells and periodic boundaries
cons.FillBoundary(m_geom.periodicity());
qmoist.FillBoundary(m_geom.periodicity());
}


0 comments on commit 2587dab

Please sign in to comment.