ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
SatAdj Class Reference

#include <ERF_SatAdj.H>

Inheritance diagram for SatAdj:
Collaboration diagram for SatAdj:

Public Member Functions

 SatAdj ()
 
virtual ~SatAdj ()=default
 
void AdvanceSatAdj (const SolverChoice &)
 
void Define (SolverChoice &sc) override
 
void Init (const amrex::MultiFab &cons_in, const amrex::BoxArray &, const amrex::Geometry &geom, const amrex::Real &dt_advance, std::unique_ptr< amrex::MultiFab > &, std::unique_ptr< amrex::MultiFab > &) override
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in) override
 
void Copy_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &solverChoice) override
 
amrex::MultiFab * Qmoist_Ptr (const int &varIdx) override
 
int Qmoist_Size () override
 
int Qstate_Moist_Size () override
 
void Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
 
void Set_RealWidth (const int real_width) override
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual int Qstate_NonMoist_Size ()
 
virtual void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const
 
virtual void GetPlotVar (const std::string &, amrex::MultiFab &) const
 
virtual void Set_dzmin (const amrex::Real)
 

Static Public Member Functions

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real NewtonIterSat (int &i, int &j, int &k, const amrex::Real &fac_cond, 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)
 

Private Types

using FabPtr = std::shared_ptr< amrex::MultiFab >
 

Private Attributes

int m_qmoist_size = 0
 
int n_qstate_moist_size = 2
 
amrex::Geometry m_geom
 
int m_real_width {0}
 
amrex::Real dt
 
amrex::Real m_fac_cond
 
amrex::Real m_rdOcp
 
bool m_do_cond
 
amrex::Array< FabPtr, MicVar_SatAdj::NumVarsmic_fab_vars
 

Member Typedef Documentation

◆ FabPtr

using SatAdj::FabPtr = std::shared_ptr<amrex::MultiFab>
private

Constructor & Destructor Documentation

◆ SatAdj()

SatAdj::SatAdj ( )
inline
47 {}

◆ ~SatAdj()

virtual SatAdj::~SatAdj ( )
virtualdefault

Member Function Documentation

◆ Advance()

void SatAdj::Advance ( const amrex::Real dt_advance,
const SolverChoice solverChoice 
)
inlineoverridevirtual

Reimplemented from NullMoist.

99  {
100  dt = dt_advance;
101 
102  this->AdvanceSatAdj(solverChoice);
103  }
void AdvanceSatAdj(const SolverChoice &)
Definition: ERF_SatAdj.cpp:8
amrex::Real dt
Definition: ERF_SatAdj.H:209
Here is the call graph for this function:

◆ AdvanceSatAdj()

void SatAdj::AdvanceSatAdj ( const SolverChoice )

Compute Precipitation-related Microphysics quantities.

