ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Private Types | Private Attributes | Static Private Attributes | List of all members
Kessler Class Reference

#include <ERF_Kessler.H>

Inheritance diagram for Kessler:
Inheritance graph
[legend]
Collaboration diagram for Kessler:
Collaboration graph
[legend]

Public Member Functions

 Kessler ()
 
virtual ~Kessler ()=default
 
void AdvanceKessler (const SolverChoice &solverChoice)
 
void Define (SolverChoice &sc) override
 
void Init (const amrex::MultiFab &cons_in, const amrex::BoxArray &grids, const amrex::Geometry &geom, const amrex::Real &dt_advance, std::unique_ptr< amrex::MultiFab > &z_phys_nd, std::unique_ptr< amrex::MultiFab > &detJ_cc) 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
 
- 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
 

Private Types

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

Private Attributes

int m_qmoist_size = 1
 
int n_qstate_moist_size = 3
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
amrex::BoxArray m_gtoe
 
amrex::Real dt
 
int nlev
 
int zlo
 
int zhi
 
int m_axis
 
amrex::Real m_fac_cond
 
amrex::Real m_fac_fus
 
amrex::Real m_fac_sub
 
amrex::Real m_gOcp
 
bool m_do_cond
 
amrex::MultiFab * m_z_phys_nd
 
amrex::MultiFab * m_detJ_cc
 
amrex::Array< FabPtr, MicVar_Kess::NumVars > mic_fab_vars
 

Static Private Attributes

static constexpr amrex::Real CFL_MAX = 0.5
 

Member Typedef Documentation

◆ FabPtr

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

Constructor & Destructor Documentation

◆ Kessler()

Kessler::Kessler ( )
inline
44 {}

◆ ~Kessler()

virtual Kessler::~Kessler ( )
virtualdefault

Member Function Documentation

◆ Advance()

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

Reimplemented from NullMoist.

99  {
100  dt = dt_advance;
101 
102  this->AdvanceKessler(solverChoice);
103  }
Kessler::dt
amrex::Real dt
Definition: ERF_Kessler.H:150
Kessler::AdvanceKessler
void AdvanceKessler(const SolverChoice &solverChoice)
Definition: ERF_Kessler.cpp:11
Here is the call graph for this function:

◆ AdvanceKessler()

void Kessler::AdvanceKessler ( const SolverChoice &  solverChoice)

Compute Precipitation-related Microphysics quantities.

12 {
13  bool do_cond = m_do_cond;
14  auto tabs = mic_fab_vars[MicVar_Kess::tabs];
15  if (solverChoice.moisture_type == MoistureType::Kessler){
16  auto dz = m_geom.CellSize(2);
17  auto domain = m_geom.Domain();
18  int k_lo = domain.smallEnd(2);
19  int k_hi = domain.bigEnd(2);
20 
21  MultiFab fz;
22  auto ba = tabs->boxArray();
23  auto dm = tabs->DistributionMap();
24  fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells
25 
26  Real dtn = dt;
27  Real coef = dtn/dz;
28 
29  // Saturation and evaporation calculations go first
30  for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
31  auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
32  auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
33  auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
34  auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi);
35  auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi);
36  auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi);
37  auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi);
38  auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
39 
40  const auto& tbx = mfi.tilebox();
41 
42  // Expose for GPU
43  Real d_fac_cond = m_fac_cond;
44 
45  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
46  {
47  qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
48  qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
49  qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
50 
51  //------- Autoconversion/accretion
52  Real qcc, auto_r, accrr;
53  Real qsat, dtqsat;
54  Real dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater;
55 
56  Real pressure = pres_array(i,j,k);
57  erf_qsatw(tabs_array(i,j,k), pressure, qsat);
58  erf_dtqsatw(tabs_array(i,j,k), pressure, dtqsat);
59 
60  if (qsat <= 0.0) {
61  amrex::Warning("qsat computed as non-positive; setting to 0.!");
62  qsat = 0.0;
63  }
64 
65  // If there is precipitating water (i.e. rain), and the cell is not saturated
66  // then the rain water can evaporate leading to extraction of latent heat, hence
67  // reducing temperature and creating negative buoyancy
68 
69  dq_clwater_to_rain = 0.0;
70  dq_rain_to_vapor = 0.0;
71  dq_vapor_to_clwater = 0.0;
72  dq_clwater_to_vapor = 0.0;
73 
74  //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
75  //Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
76  Real fac = (L_v/Cp_d)*dtqsat;
77 
78  // If water vapor content exceeds saturation value, then vapor condenses to water and latent heat is released, increasing temperature
79  if ( (qv_array(i,j,k) > qsat) && do_cond ) {
80  dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
81  }
82 
83  // If water vapor is less than the saturated value, then the cloud water can evaporate,
84  // leading to evaporative cooling and reducing temperature
85  if ( (qv_array(i,j,k) < qsat) && (qc_array(i,j,k) > 0.0) && do_cond ) {
86  dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
87  }
88 
89  if (( qp_array(i,j,k) > 0.0) && (qv_array(i,j,k) < qsat) ) {
90  Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046);
91  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)/
92  (5.4e5 + 2.55e6/(pressure*qsat))*dtn;
93  // The negative sign is to make this variable (vapor formed from evaporation)
94  // a positive quantity (as qv/qs < 1)
95  dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor});
96 
97  // Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature
98  }
99 
100  // If there is cloud water present then do accretion and autoconversion to rain
101  if (qc_array(i,j,k) > 0.0) {
102  qcc = qc_array(i,j,k);
103 
104  auto_r = 0.0;
105  if (qcc > qcw0) {
106  auto_r = alphaelq;
107  }
108 
109  accrr = 0.0;
110  accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875);
111  dq_clwater_to_rain = dtn *(accrr*qcc + auto_r*(qcc - qcw0));
112 
113  // If the amount of change is more than the amount of qc present, then dq = qc
114  dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k));
115  }
116 
117  qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor;
118  qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
119  qp_array(i,j,k) += dq_clwater_to_rain - dq_rain_to_vapor;
120 
121  Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);
122  theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);
123 
124  qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
125  qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
126  qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
127 
128  qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);
129  });
130  }
131 
132  // Precompute terminal velocity for substepping
133  for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
134  auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
135  auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
136 
137  auto fz_array = fz.array(mfi);
138  const Box& tbz = mfi.tilebox();
139 
140  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
141  {
142  Real rho_avg, qp_avg;
143 
144  if (k==k_lo) {
145  rho_avg = rho_array(i,j,k);
146  qp_avg = qp_array(i,j,k);
147  } else if (k==k_hi+1) {
148  rho_avg = rho_array(i,j,k-1);
149  qp_avg = qp_array(i,j,k-1);
150  } else {
151  rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
152  qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
153  }
154 
155  qp_avg = std::max(0.0, qp_avg);
156 
157  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
158 
159  // NOTE: Fz is the sedimentation flux from the advective operator.
160  // In the terrain-following coordinate system, the z-deriv in
161  // the divergence uses the normal velocity (Omega). However,
162  // there are no u/v components to the sedimentation velocity.
163  // Therefore, we simply end up with a division by detJ when
164  // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
165  fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
166  });
167  } // mfi
168 
169  // Compute number of substeps from maximum terminal velocity
170  Real wt_max;
171  int n_substep;
172  auto const& ma_fz_arr = fz.const_arrays();
173  GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
174  TypeList<Real>{},
175  fz, IntVect(0),
176  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
177  -> GpuTuple<Real>
178  {
179  return { ma_fz_arr[box_no](i,j,k) };
180  });
181  wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
182  n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
183  AMREX_ALWAYS_ASSERT(n_substep >= 1);
184  coef /= Real(n_substep);
185  dtn /= Real(n_substep);
186 
187  // Substep the vertical advection
188  for (int nsub(0); nsub<n_substep; ++nsub) {
189  for ( MFIter mfi(*tabs, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
190  auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
191  auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
192  auto rain_accum_array = mic_fab_vars[MicVar_Kess::rain_accum]->array(mfi);
193  auto fz_array = fz.array(mfi);
194 
195  const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};
196 
197  const Box& tbx = mfi.tilebox();
198  const Box& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));
199 
200  // Update vertical flux every substep
201  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
202  {
203  Real rho_avg, qp_avg;
204  if (k==k_lo) {
205  rho_avg = rho_array(i,j,k);
206  qp_avg = qp_array(i,j,k);
207  } else if (k==k_hi+1) {
208  rho_avg = rho_array(i,j,k-1);
209  qp_avg = qp_array(i,j,k-1);
210  } else {
211  rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
212  qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
213  }
214 
215  qp_avg = std::max(0.0, qp_avg);
216 
217  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
218 
219  // NOTE: Fz is the sedimentation flux from the advective operator.
220  // In the terrain-following coordinate system, the z-deriv in
221  // the divergence uses the normal velocity (Omega). However,
222  // there are no u/v components to the sedimentation velocity.
223  // Therefore, we simply end up with a division by detJ when
224  // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
225  fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
226 
227  if(k==k_lo){
228  rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*qp_avg*V_terminal*dtn/1000.0*1000.0; // Divide by rho_water and convert to mm
229  }
230  });
231 
232  // Update precip every substep
233  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
234  {
235  // Jacobian determinant
236  Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;
237 
238  if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0;
239  if(std::fabs(fz_array(i,j,k )) < 1e-14) fz_array(i,j,k ) = 0.0;
240  Real dq_sed = dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k)) * coef;
241  if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0;
242 
243  qp_array(i,j,k) += dq_sed;
244  qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
245  });
246  } // mfi
247  } // nsub
248  }
249 
250 
251  if (solverChoice.moisture_type == MoistureType::Kessler_NoRain) {
252  if (!do_cond) { return; }
253  // get the temperature, density, theta, qt and qc from input
254  for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
255  auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
256  auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
257  auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi);
258  auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi);
259  auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi);
260  auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi);
261 
262  const auto& box3d = mfi.tilebox();
263 
264  // Expose for GPU
265  Real d_fac_cond = m_fac_cond;
266 
267  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
268  {
269  qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
270 
271  //------- Autoconversion/accretion
272  Real qsat, dtqsat;
273  Real dq_clwater_to_vapor, dq_vapor_to_clwater;
274 
275  Real pressure = pres_array(i,j,k);
276  erf_qsatw(tabs_array(i,j,k), pressure, qsat);
277  erf_dtqsatw(tabs_array(i,j,k), pressure, dtqsat);
278 
279  // If there is precipitating water (i.e. rain), and the cell is not saturated
280  // then the rain water can evaporate leading to extraction of latent heat, hence
281  // reducing temperature and creating negative buoyancy
282 
283  dq_vapor_to_clwater = 0.0;
284  dq_clwater_to_vapor = 0.0;
285 
286  //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2));
287  //Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));
288  Real fac = (L_v/Cp_d)*dtqsat;
289 
290  // If water vapor content exceeds saturation value, then vapor condenses to water and latent heat is released, increasing temperature
291  if (qv_array(i,j,k) > qsat){
292  dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
293  }
294  // If water vapor is less than the saturated value, then the cloud water can evaporate, leading to evaporative cooling and
295  // reducing temperature
296  if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){
297  dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
298  }
299 
300  qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor;
301  qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor;
302 
303  Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);
304 
305  theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor);
306 
307  qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
308  qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
309 
310  qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);
311  });
312  }
313  }
314 }
Cp_d
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
alphaelq
constexpr amrex::Real alphaelq
Definition: ERF_Constants.H:48
qcw0
constexpr amrex::Real qcw0
Definition: ERF_Constants.H:46
L_v
constexpr amrex::Real L_v
Definition: ERF_Constants.H:16
erf_dtqsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw(amrex::Real t, amrex::Real p, amrex::Real &dtqsatw)
Definition: ERF_MicrophysicsUtils.H:178
erf_qsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:163
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Kessler::m_geom
amrex::Geometry m_geom
Definition: ERF_Kessler.H:144
Kessler::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_Kessler.H:167
Kessler::m_fac_cond
amrex::Real m_fac_cond
Definition: ERF_Kessler.H:159
Kessler::m_do_cond
bool m_do_cond
Definition: ERF_Kessler.H:163
Kessler::mic_fab_vars
amrex::Array< FabPtr, MicVar_Kess::NumVars > mic_fab_vars
Definition: ERF_Kessler.H:170
Kessler::CFL_MAX
static constexpr amrex::Real CFL_MAX
Definition: ERF_Kessler.H:138
MicVar_Kess::qp
@ qp
Definition: ERF_Kessler.H:31
MicVar_Kess::qcl
@ qcl
Definition: ERF_Kessler.H:29
MicVar_Kess::tabs
@ tabs
Definition: ERF_Kessler.H:24
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:25
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:22
MicVar_Kess::theta
@ theta
Definition: ERF_Kessler.H:23
MicVar_Kess::qt
@ qt
Definition: ERF_Kessler.H:27
MicVar_Kess::rain_accum
@ rain_accum
Definition: ERF_Kessler.H:33
MicVar_Kess::qv
@ qv
Definition: ERF_Kessler.H:28
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:911

Referenced by Advance().

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

◆ Copy_Micro_to_State()

void Kessler::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, qi, qr, qs, qg

Reimplemented from NullMoist.

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

Referenced by Update_State_Vars().

Here is the caller graph for this function:

◆ Copy_State_to_Micro()

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

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input

Reimplemented from NullMoist.

65 {
66  // Get the temperature, density, theta, qt and qp from input
67  for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
68  const auto& box3d = mfi.tilebox();
69 
70  auto states_array = cons_in.array(mfi);
71 
72  auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi);
73  auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
74  auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
75 
76  auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
77 
78  auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
79  auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi);
80  auto tabs_array = mic_fab_vars[MicVar_Kess::tabs]->array(mfi);
81  auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi);
82 
83  // Get pressure, theta, temperature, density, and qt, qp
84  ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
85  {
86  rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
87  theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
88  qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
89  qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
90  qp_array(i,j,k) = states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp);
91  qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);
92 
93  tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
94  states_array(i,j,k,RhoTheta_comp),
95  qv_array(i,j,k));
96  pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * 0.01;
97  });
98  }
99 }
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:81
getTgivenRandRTh
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
Rho_comp
#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 Kessler::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

55  {
56  m_fac_cond = lcond / sc.c_p;
57  m_fac_fus = lfus / sc.c_p;
58  m_fac_sub = lsub / sc.c_p;
59  m_gOcp = CONST_GRAV / sc.c_p;
60  m_axis = sc.ave_plane;
61  m_do_cond = (!sc.use_shoc);
62  }
lfus
constexpr amrex::Real lfus
Definition: ERF_Constants.H:67
lsub
constexpr amrex::Real lsub
Definition: ERF_Constants.H:68
lcond
constexpr amrex::Real lcond
Definition: ERF_Constants.H:66
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
Kessler::m_fac_sub
amrex::Real m_fac_sub
Definition: ERF_Kessler.H:161
Kessler::m_gOcp
amrex::Real m_gOcp
Definition: ERF_Kessler.H:162
Kessler::m_fac_fus
amrex::Real m_fac_fus
Definition: ERF_Kessler.H:160
Kessler::m_axis
int m_axis
Definition: ERF_Kessler.H:156
SolverChoice::c_p
amrex::Real c_p
Definition: ERF_DataStruct.H:865
SolverChoice::use_shoc
bool use_shoc
Definition: ERF_DataStruct.H:895
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:926

◆ Init()

void Kessler::Init ( const amrex::MultiFab &  cons_in,
const amrex::BoxArray &  grids,
const amrex::Geometry &  geom,
const amrex::Real &  dt_advance,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
std::unique_ptr< amrex::MultiFab > &  detJ_cc 
)
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.

28 {
29  dt = dt_advance;
30  m_geom = geom;
31  m_gtoe = grids;
32 
33  m_z_phys_nd = z_phys_nd.get();
34  m_detJ_cc = detJ_cc.get();
35 
36  MicVarMap.resize(m_qmoist_size);
37  MicVarMap = {MicVar_Kess::rain_accum};
38 
39  // initialize microphysics variables
40  for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) {
41  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
42  1, cons_in.nGrowVect());
43  mic_fab_vars[ivar]->setVal(0.);
44  }
45 
46  // Set class data members
47  for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
48  const auto& box3d = mfi.tilebox();
49 
50  const auto& lo = lbound(box3d);
51  const auto& hi = ubound(box3d);
52 
53  nlev = box3d.length(2);
54  zlo = lo.z;
55  zhi = hi.z;
56  }
57 }
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
Kessler::zlo
int zlo
Definition: ERF_Kessler.H:153
Kessler::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_Kessler.H:141
Kessler::zhi
int zhi
Definition: ERF_Kessler.H:153
Kessler::nlev
int nlev
Definition: ERF_Kessler.H:153
Kessler::m_z_phys_nd
amrex::MultiFab * m_z_phys_nd
Definition: ERF_Kessler.H:166
Kessler::m_gtoe
amrex::BoxArray m_gtoe
Definition: ERF_Kessler.H:147
Kessler::m_qmoist_size
int m_qmoist_size
Definition: ERF_Kessler.H:132
MicVar_Kess::NumVars
@ NumVars
Definition: ERF_Kessler.H:34
Here is the call graph for this function:

◆ Qmoist_Ptr()

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

Reimplemented from NullMoist.

107  {
108  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
109  return mic_fab_vars[MicVarMap[varIdx]].get();
110  }

◆ Qmoist_Restart_Vars()

void Kessler::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 
126  // NOTE: These are the indices to access into qmoist (not mic_fab_vars)
127  a_idx.push_back(0); a_names.push_back("RainAccum");
128  }

◆ Qmoist_Size()

int Kessler::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

113 { return Kessler::m_qmoist_size; }

◆ Qstate_Moist_Size()

int Kessler::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

116 { return Kessler::n_qstate_moist_size; }
Kessler::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_Kessler.H:135

◆ Update_Micro_Vars()

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

Reimplemented from NullMoist.

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

◆ Update_State_Vars()

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

Reimplemented from NullMoist.

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

Member Data Documentation

◆ CFL_MAX

constexpr amrex::Real Kessler::CFL_MAX = 0.5
staticconstexprprivate

◆ dt

amrex::Real Kessler::dt
private

Referenced by Advance().

◆ m_axis

int Kessler::m_axis
private

Referenced by Define().

◆ m_detJ_cc

amrex::MultiFab* Kessler::m_detJ_cc
private

◆ m_do_cond

bool Kessler::m_do_cond
private

Referenced by Define().

◆ m_fac_cond

amrex::Real Kessler::m_fac_cond
private

Referenced by Define().

◆ m_fac_fus

amrex::Real Kessler::m_fac_fus
private

Referenced by Define().

◆ m_fac_sub

amrex::Real Kessler::m_fac_sub
private

Referenced by Define().

◆ m_geom

amrex::Geometry Kessler::m_geom
private

◆ m_gOcp

amrex::Real Kessler::m_gOcp
private

Referenced by Define().

◆ m_gtoe

amrex::BoxArray Kessler::m_gtoe
private

◆ m_qmoist_size

int Kessler::m_qmoist_size = 1
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_z_phys_nd

amrex::MultiFab* Kessler::m_z_phys_nd
private

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar_Kess::NumVars> Kessler::mic_fab_vars
private

Referenced by Qmoist_Ptr().

◆ MicVarMap

amrex::Vector<int> Kessler::MicVarMap
private

Referenced by Qmoist_Ptr().

◆ n_qstate_moist_size

int Kessler::n_qstate_moist_size = 3
private

Referenced by Qstate_Moist_Size().

◆ nlev

int Kessler::nlev
private

◆ zhi

int Kessler::zhi
private

◆ zlo

int Kessler::zlo
private
The documentation for this class was generated from the following files: