Skip to content

Commit

Permalink
Fix index error in last commit and fix plotfile.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Dec 19, 2023
1 parent 1f9ffee commit 252d7e2
Show file tree
Hide file tree
Showing 9 changed files with 47 additions and 665 deletions.
54 changes: 21 additions & 33 deletions Source/IO/Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -605,78 +605,66 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
mf_comp ++;
}

// TODO: The introduction of Q3 changes these plots
// Should we make the micro model own qt and qp
// since it knows whether qt = qv + qc or qt = qv + qc + qi?
//-----------------------------------------------------------
// NOTE: Protect against accessing non-existent data
if (use_moisture) {
int q_size = qmoist[lev].size();
MultiFab qv_mf(*(qmoist[lev][0]), make_alias, 0, 1);

if (containerHasElement(plot_var_names, "qt") && (q_size >= 3))
if (containerHasElement(plot_var_names, "qt") && (q_size >= 1))
{
MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1);
MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0);
MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0);
MultiFab qt_mf(*(qmoist[lev][0]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qt_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qp") && (q_size >= 6))
if (containerHasElement(plot_var_names, "qv") && (q_size >= 2))
{
MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1);
MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1);
MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0);
MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0);
MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0);
MultiFab qv_mf(*(qmoist[lev][1]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qv"))
if (containerHasElement(plot_var_names, "qc") && (q_size >= 3))
{
MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
MultiFab qc_mf(*(qmoist[lev][2]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qc") && (q_size >= 2))
if (containerHasElement(plot_var_names, "qi") && (q_size >= 4))
{
MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0);
MultiFab qi_mf(*(qmoist[lev][3]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qi") && (q_size >= 3))
if (containerHasElement(plot_var_names, "qp") && (q_size >= 5))
{
MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0);
MultiFab qp_mf(*(qmoist[lev][4]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qp_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4))
if (containerHasElement(plot_var_names, "qrain") && (q_size >= 6))
{
MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1);
MultiFab qr_mf(*(qmoist[lev][5]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 5))
if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 7))
{
MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1);
MultiFab qs_mf(*(qmoist[lev][6]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0);
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 6))
if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 8))
{
MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1);
MultiFab qg_mf(*(qmoist[lev][7]), make_alias, 0, 1);
MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0);
mf_comp += 1;
}
}
//-----------------------------------------------------------

#ifdef ERF_USE_PARTICLES
if (containerHasElement(plot_var_names, "tracer_particle_count"))
Expand Down
8 changes: 5 additions & 3 deletions Source/Microphysics/FastEddy/FastEddy.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,11 @@ namespace MicVar_FE {
// independent variables
qv = 0,
qc,
qt,
rho, // density
theta, // liquid/ice water potential temperature
tabs, // temperature
rho, // density
pres, // pressure
NumVars
};
}
Expand Down Expand Up @@ -114,8 +116,8 @@ public:
Qstate_Size () { return FastEddy::m_qstate_size; }

private:
// Number of qmoist variables
int m_qmoist_size = 2;
// Number of qmoist variables (qt, qv, qc)
int m_qmoist_size = 3;

// Number of qstate variables
int m_qstate_size = 2;
Expand Down
2 changes: 1 addition & 1 deletion Source/Microphysics/FastEddy/FastEddy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ void FastEddy::AdvanceFE ()
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);
auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi);

const auto& box3d = mfi.tilebox();

Expand Down
98 changes: 9 additions & 89 deletions Source/Microphysics/FastEddy/Init_FE.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

#include <AMReX_GpuContainers.H>
#include "Microphysics.H"
#include "IndexDefines.H"
Expand Down Expand Up @@ -29,7 +28,7 @@ void FastEddy::Init (const MultiFab& cons_in,
m_gtoe = grids;

MicVarMap.resize(m_qmoist_size);
MicVarMap = {MicVar_FE::qv, MicVar_FE::qc};
MicVarMap = {MicVar_FE::qt, MicVar_FE::qv, MicVar_FE::qc};

// initialize microphysics variables
for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) {
Expand Down Expand Up @@ -65,13 +64,14 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in)

auto states_array = cons_in.array(mfi);

auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
auto qt_array = mic_fab_vars[MicVar_FE::qt]->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_Kess::rho]->array(mfi);
auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi);
auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi);
auto pres_array = mic_fab_vars[MicVar::pres]->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 tabs_array = mic_fab_vars[MicVar_FE::tabs]->array(mfi);
auto pres_array = mic_fab_vars[MicVar_FE::pres]->array(mfi);

// Get pressure, theta, temperature, density, and qt, qp
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
Expand All @@ -80,6 +80,7 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in)
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);
qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);

tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
states_array(i,j,k,RhoTheta_comp),
Expand All @@ -88,84 +89,3 @@ void FastEddy::Copy_State_to_Micro (const MultiFab& cons_in)
});
}
}



#if 0
#include <AMReX_GpuContainers.H>
#include "Microphysics.H"
#include "IndexDefines.H"
#include "PlaneAverage.H"
#include "EOS.H"
#include "TileNoZ.H"

using namespace amrex;

/**
* Initializes the Microphysics module.
*
* @param[in] cons_in Conserved variables input
* @param[in] qc_in Cloud variables input
* @param[in,out] qv_in Vapor variables input
* @param[in] qi_in Ice variables input
* @param[in] grids The boxes on which we will evolve the solution
* @param[in] geom Geometry associated with these MultiFabs and grids
* @param[in] dt_advance Timestep for the advance
*/
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;

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

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

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

const auto& box3d = mfi.tilebox();

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

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

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

// 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));
});
}
}
#endif
2 changes: 1 addition & 1 deletion Source/Microphysics/Kessler/Diagnose_Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ void Kessler::Diagnose ()
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;;
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);
});
Expand Down
Loading

0 comments on commit 252d7e2

Please sign in to comment.