9 {
10  if (!m_do_cond) { return; }
11 
13 
14  // Expose for GPU
15  Real d_fac_cond = m_fac_cond;
16  Real rdOcp = m_rdOcp;
17 
18  auto domain = m_geom.Domain();
19  int i_lo = domain.smallEnd(0);
20  int i_hi = domain.bigEnd(0);
21  int j_lo = domain.smallEnd(1);
22  int j_hi = domain.bigEnd(1);
23 
24  // get the temperature, density, theta, qt and qc from input
25  for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
26 
27  auto tbx = mfi.tilebox();
28  if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); }
29  if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); }
30  if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); }
31  if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); }
32 
33  auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi);
34  auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi);
35  auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi);
36  auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi);
37  auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi);
38 
39  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
40  {
41  //------- Evaporation/condensation
42  Real qsat;
43  erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat);
44 
45  // There is enough moisture to drive to equilibrium
46  if ((qv_array(i,j,k)+qc_array(i,j,k)) > qsat) {
47  Real qvprev = qv_array(i,j,k);
48  Real qcprev = qc_array(i,j,k);
49 
50  // clip qc but maintain total water
51  if (qc_array(i,j,k) < 0) {
52  qv_array(i,j,k) += qc_array(i,j,k);
53  qc_array(i,j,k) = 0.0;
54  }
55 
56  // Update temperature
57  tabs_array(i,j,k) = NewtonIterSat(i, j, k ,
58  d_fac_cond, tabs_array, pres_array,
59  qv_array , qc_array );
60 
61  Real qsatnew;
62  erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatnew);
63 
64  AMREX_ASSERT(std::abs(qv_array(i,j,k)-qsatnew) < 1e-12);
65 
66  amrex::ignore_unused(qvprev);
67  amrex::ignore_unused(qcprev);
68  AMREX_ASSERT(std::abs(qv_array(i,j,k)+qc_array(i,j,k)-qvprev-qcprev) < 1e-14);
69 
70  // Update theta (constant pressure)
71  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
72 
73  //
74  // We cannot blindly relax to qsat, but we can convert qc/qi -> qv.
75  // The concept here is that if we put all the moisture into qv and modify
76  // the temperature, we can then check if qv > qsat occurs (for final T/P/qv).
77  // If the reduction in T/qsat and increase in qv does trigger the
78  // aforementioned condition, we can do Newton iteration to drive qv = qsat.
79  //
80  } else {
81  // Changes in each component
82  Real delta_qc = qc_array(i,j,k);
83 
84  // Partition the change in non-precipitating q
85  qv_array(i,j,k) += qc_array(i,j,k);
86  qc_array(i,j,k) = 0.0;
87 
88  // Update temperature (endothermic since we evap/sublime)
89  tabs_array(i,j,k) -= d_fac_cond * delta_qc;
90 
91  // Update theta
92  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
93 
94  // Verify assumption that qv > qsat does not occur
95  erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat);
96  if (qv_array(i,j,k) > qsat) {
97  Real qvprev = qv_array(i,j,k);
98  Real qcprev = qc_array(i,j,k);
99  Real Tprev = tabs_array(i,j,k);
100 
101  // Update temperature
102  tabs_array(i,j,k) = NewtonIterSat(i, j, k ,
103  d_fac_cond , tabs_array, pres_array,
104  qv_array , qc_array );
105 
106  Real qsatnew;
107  erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatnew);
108  amrex::ignore_unused(qvprev);
109  amrex::ignore_unused(qcprev);
110  amrex::ignore_unused(Tprev);
111  AMREX_ASSERT(qv_array(i,j,k) < qvprev);
112  AMREX_ASSERT(qc_array(i,j,k) > qcprev);
113  AMREX_ASSERT(tabs_array(i,j,k) > Tprev);
114  AMREX_ASSERT(std::abs(qv_array(i,j,k)-qsatnew) < 1e-14);
115  AMREX_ASSERT(std::abs(qv_array(i,j,k)+qc_array(i,j,k)-qvprev-qcprev) < 1e-14);
116 
117  // Update theta
118  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
119 
120  }
121  }
122  });
123  }
124 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenTandP(const amrex::Real T, const amrex::Real P, const amrex::Real rdOcp)
Definition: ERF_EOS.H:18
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:166
amrex::Real Real
Definition: ERF_ShocInterface.H:19
amrex::Real m_fac_cond
Definition: ERF_SatAdj.H:212
amrex::Geometry m_geom
Definition: ERF_SatAdj.H:203
amrex::Array< FabPtr, MicVar_SatAdj::NumVars > mic_fab_vars
Definition: ERF_SatAdj.H:217
amrex::Real m_rdOcp
Definition: ERF_SatAdj.H:213
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real NewtonIterSat(int &i, int &j, int &k, const amrex::Real &fac_cond, 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)
Definition: ERF_SatAdj.H:133
bool m_do_cond
Definition: ERF_SatAdj.H:214
int m_real_width
Definition: ERF_SatAdj.H:206
@ tabs
Definition: ERF_Kessler.H:24
@ theta
Definition: ERF_SatAdj.H:31
@ pres
Definition: ERF_SatAdj.H:33
@ qv
Definition: ERF_SatAdj.H:35
@ tabs
Definition: ERF_SatAdj.H:32
@ qc
Definition: ERF_SatAdj.H:36

Referenced by Advance().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Copy_Micro_to_State()

void SatAdj::Copy_Micro_to_State ( amrex::MultiFab &  cons_in)
overridevirtual

Updates conserved and microphysics variables in the provided MultiFabs from the internal MultiFabs that store Microphysics module data.

Parameters
[out]consConserved variables
[out]qmoistqv, qc

Reimplemented from NullMoist.

13 {
14  // Get the temperature, density, theta, qt and qp from input
15  for (MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
16  const auto& tbx = mfi.tilebox();
17 
18  auto states_arr = cons.array(mfi);
19 
20  auto rho_arr = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi);
21  auto theta_arr = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi);
22  auto qv_arr = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi);
23  auto qc_arr = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi);
24 
25  // get potential total density, temperature, qt, qp
26  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
27  {
28  states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
29  states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k);
30  states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k);
31  });
32  }
33 
34  // Fill interior ghost cells and periodic boundaries
35  cons.FillBoundary(m_geom.periodicity());
36 }
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
@ rho
Definition: ERF_SatAdj.H:30
@ cons
Definition: ERF_IndexDefines.H:140

Referenced by Update_State_Vars().

Here is the caller graph for this function:

◆ Copy_State_to_Micro()

void SatAdj::Copy_State_to_Micro ( const amrex::MultiFab &  cons_in)
overridevirtual

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input

Reimplemented from NullMoist.

41 {
42  // Get the temperature, density, theta, qt and qp from input
43  for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
44  const auto& tbx = mfi.tilebox();
45 
46  auto states_array = cons_in.array(mfi);
47 
48  auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi);
49  auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi);
50 
51  auto rho_array = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi);
52  auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi);
53  auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi);
54  auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi);
55 
56  // Get pressure, theta, temperature, density, and qt, qp
57  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
58  {
59  rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
60  theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
61  qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
62  qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
63 
64  tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
65  states_array(i,j,k,RhoTheta_comp),
66  qv_array(i,j,k));
67 
68  // Pressure in [mbar] for qsat evaluation
69  pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * 0.01;
70  });
71  }
72 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:81
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:46
#define Rho_comp
Definition: ERF_IndexDefines.H:36

Referenced by Update_Micro_Vars().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Define()

void SatAdj::Define ( SolverChoice sc)
inlineoverridevirtual

Reimplemented from NullMoist.

58  {
59  m_fac_cond = lcond / sc.c_p;
60  m_rdOcp = sc.rdOcp;
61  m_do_cond = (!sc.use_shoc);
62  }
constexpr amrex::Real lcond
Definition: ERF_Constants.H:66
amrex::Real rdOcp
Definition: ERF_DataStruct.H:1124
amrex::Real c_p
Definition: ERF_DataStruct.H:1123
bool use_shoc
Definition: ERF_DataStruct.H:1155

◆ Init()

void SatAdj::Init ( const amrex::MultiFab &  cons_in,
const amrex::BoxArray &  ,
const amrex::Geometry &  geom,
const amrex::Real dt_advance,
std::unique_ptr< amrex::MultiFab > &  ,
std::unique_ptr< amrex::MultiFab > &   
)
overridevirtual

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input
[in]qc_inCloud variables input
[in,out]qv_inVapor variables input
[in]qi_inIce variables input
[in]gridsThe boxes on which we will evolve the solution
[in]geomGeometry associated with these MultiFabs and grids
[in]dt_advanceTimestep for the advance

Reimplemented from NullMoist.

23 {
24  dt = dt_advance;
25  m_geom = geom;
26 
27  // initialize microphysics variables
28  for (auto ivar = 0; ivar < MicVar_SatAdj::NumVars; ++ivar) {
29  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
30  1, cons_in.nGrowVect());
31  mic_fab_vars[ivar]->setVal(0.);
32  }
33 }
@ NumVars
Definition: ERF_SatAdj.H:37

