ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_SuperDropletPCMassChange.H
Go to the documentation of this file.
1 #ifndef SUPERDROPLET_PC_MASSCHANGE_H_
2 #define SUPERDROPLET_PC_MASSCHANGE_H_
3 
4 #include <cmath>
5 #include <AMReX_GpuPrint.H>
6 #include "ERF_Constants.H"
7 
8 #ifdef ERF_USE_PARTICLES
9 
10 /*! \brief Namespace with classes and functions for condensation/evaporation
11  * (liquid <--> vapour) */
12 namespace SDMassChangeUtils_LV
13 {
14  /*! \brief Ventilation factor as a function of droplet radius (Bayley et al., 2025, Eq. 3) */
15  template <typename RT /*!< real-type */ >
16  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
17  RT ventilationFactor ( const RT a_R, /*!< radius */
18  const RT a_alpha1, /*!< fit coefficient alpha_1 */
19  const RT a_beta1, /*!< fit exponent beta_1 */
20  const RT a_alpha2, /*!< fit coefficient alpha_2 */
21  const RT a_beta2, /*!< fit exponent beta_2 */
22  const RT a_fcap /*!< maximum ventilation factor */ ) noexcept
23  {
24  RT x1 = a_alpha1 * std::pow(a_R, a_beta1);
25  RT x2 = a_alpha2 * std::pow(a_R, a_beta2);
26  RT fv = RT(1.0) + (x1*x2) / (x1 + x2);
27  return std::min(fv, a_fcap);
28  }
29 
30 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
31  /*! \brief Phase change equation (in terms of R) */
32  template <typename RT /*!< real-type */ >
33  struct dRdt
34  {
35  RT L; /*!< latent heat of vaporization (condensate) */
36  RT K; /*!< thermal conductivity (condensate) */
37  RT Rv; /*!< gas constant of air with vapour */
38  RT rho_l; /*!< density of condensate */
39 
40  /*! \brief Right-hand-side of the phase change ODE */
41  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
42  RT operator () ( const RT a_R, /*!< radius */
43  const RT a_S, /*!< saturation ratio */
44  const RT a_T, /*!< temperature */
45  const RT a_e_s, /*!< saturation pressure */
46  const RT a_D, /*!< mol. diffusion coeff */
47  const RT a_a, /*!< curvature coeff */
48  const RT a_b, /*!< solute coeff */
49  const RT a_N_s, /*!< total solute moles */
50  const RT a_f_v /*!< ventilation factor */ ) const noexcept
51  {
52  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
53  RT Kn = lambda_v/a_R;
54  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
55 
56  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
57  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
58 
59  RT R_inv = RT(1.0)/a_R;
60  RT R_inv_cubed = R_inv*R_inv*R_inv;
61 
62  RT alpha = (a_S-RT(1.0)) / (F_k + F_d);
63  RT retval = alpha*R_inv;
64 
65  RT beta = -(a_a/a_T) / (F_k + F_d);
66  retval += beta*R_inv*R_inv;
67 
68  RT gamma = a_b * a_N_s / (F_k + F_d);
69  retval += gamma*R_inv_cubed*R_inv;
70 
71  return a_f_v * retval;
72  }
73 
74  };
75 #endif
76 
77  /*! \brief Phase change equation (in terms of R^2) */
78  template <typename RT /*!< real-type */ >
79  struct dRsqdt
80  {
81  RT L; /*!< latent heat of vaporization (condensate) */
82  RT K; /*!< thermal conductivity (condensate) */
83  RT Rv; /*!< gas constant of air with vapour */
84  RT rho_l; /*!< density of condensate */
85 
86  /*! \brief Right-hand-side of the phase change ODE */
87  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
88  RT rhs_func ( const RT a_R_sq, /*!< radius squared */
89  const RT a_S, /*!< saturation ratio */
90  const RT a_T, /*!< temperature */
91  const RT a_e_s, /*!< saturation pressure */
92  const RT a_D, /*!< mol. diffusion coeff */
93  const RT a_a, /*!< curvature coeff */
94  const RT a_b, /*!< solute coeff */
95  const RT a_N_s, /*!< total solute moles */
96  const RT a_f_v /*!< ventilation factor */ ) const noexcept
97  {
98  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
99  RT Kn = lambda_v/std::sqrt(a_R_sq);
100  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
101 
102  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
103  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
104 
105  RT R_inv = RT(1.0)/std::sqrt(a_R_sq);
106  RT R_inv_cubed = R_inv*R_inv*R_inv;
107 
108  RT alpha = RT(2.0) * (a_S-RT(1.0)) / (F_k + F_d);
109  RT retval = alpha;
110 
111  RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d);
112  retval += beta*R_inv;
113 
114  RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d);
115  retval += gamma*R_inv_cubed;
116 
117  return a_f_v * retval;
118  }
119 
120  /*! \brief Jacobian of right-hand-side of the phase change ODE */
121  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
122  RT rhs_jac ( const RT a_R_sq, /*!< radius squared */
123  const RT a_T, /*!< temperature */
124  const RT a_e_s, /*!< saturation pressure */
125  const RT a_D, /*!< mol. diffusion coeff */
126  const RT a_a, /*!< curvature coeff */
127  const RT a_b, /*!< solute coeff */
128  const RT a_N_s, /*!< total solute moles */
129  const RT a_f_v /*!< ventilation factor */ ) const noexcept
130  {
131  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
132  RT Kn = lambda_v/std::sqrt(a_R_sq);
133  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
134 
135  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
136  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
137 
138  RT R_inv = RT(1.0)/std::sqrt(a_R_sq);
139  RT R_inv_3 = R_inv*R_inv*R_inv;
140  RT R_inv_5 = R_inv_3*R_inv*R_inv;
141 
142  RT retval = RT(0.0);
143 
144  RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d);
145  retval -= myhalf * beta*R_inv_3;
146 
147  RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d);
148  retval -= myhalf * three*gamma*R_inv_5;
149 
150  return a_f_v * retval;
151  }
152 
153  };
154 
155  /*! \brief Scalar Newton solver for phase change equation
156  *
157  * Solves the following nonlinear equation:
158  * mu * u - F(u) - R = 0,
159  * where:
160  * u: solution variable
161  * mu: constant
162  * R: right-hand-side (constant)
163  * F(u): function
164  */
165  template<typename NE /*!< Nonlinear equation */, typename RT /*!< real-type */>
166  struct NewtonSolver
167  {
168  const NE m_ne; /*!< nonlinear equation */
169 
170  RT m_rtol; /*!< relative tolerance */
171  RT m_atol; /*!< absolute tolerance */
172  RT m_stol; /*!< step size tolerance */
173  int m_maxits; /*!< max number of iterations */
174 
175  /*! \brief solve the nonlinear equation */
176  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
177  void operator() ( RT& a_u, /*!< solution variable */
178  RT& a_r, /*!< right-hand-side */
179  const RT& a_mu, /*!< mu */
180  const RT& a_S, /*!< saturation ratio */
181  const RT& a_T, /*!< temperature */
182  const RT& a_e_s, /*!< saturation pressure */
183  const RT& a_D, /*!< mol. diff. coeff */
184  const RT& a_a, /*!< curvature coeff */
185  const RT& a_b, /*!< solute coeff */
186  const RT& a_N_s, /*!< total solute moles */
187  const RT& a_f_v, /*!< ventilation factor */
188  RT& a_res_norm_a, /*!< absolute norm at exit */
189  RT& a_res_norm_r /*!< relative norm at exit */,
190  bool& a_converged /*!< convergence status at exit */
191  ) const
192  {
193  a_converged = false;
194  RT res_norm0 = RT(0.0);
195 
196  for (int k = 0; k < m_maxits; k++) {
197  RT residual = a_mu * a_u
198  - ( a_r
199  + m_ne.rhs_func( a_u, a_S, a_T, a_e_s, a_D, a_a, a_b, a_N_s, a_f_v ) );
200  a_res_norm_a = std::sqrt(residual*residual);
201 
202  if (k == 0) {
203  if (a_res_norm_a > 0) {
204  res_norm0 = a_res_norm_a;
205  } else {
206  res_norm0 = RT(1.0);
207  }
208  }
209  a_res_norm_r = a_res_norm_a / res_norm0;
210 
211  if (a_res_norm_a <= m_atol) {
212  a_converged = true;
213  break;
214  }
215  if (a_res_norm_r <= m_rtol) {
216  a_converged = true;
217  break;
218  }
219  if (!std::isfinite(a_res_norm_a)) {
220  a_converged = false;
221  break;
222  }
223 
224  RT slope = a_mu - m_ne.rhs_jac( a_u, a_T, a_e_s, a_D, a_a, a_b, a_N_s, a_f_v );
225  RT du = RT(0.0);
226  du = - residual / slope;
227 
228  RT du_norm = std::sqrt(du*du);
229  RT u_norm = std::sqrt(a_u*a_u);
230  if (du_norm/u_norm <= m_stol) {
231  a_converged = true;
232  break;
233  }
234 
235  a_u += du;
236  if (a_u <= RT(0.0)) {
237  a_converged = false;
238  break;
239  }
240  }
241  }
242  };
243 
244  /*! \brief Implicit and explicit time integrators for the phase change equation */
245  template< typename ODE /*!< ODE */,
246  typename NewtonSolver /*!< Newton solver */,
247  /**
248  * @brief Implicit and explicit time integrators for the phase change equation.
249  *
250  * @tparam ODE ODE type.
251  * @tparam NewtonSolver Newton solver type.
252  * @tparam RT Real-type.
253  */
254  typename RT /*!< real-type */ >
255  struct TI
256  {
257  const ODE m_ode; /*!< ODE */
258  const NewtonSolver m_newton; /*!< Newton solver */
259 
260  RT m_t_final; /*!< final time */
261  RT m_max_steps; /*!< max number of timesteps */
262  RT m_S; /*!< saturation ratio */
263  RT m_T; /*!< temperature */
264  RT m_e_s; /*!< saturation pressure */
265  RT m_D; /*!< mol. diff. coeff */
266  RT m_a; /*!< coefficient of curvature */
267  RT m_b; /*!< coefficient of solute */
268  RT m_N_s; /*!< total solute moles */
269 
270  RT m_cfl; /*!< CFL */
271  RT m_atol; /*!< absolute tolerance (for adaptive dt) */
272  RT m_rtol; /*!< absolute tolerance (for adaptive dt) */
273  RT m_stol; /*!< solution update tolerance for exit due to steady state */
274 
275  bool m_adapt_dt; /*!< use error-based adaptive dt? */
276  bool m_verbose; /*!< verbosity */
277 
278  bool m_ventilation; /*!< include ventilation factor? */
279  RT m_vent_alpha1; /*!< ventilation factor fit coefficient alpha_1 */
280  RT m_vent_beta1; /*!< ventilation factor fit exponent beta_1 */
281  RT m_vent_alpha2; /*!< ventilation factor fit coefficient alpha_2 */
282  RT m_vent_beta2; /*!< ventilation factor fit exponent beta_2 */
283  RT m_vent_fcap; /*!< maximum ventilation factor */
284 
285  /*! \brief Ventilation factor at the given radius-squared, frozen over a
286  * step to keep the implicit Jacobian exact (unity if disabled) */
287  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
288  RT ventFactor (const RT& a_u) const
289  {
290  if (!m_ventilation) { return RT(1.0); }
291  return ventilationFactor( std::sqrt(a_u),
292  m_vent_alpha1, m_vent_beta1,
293  m_vent_alpha2, m_vent_beta2,
294  m_vent_fcap );
295  }
296 
297  /*! \brief Compute timestep from stiffness estimate */
298  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
299  RT computeTimestep (const RT& a_u) const
300  {
301  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
302  return m_cfl / std::sqrt(tau*tau);
303  }
304 
305  /*! \brief Compute stiffness estimate (tau) */
306  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
307  RT computeTau (const RT& a_u) const
308  {
309  return m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
310  }
311 
312  /*! \brief Limit timestep to not exceed final time */
313  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
314  RT limitTimestep (RT dt, RT cur_time) const
315  {
316  if ((cur_time + dt) > m_t_final) {
317  dt = m_t_final - cur_time;
318  }
319  return dt;
320  }
321 
322  /*! \brief Check if timestep is too small to continue */
323  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
324  bool isTimestepTooSmall (RT dt, RT tau) const
325  {
326  return (dt < (RT(1.0e-12) * m_cfl / std::sqrt(tau*tau)))
327  && (dt < (RT(1.0e-12) * m_t_final));
328  }
329 
330  /*! \brief Evaluate ODE right-hand side */
331  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
332  RT evalRHS (const RT& a_u) const
333  {
334  return m_ode.rhs_func(a_u, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
335  }
336 
337  /*! \brief Print verbose step information */
338  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
339  void printStepInfo (RT cur_time, RT dt, RT tau, RT radius, RT snorm) const
340  {
341  if (m_verbose) {
342  AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
343  cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm);
344  }
345  }
346 
347  /*! \brief Print verbose step info with Newton solver details */
348  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
349  void printStepInfoNewton (RT cur_time, RT dt, RT tau, RT radius, RT snorm,
350  RT res_norm_a, RT res_norm_r, bool converged) const
351  {
352  if (m_verbose) {
353  AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
354  cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm);
355  AMREX_DEVICE_PRINTF(" norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
356  res_norm_a, res_norm_r, (converged ? "yes" : "no"));
357  }
358  }
359 
360  /*! \brief 3rd-order, 4-stage Bogacki-Shampine explicit RK method
361  * with 2nd order embedded method */
362  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
363  void rk3bs ( RT& a_u, /*!< solution */
364  bool& a_success /*!< success/failure flag */ ) const
365  {
366  RT cur_time = RT(0.0);
367  a_success = true;
368 
369  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
370  RT dt = m_cfl / std::sqrt(tau*tau);
371 
372  RT dt_new = dt;
373  RT a_u_old = a_u;
374 
375  int n_step = 0;
376  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
377 
378  if (!m_adapt_dt) {
379  tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
380  dt = m_cfl / std::sqrt(tau*tau);
381  } else {
382  dt = dt_new;
383  }
384 
385  if ((cur_time + dt) > m_t_final) {
386  dt = m_t_final - cur_time;
387  }
388  if (!std::isfinite(dt)) {
389  a_success = false;
390  break;
391  }
392 
393  RT u_new = RT(0.0);
394  bool step_success = false;
395  while (!step_success) {
396 
397  if ( (dt < (RT(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
398  && (dt < (RT(1.0e-12)*m_t_final)) ) {
399  break;
400  }
401 
402  RT u1 = a_u;
403  if (u1 <= 0) {
404  dt *= myhalf;
405  continue;
406  }
407  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
408 
409  RT u2 = a_u + myhalf*dt*f1;
410  if (u2 <= 0) {
411  dt *= myhalf;
412  continue;
413  }
414  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
415 
416  RT u3 = a_u + RT(0.75)*dt*f2;
417  if (u3 <= 0) {
418  dt *= myhalf;
419  continue;
420  }
421  RT f3 = m_ode.rhs_func(u3, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
422 
423  RT u4 = a_u + (RT(1.0)/RT(9.0))*dt * (RT(2.0)*f1 + three*f2 + RT(4.0)*f3);
424  if (u4 <= 0) {
425  dt *= myhalf;
426  continue;
427  }
428  RT f4 = m_ode.rhs_func(u4, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u));
429 
430  u_new = u4;
431 
432  if (m_adapt_dt) {
433  RT u_embed = a_u + (RT(1.0)/RT(24.0))*dt * (RT(7.0)*f1 + RT(6.0)*f2 + RT(8.0)*f3 + three*f4);
434  RT err = std::sqrt((u_new-u_embed)*(u_new-u_embed));
435  RT tol = m_atol + m_rtol * std::max(a_u, a_u_old);
436  RT E = err / tol;
437  dt_new = dt / std::cbrt(E);
438  }
439 
440  if (std::isfinite(u_new)) {
441  if (u_new > 0) {
442  step_success = true;
443  break;
444  }
445  }
446  dt *= myhalf;
447  }
448 
449  if (step_success) {
450 
451  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
452  a_u_old = a_u;
453  a_u = u_new;
454  cur_time += dt;
455 
456  if (m_verbose) {
457  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
458  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
459  }
460  if (snorm < m_stol) {
461  break;
462  }
463 
464  } else {
465 
466  a_success = false;
467  break;
468 
469  }
470 
471  n_step++;
472  }
473 
474  return;
475  }
476 
477  /*! \brief 4th-order, 4-stage explicit Runge-Kutta method */
478  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
479  void rk4 ( RT& a_u, /*!< solution */
480  bool& a_success /*!< success/failure flag */ ) const
481  {
482  RT cur_time = RT(0.0);
483  a_success = true;
484 
485  int n_step = 0;
486  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
487 
488  RT tau = computeTau(a_u);
489  RT dt = m_cfl / std::sqrt(tau*tau);
490  dt = limitTimestep(dt, cur_time);
491  if (!std::isfinite(dt)) {
492  a_success = false;
493  break;
494  }
495 
496  RT u_new = RT(0.0);
497  bool step_success = false;
498  while (!step_success) {
499 
500  if (isTimestepTooSmall(dt, tau)) { break; }
501 
502  RT u1 = a_u;
503  if (u1 <= 0) { dt *= myhalf; continue; }
504  RT f1 = evalRHS(u1);
505 
506  RT u2 = a_u + myhalf*dt*f1;
507  if (u2 <= 0) { dt *= myhalf; continue; }
508  RT f2 = evalRHS(u2);
509 
510  RT u3 = a_u + myhalf*dt*f2;
511  if (u3 <= 0) { dt *= myhalf; continue; }
512  RT f3 = evalRHS(u3);
513 
514  RT u4 = a_u + dt*f3;
515  if (u4 <= 0) { dt *= myhalf; continue; }
516  RT f4 = evalRHS(u4);
517 
518  u_new = a_u + dt*(f1 + RT(2.0)*f2 + RT(2.0)*f3 + f4)/RT(6.0);
519 
520  if (std::isfinite(u_new) && (u_new > 0)) {
521  step_success = true;
522  break;
523  }
524  dt *= myhalf;
525  }
526 
527  if (step_success) {
528  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
529  a_u = u_new;
530  cur_time += dt;
531  printStepInfo(cur_time, dt, tau, std::sqrt(a_u), snorm);
532  if (snorm < m_stol) { break; }
533  } else {
534  a_success = false;
535  break;
536  }
537 
538  n_step++;
539  }
540 
541  return;
542  }
543 
544  /*! \brief 1st order implicit backward Euler method */
545  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
546  void be ( RT& a_u, /*!< solution */
547  bool& a_success /*!< success/failure flag */ ) const
548  {
549  RT cur_time = RT(0.0);
550  a_success = true;
551 
552  int n_step = 0;
553  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
554 
555  RT tau = computeTau(a_u);
556  RT dt = m_cfl / std::sqrt(tau*tau);
557  dt = limitTimestep(dt, cur_time);
558  if (!std::isfinite(dt)) {
559  a_success = false;
560  break;
561  }
562 
563 #ifdef AMREX_USE_FLOAT
564  RT res_norm_a = std::numeric_limits<float>::max();
565  RT res_norm_r = std::numeric_limits<float>::max();
566 #else
567  RT res_norm_a = std::numeric_limits<double>::max();
568  RT res_norm_r = std::numeric_limits<double>::max();
569 #endif
570  bool converged = false;
571 
572  RT u_new = RT(0.0);
573  bool step_success = false;
574  while (!step_success) {
575 
576  if (isTimestepTooSmall(dt, tau)) { break; }
577 
578  RT mu = RT(1.0) / dt;
579  RT rhs = mu * a_u;
580  u_new = a_u;
581  m_newton(u_new, rhs, mu,
582  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
583  res_norm_a, res_norm_r, converged);
584 
585  if (converged && std::isfinite(u_new) && (u_new > 0)) {
586  step_success = true;
587  break;
588  }
589  dt *= myhalf;
590  }
591 
592  if (step_success) {
593  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
594  a_u = u_new;
595  cur_time += dt;
596  printStepInfoNewton(cur_time, dt, tau, std::sqrt(a_u), snorm,
597  res_norm_a, res_norm_r, converged);
598  if (snorm < m_stol) { break; }
599  } else {
600  a_success = false;
601  break;
602  }
603 
604  n_step++;
605  }
606 
607  return;
608  }
609 
610  /*! \brief 2nd-order implicit Crank-Nicolson method */
611  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
612  void cn ( RT& a_u, /*!< solution */
613  bool& a_success /*!< success/failure flag */ ) const
614  {
615  RT cur_time = RT(0.0);
616  a_success = true;
617 
618  int n_step = 0;
619  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
620 
621  RT tau = computeTau(a_u);
622  RT dt = m_cfl / std::sqrt(tau*tau);
623  dt = limitTimestep(dt, cur_time);
624  if (!std::isfinite(dt)) {
625  a_success = false;
626  break;
627  }
628 
629  RT res_norm_a = std::numeric_limits<RT>::max();
630  RT res_norm_r = std::numeric_limits<RT>::max();
631  bool converged = false;
632 
633  RT u_new = RT(0.0);
634  bool step_success = false;
635  while (!step_success) {
636 
637  if (isTimestepTooSmall(dt, tau)) { break; }
638 
639  RT mu = RT(1.0) / (myhalf*dt);
640  RT u1 = a_u;
641  RT f1 = evalRHS(u1);
642 
643  RT u2 = u1;
644  RT rhs = mu * (a_u + myhalf*dt*f1);
645  m_newton(u2, rhs, mu,
646  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
647  res_norm_a, res_norm_r, converged);
648  if (u2 <= 0) { dt *= myhalf; continue; }
649  RT f2 = evalRHS(u2);
650 
651  u_new = a_u + myhalf * dt * (f1 + f2);
652 
653  if (converged && std::isfinite(u_new) && (u_new > 0)) {
654  step_success = true;
655  break;
656  }
657  dt *= myhalf;
658  }
659 
660  if (step_success) {
661  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
662  a_u = u_new;
663  cur_time += dt;
664  printStepInfoNewton(cur_time, dt, tau, std::sqrt(a_u), snorm,
665  res_norm_a, res_norm_r, converged);
666  if (snorm < m_stol) { break; }
667  } else {
668  a_success = false;
669  break;
670  }
671 
672  n_step++;
673  }
674 
675  return;
676  }
677 
678  /*! \brief 2nd-order, 2-stage diagonally-implicit RK method */
679  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
680  void dirk212 ( RT& a_u, /*!< solution */
681  bool& a_success /*!< success/failure flag */ ) const
682  {
683  RT cur_time = RT(0.0);
684  a_success = true;
685 
686  int n_step = 0;
687  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
688 
689  RT tau = computeTau(a_u);
690  RT dt = m_cfl / std::sqrt(tau*tau);
691  dt = limitTimestep(dt, cur_time);
692  if (!std::isfinite(dt)) {
693  a_success = false;
694  break;
695  }
696 
697  bool converged = false;
698  RT res_norm_a = std::numeric_limits<RT>::max();
699  RT res_norm_r = std::numeric_limits<RT>::max();
700 
701  RT u_new = RT(0.0);
702  bool step_success = false;
703  while (!step_success) {
704 
705  if (isTimestepTooSmall(dt, tau)) {
706  break;
707  }
708  RT mu = RT(1.0) / dt;
709 
710  converged = true;
711  res_norm_a = RT(0.0);
712  res_norm_r = RT(0.0);
713 
714  RT u1 = a_u;
715  {
716  RT rhs = mu * a_u;
717  RT res_norm_a_i = std::numeric_limits<RT>::max();
718  RT res_norm_r_i = std::numeric_limits<RT>::max();
719  bool converged_i = false;
720  m_newton( u1, rhs, mu,
721  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
722  res_norm_a_i, res_norm_r_i, converged_i );
723  converged = converged && converged_i;
724  res_norm_a = std::max(res_norm_a, res_norm_a_i);
725  res_norm_r = std::max(res_norm_r, res_norm_r_i);
726  }
727  if (u1 <= 0) {
728  dt *= myhalf;
729  continue;
730  }
731  RT f1 = evalRHS(u1);
732 
733  RT u2 = u1;
734  {
735  RT rhs = mu * (a_u - dt*f1);
736  RT res_norm_a_i = std::numeric_limits<RT>::max();
737  RT res_norm_r_i = std::numeric_limits<RT>::max();
738  bool converged_i = false;
739  m_newton( u2, rhs, mu,
740  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s, ventFactor(a_u),
741  res_norm_a_i, res_norm_r_i, converged_i );
742  converged = converged && converged_i;
743  res_norm_a = std::max(res_norm_a, res_norm_a_i);
744  res_norm_r = std::max(res_norm_r, res_norm_r_i);
745  }
746  if (u2 <= 0) {
747  dt *= myhalf;
748  continue;
749  }
750  RT f2 = evalRHS(u2);
751 
752  u_new = a_u + myhalf * dt * (f1 + f2);
753 
754  if (converged) {
755  if (std::isfinite(u_new)) {
756  if (u_new > 0) {
757  step_success = true;
758  break;
759  }
760  }
761  }
762  dt *= myhalf;
763  }
764 
765  if (step_success) {
766 
767  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
768  a_u = u_new;
769  cur_time += dt;
770 
771  printStepInfoNewton(cur_time, dt, tau, std::sqrt(a_u), snorm,
772  res_norm_a, res_norm_r, converged);
773  if (snorm < m_stol) {
774  break;
775  }
776 
777  } else {
778 
779  a_success = false;
780  break;
781 
782  }
783 
784  n_step++;
785  }
786 
787  return;
788  }
789 
790  };
791 }
792 
793 /*! \brief Namespace with classes and functions for deposition/sublimation
794  * (solid <--> vapour) */
795 namespace SDMassChangeUtils_SV
796 {
797  /*! \brief Phase change equation (in terms of mass) */
798  template <typename RT /*!< real-type */ >
799  struct dMdt
800  {
801  RT L; /*!< latent heat of vaporization (condensate) */
802  RT K; /*!< thermal conductivity */
803  RT Rv; /*!< gas constant of air with vapour */
804  RT rho_ice; /*!< True density of ice */
805 
806  /*! \brief Compute growth ratio */
807  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
808  RT growthRatio ( const RT a_T /*!< temperature */ ) const noexcept
809  {
810  // Table of the inherent growth ratio (Chen and Lamb 1994)
811  // T(n) = (1.375-n)/4. Here, T is the temperature in [degreeC], and n is the index.
812  static constexpr amrex::Real tb_habit[121] = {
813  amrex::Real(1.000000e+00), amrex::Real(9.799412e-01), amrex::Real(9.600636e-01), amrex::Real(9.399397e-01), amrex::Real(9.200258e-01),
814  amrex::Real(9.001191e-01), amrex::Real(8.800351e-01), amrex::Real(8.570378e-01), amrex::Real(8.340652e-01), amrex::Real(8.109611e-01),
815  amrex::Real(7.830689e-01), amrex::Real(7.550922e-01), amrex::Real(7.030723e-01), amrex::Real(5.370318e-01), amrex::Real(4.677351e-01),
816  amrex::Real(5.248075e-01), amrex::Real(6.309573e-01), amrex::Real(8.128305e-01), amrex::Real(1.096478e+00), amrex::Real(1.479108e+00),
817  amrex::Real(1.905461e+00), amrex::Real(2.089296e+00), amrex::Real(2.290868e+00), amrex::Real(2.398833e+00), amrex::Real(2.454709e+00),
818  amrex::Real(2.426610e+00), amrex::Real(2.371374e+00), amrex::Real(2.290868e+00), amrex::Real(2.137962e+00), amrex::Real(1.995262e+00),
819  amrex::Real(1.862087e+00), amrex::Real(1.737801e+00), amrex::Real(1.621810e+00), amrex::Real(1.513561e+00), amrex::Real(1.396368e+00),
820  amrex::Real(1.288250e+00), amrex::Real(1.188502e+00), amrex::Real(1.096478e+00), amrex::Real(1.000000e+00), amrex::Real(9.225714e-01),
821  amrex::Real(8.511380e-01), amrex::Real(7.852356e-01), amrex::Real(7.244360e-01), amrex::Real(6.683439e-01), amrex::Real(6.165950e-01),
822  amrex::Real(5.754399e-01), amrex::Real(5.370318e-01), amrex::Real(5.011872e-01), amrex::Real(4.677351e-01), amrex::Real(4.365158e-01),
823  amrex::Real(4.073803e-01), amrex::Real(3.801894e-01), amrex::Real(3.548134e-01), amrex::Real(3.311311e-01), amrex::Real(3.162278e-01),
824  amrex::Real(3.019952e-01), amrex::Real(2.917427e-01), amrex::Real(2.851018e-01), amrex::Real(2.818383e-01), amrex::Real(2.786121e-01),
825  amrex::Real(2.754229e-01), amrex::Real(2.786121e-01), amrex::Real(2.818383e-01), amrex::Real(2.851018e-01), amrex::Real(2.917427e-01),
826  amrex::Real(2.985383e-01), amrex::Real(3.090295e-01), amrex::Real(3.198895e-01), amrex::Real(3.311311e-01), amrex::Real(3.467369e-01),
827  amrex::Real(3.672823e-01), amrex::Real(3.935501e-01), amrex::Real(4.265795e-01), amrex::Real(4.570882e-01), amrex::Real(4.897788e-01),
828  amrex::Real(5.248075e-01), amrex::Real(5.623413e-01), amrex::Real(6.095369e-01), amrex::Real(6.606934e-01), amrex::Real(7.161434e-01),
829  amrex::Real(7.852356e-01), amrex::Real(8.609938e-01), amrex::Real(9.549926e-01), amrex::Real(1.047129e+00), amrex::Real(1.148154e+00),
830  amrex::Real(1.258925e+00), amrex::Real(1.380384e+00), amrex::Real(1.496236e+00), amrex::Real(1.603245e+00), amrex::Real(1.698244e+00),
831  amrex::Real(1.778279e+00), amrex::Real(1.840772e+00), amrex::Real(1.883649e+00), amrex::Real(1.905461e+00), amrex::Real(1.905461e+00),
832  amrex::Real(1.883649e+00), amrex::Real(1.862087e+00), amrex::Real(1.840772e+00), amrex::Real(1.798871e+00), amrex::Real(1.737801e+00),
833  amrex::Real(1.698244e+00), amrex::Real(1.640590e+00), amrex::Real(1.584893e+00), amrex::Real(1.549173e+00), amrex::Real(1.513910e+00),
834  amrex::Real(1.476046e+00), amrex::Real(1.452112e+00), amrex::Real(1.428894e+00), amrex::Real(1.412863e+00), amrex::Real(1.393157e+00),
835  amrex::Real(1.376892e+00), amrex::Real(1.361131e+00), amrex::Real(1.348963e+00), amrex::Real(1.336903e+00), amrex::Real(1.327089e+00),
836  amrex::Real(1.317953e+00), amrex::Real(1.308881e+00), amrex::Real(1.302867e+00), amrex::Real(1.293898e+00), amrex::Real(1.287953e+00),
837  amrex::Real(1.279087e+00)
838  };
839  auto T_degC = a_T - RT(273.15); // [K] => [degC]
840  int n = std::min(std::max(static_cast<int>(std::round(RT(1.375)-RT(4.0)*T_degC)),1),121);
841  return static_cast<RT>(tb_habit[n-1]);
842 
843  }
844 
845  /*! \brief Compute F_k + F_d */
846  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
847  RT Fk_plus_Fd ( const RT a_T, /*!< temperature */
848  const RT a_e_s, /*!< saturation pressure */
849  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
850  {
851  auto F_k = (L/(Rv*a_T) - RT(1.0)) * (L / (K*a_T));
852  auto F_d = (Rv*a_T) / (a_D*a_e_s);
853  return (F_k + F_d);
854  }
855 
856  /*! \brief Compute dynamic viscosity */
857  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
858  RT viscCoeff ( const RT a_T /*!< temperature */ ) const noexcept
859  {
860  auto T_degC = a_T - RT(273.15); // [K] => [degC]
861  RT visc_coeff = RT(0.0);
862  if( T_degC >= RT(0.0) ) {
863  visc_coeff = ( RT(1.7180) + RT(4.9E-3)*T_degC ) * RT(1.E-5);
864  } else {
865  visc_coeff = ( RT(1.7180) + RT(4.9E-3)*T_degC -RT(1.2E-5)*T_degC*T_degC ) * RT(1.E-5);
866  }
867  return visc_coeff;
868  }
869 
870  /*! \brief Compute capacitance */
871  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
872  RT capacitance( const RT a_a, /*!< equatorial radius */
873  const RT a_c /*!< polar radius */ ) const noexcept
874  {
875  auto phi = a_c / a_a;
876  RT capaci = RT(0.0);
877  if(phi >= RT(1.001)) {
878  auto ecent = std::sqrt(RT(1.0) - RT(1.0)/(phi*phi));
879  capaci = ecent*a_c / std::log((RT(1.0)+ecent)*phi);
880  } else if (phi <= RT(0.999)) {
881  auto ecent = std::sqrt(RT(1.0) - phi*phi);
882  capaci = ecent*a_a / std::asin(ecent);
883  } else {
884  capaci = (RT(2.0)*a_a + a_c)/three;
885  }
886  return capaci;
887  }
888 
889  /*! \brief Compute ventilation effect */
890  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
891  RT ventX ( const RT a_a, /*!< equatorial radius */
892  const RT a_c, /*!< polar radius */
893  const RT a_vterm,/*!< terminal velocity */
894  const RT a_T, /*!< temperature */
895  const RT a_rhom, /*!< moist density */
896  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
897  {
898  auto visc_coeff = viscCoeff(a_T);
899  auto maxdim = RT(2.0) * std::max(a_a, a_c);
900  auto n_re = a_rhom * maxdim * a_vterm / visc_coeff;
901  auto n_sc = visc_coeff / (a_rhom*a_D);
902  auto X_vent = std::cbrt(n_sc) * std::sqrt(n_re);
903  return X_vent;
904  }
905 
906  /*! \brief Radius change */
907  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
908  RT growthRatioStar ( const RT a_dmass, /*!< delta mass */
909  const RT a_a, /*!< equatorial radius */
910  const RT a_c, /*!< polar radius */
911  const RT a_vterm, /*!< terminal velocity */
912  const RT a_rhom, /*!< moist air density */
913  const RT a_T, /*!< temperature */
914  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
915  {
916  auto capaci = capacitance(a_a, a_c);
917  auto X_vent = ventX(a_a, a_c, a_vterm, a_T, a_rhom, a_D);
918  RT fratio = RT(0.0);
919  if(X_vent <= RT(1.0)) {
920  fratio = (RT(1.0) + RT(0.14)*X_vent*X_vent*std::sqrt(a_c/capaci))
921  / (RT(1.0) + RT(0.14)*X_vent*X_vent*std::sqrt(a_a/capaci));
922  } else {
923  fratio = (RT(0.86) + RT(0.28)*X_vent*std::sqrt(a_c/capaci))
924  / (RT(0.86) + RT(0.28)*X_vent*std::sqrt(a_a/capaci));
925  }
926 
927  auto gr = growthRatio(a_T);
928  if ((RT(2.0)*std::max(a_a,a_c) < RT(1.0e-5)) || (a_dmass <= RT(0.0))) { gr = RT(1.0); }
929  auto gr_star = gr * fratio;
930  if((a_dmass > RT(0.0)) && (a_c/a_a > RT(40.0))) { gr_star = std::min(gr_star,RT(1.0)); }
931 
932  return gr_star;
933  }
934 
935  /*! \brief Volume change */
936  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
937  RT dVolume ( const RT a_dmass, /*!< delta mass */
938  const RT a_a, /*!< equatorial radius */
939  const RT a_c, /*!< polar radius */
940  const RT a_rho, /*!< ice density of particle */
941  const RT a_S, /*!< saturation ratio */
942  const RT a_T, /*!< temperature */
943  const RT a_e_s, /*!< saturation pressure */
944  const RT a_es_ratio_wi, /*!< ratio of saturation pressures of water to ice */
945  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
946  {
947  RT rho_i_ds = a_rho;
948  auto gr = growthRatio(a_T);
949  if (a_dmass > RT(0.0)) {
950  if ((gr < RT(1.0)) && (a_a < RT(1.0e-4))) {
951  rho_i_ds = rho_ice;
952  } else {
953  if (RT(2.0)*std::max(a_a,a_c) < RT(1.0e-5)) { gr = RT(1.0); }
954  auto ex_vap_dens = RT(1000.0) * (std::min(a_S,a_es_ratio_wi) - RT(1.0)) / (Fk_plus_Fd(a_T,a_e_s,a_D) * a_D);
955  rho_i_ds = rho_ice * std::exp(-three*std::max(ex_vap_dens-RT(0.05),RT(0.0))/gr);
956  }
957  }
958 
959  auto dvol = a_dmass / rho_i_ds;
960  return dvol;
961  }
962 
963  /*! \brief Radius change */
964  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
965  RT dLogRadius ( const RT a_grstar, /*!< growth ratio */
966  const RT a_dlogvol /*!< delta log volume */ ) const noexcept
967  {
968  return a_dlogvol / (a_grstar+RT(2.0));
969  }
970 
971  /*! \brief Right-hand-side of the phase change ODE */
972  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
973  RT rhs_func ( const RT a_a, /*!< equatorial radius */
974  const RT a_c, /*!< polar radius */
975  const RT a_vterm,/*!< terminal velocity */
976  const RT a_S, /*!< saturation ratio */
977  const RT a_rhom, /*!< moist density */
978  const RT a_T, /*!< temperature */
979  const RT a_e_s, /*!< saturation pressure */
980  const RT a_D /*!< mol. diffusion coeff */ ) const noexcept
981  {
982  auto capaci = capacitance(a_a, a_c);
983  auto X_vent = ventX(a_a, a_c, a_vterm, a_T, a_rhom, a_D);
984 
985  RT f_vent = RT(0.0);
986  if (X_vent <= RT(1.0)) { f_vent = RT(1.0) + RT(0.14)*X_vent*X_vent; }
987  else { f_vent = RT(0.86) + RT(0.28)*X_vent; }
988 
989  auto retval = RT(4.0)*PI * capaci * (a_S-RT(1.0)) * f_vent / Fk_plus_Fd(a_T, a_e_s, a_D);
990  return retval;
991  }
992 
993  /*! \brief Ice melt rate dm_ice/dt (Pruppacher-Klett; Seifert-Beheng 2006 Eq. 72).
994  * Negative when melting. Heat conducted from warm air plus latent heat of the
995  * vapour flux to the 0 degC wet surface, divided by the latent heat of fusion.
996  * Reuses capacitance and ventilation. For this functor built with the water
997  * material, L is the latent heat of vaporization, K the air thermal
998  * conductivity, and Rv the vapour gas constant. */
999  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1000  RT meltRate ( const RT a_a, /*!< equatorial radius */
1001  const RT a_c, /*!< polar radius */
1002  const RT a_vterm, /*!< terminal velocity */
1003  const RT a_T, /*!< temperature */
1004  const RT a_e_inf, /*!< ambient vapour pressure */
1005  const RT a_e_sat_wT0, /*!< sat. pressure over water at melting point */
1006  const RT a_rhom, /*!< moist density */
1007  const RT a_D, /*!< vapour diffusion coeff */
1008  const RT a_L_f /*!< latent heat of fusion */ ) const noexcept
1009  {
1010  auto capaci = capacitance(a_a, a_c);
1011  auto X_vent = ventX(a_a, a_c, a_vterm, a_T, a_rhom, a_D);
1012 
1013  RT f_vent = RT(0.0);
1014  if (X_vent <= RT(1.0)) { f_vent = RT(1.0) + RT(0.14)*X_vent*X_vent; }
1015  else { f_vent = RT(0.86) + RT(0.28)*X_vent; }
1016 
1017  RT T0 = RT(tmelt);
1018  RT heat = K * f_vent * (a_T - T0)
1019  + (L*a_D/Rv) * (a_e_inf/a_T - a_e_sat_wT0/T0) * f_vent;
1020  return -(RT(4.0)*PI*capaci/a_L_f) * heat;
1021  }
1022 
1023  };
1024 
1025  /*! \brief Implicit and explicit time integrators for the phase change equation */
1026  template< typename ODE /*!< ODE */,
1027  typename RT /*!< real-type */ >
1028  struct TI
1029  {
1030  const ODE m_ode; /*!< ODE */
1031 
1032  RT m_t_final; /*!< final time */
1033  mutable RT m_dt;/*!< dt */
1034  RT m_max_steps; /*!< max number of timesteps */
1035 
1036  RT m_a; /*!< equatorial radius */
1037  RT m_c; /*!< polar radius */
1038  RT m_vt; /*!< terminal velocity */
1039  RT m_S; /*!< saturation ratio */
1040  RT m_rhom; /*!< moist density */
1041  RT m_T; /*!< temperature */
1042  RT m_e_s; /*!< saturation pressure */
1043  RT m_D; /*!< mol. diff. coeff */
1044 
1045  bool m_verbose; /*!< verbosity */
1046 
1047  /*! \brief 1st-order, forward Euler method */
1048  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1049  void fe ( RT& a_u, /*!< solution */
1050  bool& a_success /*!< success/failure flag */ ) const
1051  {
1052  RT cur_time = RT(0.0);
1053  a_success = true;
1054 
1055  int n_step = 0;
1056  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
1057 
1058  if ((cur_time + m_dt) > m_t_final) {
1059  m_dt = m_t_final - cur_time;
1060  }
1061 
1062  RT rhs = m_ode.rhs_func(m_a, m_c, m_vt, m_S, m_rhom, m_T, m_e_s, m_D);
1063  RT u_new = a_u + m_dt*rhs;
1064  if (!std::isfinite(u_new)) {
1065  a_success = false;
1066  break;
1067  }
1068 
1069  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
1070  a_u = u_new;
1071  cur_time += m_dt;
1072 
1073  if (m_verbose) {
1074  AMREX_DEVICE_PRINTF( "Time %1.2e, m_dt = %1.2e, mass = %1.4e, snorm = %1.1e\n",
1075  cur_time, m_dt, a_u, snorm);
1076  }
1077 
1078  n_step++;
1079  }
1080 
1081  return;
1082  }
1083 
1084  };
1085 }
1086 
1087 #endif
1088 #endif
ERF_Constants.H
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
tmelt
constexpr amrex::Real tmelt
Definition: ERF_Constants.H:130
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
gamma
amrex::Real gamma
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:9
beta
amrex::Real beta
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:10
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
MORRInd::mu
@ mu
Definition: ERF_AdvanceMorrison.cpp:92
mp_wdm6::alpha
real(kind=kind_phys), parameter, private alpha
Definition: ERF_module_mp_wdm6.F90:62