◆ NewtonIterSat()

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real SatAdj::NewtonIterSat ( int &  i,
int &  j,
int &  k,
const amrex::Real fac_cond,
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 
)
inlinestatic
139  {
140  // Solution tolerance
141  amrex::Real tol = 1.0e-8;
142 
143  // Saturation moisture fractions
144  amrex::Real qsat, dqsat;
145 
146  // Newton iteration vars
147  int niter;
148  amrex::Real fff, dfff;
149  amrex::Real dtabs, delta_qv;
150 
151  // Initial guess for temperature & pressure
152  amrex::Real tabs = tabs_array(i,j,k);
153  amrex::Real pres = pres_array(i,j,k);
154 
155  niter = 0;
156  dtabs = 1;
157 
158  //==================================================
159  // Newton iteration to qv=qsat (cloud phase only)
160  //==================================================
161  do {
162  // Saturation moisture fractions
163  erf_qsatw(tabs, pres, qsat);
164  erf_dtqsatw(tabs, pres, dqsat);
165 
166  // Function for root finding:
167  // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
168  fff = -tabs + tabs_array(i,j,k) + fac_cond*(qv_array(i,j,k) - qsat);
169 
170  // Derivative of function (T_new iterated on)
171  dfff = -1.0 - fac_cond*dqsat;
172 
173  // Update the temperature
174  dtabs = -fff/dfff;
175  tabs += dtabs;
176 
177  // Update iteration
178  niter = niter+1;
179  } while(std::abs(dtabs) > tol && niter < 20);
180 
181  // Update qsat from last iteration (dq = dq/dt * dt)
182  qsat += dqsat*dtabs;
183 
184  // Changes in each component
185  delta_qv = qv_array(i,j,k) - qsat;
186 
187  // Partition the change in non-precipitating q
188  qv_array(i,j,k) = qsat;
189  qc_array(i,j,k) += delta_qv;
190 
191  // Return to temperature
192  return tabs;
193  }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw(amrex::Real t, amrex::Real p, amrex::Real &dtqsatw)
Definition: ERF_MicrophysicsUtils.H:181
Here is the call graph for this function:

◆ Qmoist_Ptr()

amrex::MultiFab* SatAdj::Qmoist_Ptr ( const int &  varIdx)
inlineoverridevirtual

Reimplemented from NullMoist.

107  {
108  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
109  return nullptr;
110  }
int m_qmoist_size
Definition: ERF_SatAdj.H:197

◆ Qmoist_Restart_Vars()

void SatAdj::Qmoist_Restart_Vars ( const SolverChoice ,
std::vector< int > &  a_idx,
std::vector< std::string > &  a_names 
) const
inlineoverridevirtual

Reimplemented from NullMoist.

122  {
123  a_idx.clear();
124  a_names.clear();
125  }

◆ Qmoist_Size()

int SatAdj::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

113 { return SatAdj::m_qmoist_size; }

◆ Qstate_Moist_Size()

int SatAdj::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

116 { return SatAdj::n_qstate_moist_size; }
int n_qstate_moist_size
Definition: ERF_SatAdj.H:200

◆ Set_RealWidth()

void SatAdj::Set_RealWidth ( const int  real_width)
inlineoverridevirtual

Reimplemented from NullMoist.

128 { m_real_width = real_width; }

◆ Update_Micro_Vars()

void SatAdj::Update_Micro_Vars ( amrex::MultiFab &  cons_in)
inlineoverridevirtual

Reimplemented from NullMoist.

84  {
85  this->Copy_State_to_Micro(cons_in);
86  }
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitSatAdj.cpp:40
Here is the call graph for this function:

◆ Update_State_Vars()

void SatAdj::Update_State_Vars ( amrex::MultiFab &  cons_in)
inlineoverridevirtual

Reimplemented from NullMoist.

91  {
92  this->Copy_Micro_to_State(cons_in);
93  }
void Copy_Micro_to_State(amrex::MultiFab &cons_in) override
Definition: ERF_UpdateSatAdj.cpp:12
Here is the call graph for this function:

Member Data Documentation

◆ dt

amrex::Real SatAdj::dt
private

Referenced by Advance().

◆ m_do_cond

bool SatAdj::m_do_cond
private

Referenced by Define().

◆ m_fac_cond

amrex::Real SatAdj::m_fac_cond
private

Referenced by Define().

◆ m_geom

amrex::Geometry SatAdj::m_geom
private

◆ m_qmoist_size

int SatAdj::m_qmoist_size = 0
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_rdOcp

amrex::Real SatAdj::m_rdOcp
private

Referenced by Define().

◆ m_real_width

int SatAdj::m_real_width {0}
private

Referenced by Set_RealWidth().

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar_SatAdj::NumVars> SatAdj::mic_fab_vars
private

◆ n_qstate_moist_size

int SatAdj::n_qstate_moist_size = 2
private

Referenced by Qstate_Moist_Size().


The documentation for this class was generated from the following files: