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 "ERF_Constants.H"
6 
7 #ifdef ERF_USE_PARTICLES
8 
9 /*! \brief Namespace with classes and functions for condensation/evaporation
10  * (liquid <--> vapour) */
11 namespace SDMassChangeUtils_LV
12 {
13 #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS
14  /*! \brief Phase change equation (in terms of R) */
15  template <typename RT /*!< real-type */ >
16  struct dRdt
17  {
18  RT L; /*!< latent heat of vaporization (condensate) */
19  RT K; /*!< thermal conductivity (condensate) */
20  RT Rv; /*!< gas constant of air with vapour */
21  RT rho_l; /*!< density of condensate */
22 
23  /*! \brief Right-hand-side of the phase change ODE */
24  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
25  RT operator () ( const RT a_R, /*!< radius */
26  const RT a_S, /*!< saturation ratio */
27  const RT a_T, /*!< temperature */
28  const RT a_e_s, /*!< saturation pressure */
29  const RT a_D, /*!< mol. diffusion coeff */
30  const RT a_a, /*!< curvature coeff */
31  const RT a_b, /*!< solute coeff */
32  const RT a_N_s /*!< total solute moles */ ) const noexcept
33  {
34  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
35  RT Kn = lambda_v/a_R;
36  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
37 
38  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
39  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
40 
41  RT R_inv = RT(1.0)/a_R;
42  RT R_inv_cubed = R_inv*R_inv*R_inv;
43 
44  RT alpha = (a_S-RT(1.0)) / (F_k + F_d);
45  RT retval = alpha*R_inv;
46 
47  RT beta = -(a_a/a_T) / (F_k + F_d);
48  retval += beta*R_inv*R_inv;
49 
50  RT gamma = a_b * a_N_s / (F_k + F_d);
51  retval += gamma*R_inv_cubed*R_inv;
52 
53  return retval;
54  }
55 
56  };
57 #endif
58 
59  /*! \brief Phase change equation (in terms of R^2) */
60  template <typename RT /*!< real-type */ >
61  struct dRsqdt
62  {
63  RT L; /*!< latent heat of vaporization (condensate) */
64  RT K; /*!< thermal conductivity (condensate) */
65  RT Rv; /*!< gas constant of air with vapour */
66  RT rho_l; /*!< density of condensate */
67 
68  /*! \brief Right-hand-side of the phase change ODE */
69  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
70  RT rhs_func ( const RT a_R_sq, /*!< radius squared */
71  const RT a_S, /*!< saturation ratio */
72  const RT a_T, /*!< temperature */
73  const RT a_e_s, /*!< saturation pressure */
74  const RT a_D, /*!< mol. diffusion coeff */
75  const RT a_a, /*!< curvature coeff */
76  const RT a_b, /*!< solute coeff */
77  const RT a_N_s /*!< total solute moles */ ) const noexcept
78  {
79  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
80  RT Kn = lambda_v/std::sqrt(a_R_sq);
81  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
82 
83  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
84  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
85 
86  RT R_inv = RT(1.0)/std::sqrt(a_R_sq);
87  RT R_inv_cubed = R_inv*R_inv*R_inv;
88 
89  RT alpha = RT(2.0) * (a_S-RT(1.0)) / (F_k + F_d);
90  RT retval = alpha;
91 
92  RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d);
93  retval += beta*R_inv;
94 
95  RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d);
96  retval += gamma*R_inv_cubed;
97 
98  return retval;
99  }
100 
101  /*! \brief Jacobian of right-hand-side of the phase change ODE */
102  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
103  RT rhs_jac ( const RT a_R_sq, /*!< radius squared */
104  const RT a_T, /*!< temperature */
105  const RT a_e_s, /*!< saturation pressure */
106  const RT a_D, /*!< mol. diffusion coeff */
107  const RT a_a, /*!< curvature coeff */
108  const RT a_b, /*!< solute coeff */
109  const RT a_N_s /*!< total solute moles */ ) const noexcept
110  {
111  RT lambda_v = 2*a_D/std::sqrt(8*a_T*Rv/PI);
112  RT Kn = lambda_v/std::sqrt(a_R_sq);
113  RT dcf = (1+Kn) / (1+2*Kn*(1+Kn));
114 
115  RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T));
116  RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s);
117 
118  RT R_inv = RT(1.0)/std::sqrt(a_R_sq);
119  RT R_inv_3 = R_inv*R_inv*R_inv;
120  RT R_inv_5 = R_inv_3*R_inv*R_inv;
121 
122  RT retval = RT(0.0);
123 
124  RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d);
125  retval -= myhalf * beta*R_inv_3;
126 
127  RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d);
128  retval -= myhalf * three*gamma*R_inv_5;
129 
130  return retval;
131  }
132 
133  };
134 
135  /*! \brief Scalar Newton solver for phase change equation
136  *
137  * Solves the following nonlinear equation:
138  * mu * u - F(u) - R = 0,
139  * where:
140  * u: solution variable
141  * mu: constant
142  * R: right-hand-side (constant)
143  * F(u): function
144  */
145  template<typename NE /*!< Nonlinear equation */, typename RT /*!< real-type */>
146  struct NewtonSolver
147  {
148  const NE m_ne; /*!< nonlinear equation */
149 
150  RT m_rtol; /*!< relative tolerance */
151  RT m_atol; /*!< absolute tolerance */
152  RT m_stol; /*!< step size tolerance */
153  int m_maxits; /*!< max number of iterations */
154 
155  /*! \brief solve the nonlinear equation */
156  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
157  void operator() ( RT& a_u, /*!< solution variable */
158  RT& a_r, /*!< right-hand-side */
159  const RT& a_mu, /*!< mu */
160  const RT& a_S, /*!< saturation ratio */
161  const RT& a_T, /*!< temperature */
162  const RT& a_e_s, /*!< saturation pressure */
163  const RT& a_D, /*!< mol. diff. coeff */
164  const RT& a_a, /*!< curvature coeff */
165  const RT& a_b, /*!< solute coeff */
166  const RT& a_N_s, /*!< total solute moles */
167  RT& a_res_norm_a, /*!< absolute norm at exit */
168  RT& a_res_norm_r /*!< relative norm at exit */,
169  bool& a_converged /*!< convergence status at exit */
170  ) const
171  {
172  a_converged = false;
173  RT res_norm0 = RT(0.0);
174 
175  for (int k = 0; k < m_maxits; k++) {
176  RT residual = a_mu * a_u
177  - ( a_r
178  + m_ne.rhs_func( a_u, a_S, a_T, a_e_s, a_D, a_a, a_b, a_N_s ) );
179  a_res_norm_a = std::sqrt(residual*residual);
180 
181  if (k == 0) {
182  if (a_res_norm_a > 0) {
183  res_norm0 = a_res_norm_a;
184  } else {
185  res_norm0 = RT(1.0);
186  }
187  }
188  a_res_norm_r = a_res_norm_a / res_norm0;
189 
190  if (a_res_norm_a <= m_atol) {
191  a_converged = true;
192  break;
193  }
194  if (a_res_norm_r <= m_rtol) {
195  a_converged = true;
196  break;
197  }
198  if (!std::isfinite(a_res_norm_a)) {
199  a_converged = false;
200  break;
201  }
202 
203  RT slope = a_mu - m_ne.rhs_jac( a_u, a_T, a_e_s, a_D, a_a, a_b, a_N_s );
204  RT du = RT(0.0);
205  du = - residual / slope;
206 
207  RT du_norm = std::sqrt(du*du);
208  RT u_norm = std::sqrt(a_u*a_u);
209  if (du_norm/u_norm <= m_stol) {
210  a_converged = true;
211  break;
212  }
213 
214  a_u += du;
215  if (a_u <= RT(0.0)) {
216  a_converged = false;
217  break;
218  }
219  }
220  }
221  };
222 
223  /*! \brief Implicit and explicit time integrators for the phase change equation */
224  template< typename ODE /*!< ODE */,
225  typename NewtonSolver /*!< Newton solver */,
226  typename RT /*!< real-type */ >
227  struct TI
228  {
229  const ODE m_ode; /*!< ODE */
230  const NewtonSolver m_newton; /*!< Newton solver */
231 
232  RT m_t_final; /*!< final time */
233  RT m_max_steps; /*!< max number of timesteps */
234  RT m_S; /*!< saturation ratio */
235  RT m_T; /*!< temperature */
236  RT m_e_s; /*!< saturation pressure */
237  RT m_D; /*!< mol. diff. coeff */
238  RT m_a; /*!< coefficient of curvature */
239  RT m_b; /*!< coefficient of solute */
240  RT m_N_s; /*!< total solute moles */
241 
242  RT m_cfl; /*!< CFL */
243  RT m_atol; /*!< absolute tolerance (for adaptive dt) */
244  RT m_rtol; /*!< absolute tolerance (for adaptive dt) */
245  RT m_stol; /*!< solution update tolerance for exit due to steady state */
246 
247  bool m_adapt_dt; /*!< use error-based adaptive dt? */
248  bool m_verbose; /*!< verbosity */
249 
250  /*! \brief Compute timestep from stiffness estimate */
251  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
252  RT computeTimestep (const RT& a_u) const
253  {
254  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
255  return m_cfl / std::sqrt(tau*tau);
256  }
257 
258  /*! \brief Compute stiffness estimate (tau) */
259  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
260  RT computeTau (const RT& a_u) const
261  {
262  return m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
263  }
264 
265  /*! \brief Limit timestep to not exceed final time */
266  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
267  RT limitTimestep (RT dt, RT cur_time) const
268  {
269  if ((cur_time + dt) > m_t_final) {
270  dt = m_t_final - cur_time;
271  }
272  return dt;
273  }
274 
275  /*! \brief Check if timestep is too small to continue */
276  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
277  bool isTimestepTooSmall (RT dt, RT tau) const
278  {
279  return (dt < (RT(1.0e-12) * m_cfl / std::sqrt(tau*tau)))
280  && (dt < (RT(1.0e-12) * m_t_final));
281  }
282 
283  /*! \brief Evaluate ODE right-hand side */
284  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
285  RT evalRHS (const RT& a_u) const
286  {
287  return m_ode.rhs_func(a_u, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
288  }
289 
290  /*! \brief Print verbose step information */
291  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
292  void printStepInfo (RT cur_time, RT dt, RT tau, RT radius, RT snorm) const
293  {
294  if (m_verbose) {
295  AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
296  cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm);
297  }
298  }
299 
300  /*! \brief Print verbose step info with Newton solver details */
301  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
302  void printStepInfoNewton (RT cur_time, RT dt, RT tau, RT radius, RT snorm,
303  RT res_norm_a, RT res_norm_r, bool converged) const
304  {
305  if (m_verbose) {
306  AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
307  cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm);
308  AMREX_DEVICE_PRINTF(" norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
309  res_norm_a, res_norm_r, (converged ? "yes" : "no"));
310  }
311  }
312 
313  /*! \brief 3rd-order, 4-stage Bogacki-Shampine explicit RK method
314  * with 2nd order embedded method */
315  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
316  void rk3bs ( RT& a_u, /*!< solution */
317  bool& a_success /*!< success/failure flag */ ) const
318  {
319  RT cur_time = RT(0.0);
320  a_success = true;
321 
322  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
323  RT dt = m_cfl / std::sqrt(tau*tau);
324 
325  RT dt_new = dt;
326  RT a_u_old = a_u;
327 
328  int n_step = 0;
329  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
330 
331  if (!m_adapt_dt) {
332  tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
333  dt = m_cfl / std::sqrt(tau*tau);
334  } else {
335  dt = dt_new;
336  }
337 
338  if ((cur_time + dt) > m_t_final) {
339  dt = m_t_final - cur_time;
340  }
341  if (!std::isfinite(dt)) {
342  a_success = false;
343  break;
344  }
345 
346  RT u_new = zero;
347  bool step_success = false;
348  while (!step_success) {
349 
350  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
351  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
352  break;
353  }
354 
355  RT u1 = a_u;
356  if (u1 <= 0) {
357  dt *= myhalf;
358  continue;
359  }
360  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
361 
362  RT u2 = a_u + myhalf*dt*f1;
363  if (u2 <= 0) {
364  dt *= myhalf;
365  continue;
366  }
367  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
368 
369  RT u3 = a_u + amrex::Real(0.75)*dt*f2;
370  if (u3 <= 0) {
371  dt *= myhalf;
372  continue;
373  }
374  RT f3 = m_ode.rhs_func(u3, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
375 
376  RT u4 = a_u + (one/amrex::Real(9.0))*dt * (two*f1 + three*f2 + amrex::Real(4.0)*f3);
377  if (u4 <= 0) {
378  dt *= myhalf;
379  continue;
380  }
381  RT f4 = m_ode.rhs_func(u4, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
382 
383  u_new = u4;
384 
385  if (m_adapt_dt) {
386  RT u_embed = a_u + (one/amrex::Real(24.0))*dt * (amrex::Real(7.0)*f1 + amrex::Real(6.0)*f2 + amrex::Real(8.0)*f3 + three*f4);
387  RT err = std::sqrt((u_new-u_embed)*(u_new-u_embed));
388  RT tol = m_atol + m_rtol * std::max(a_u, a_u_old);
389  RT E = err / tol;
390  dt_new = dt * std::exp((one/3)*std::log(one/E));
391  }
392 
393  if (std::isfinite(u_new)) {
394  if (u_new > 0) {
395  step_success = true;
396  break;
397  }
398  }
399  dt *= myhalf;
400  }
401 
402  if (step_success) {
403 
404  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
405  a_u_old = a_u;
406  a_u = u_new;
407  cur_time += dt;
408 
409  if (m_verbose) {
410  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
411  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
412  }
413  if (snorm < m_stol) {
414  break;
415  }
416 
417  } else {
418 
419  a_success = false;
420  break;
421 
422  }
423 
424  n_step++;
425  }
426 
427  return;
428  }
429 
430  /*! \brief 4th-order, 4-stage explicit Runge-Kutta method */
431  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
432  void rk4 ( RT& a_u, /*!< solution */
433  bool& a_success /*!< success/failure flag */ ) const
434  {
435  RT cur_time = zero;
436  a_success = true;
437 
438  int n_step = 0;
439  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
440 
441  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
442  RT dt = m_cfl / std::sqrt(tau*tau);
443  if ((cur_time + dt) > m_t_final) {
444  dt = m_t_final - cur_time;
445  }
446  if (!std::isfinite(dt)) {
447  a_success = false;
448  break;
449  }
450 
451  RT u_new = zero;
452  bool step_success = false;
453  while (!step_success) {
454 
455  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
456  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
457  break;
458  }
459 
460  RT u1 = a_u;
461  if (u1 <= 0) {
462  dt *= myhalf;
463  continue;
464  }
465  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
466 
467  RT u2 = a_u + myhalf*dt*f1;
468  if (u2 <= 0) {
469  dt *= myhalf;
470  continue;
471  }
472  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
473 
474  RT u3 = a_u + myhalf*dt*f2;
475  if (u3 <= 0) {
476  dt *= myhalf;
477  continue;
478  }
479  RT f3 = m_ode.rhs_func(u3, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
480 
481  RT u4 = a_u + one*dt*f3;
482  if (u4 <= 0) {
483  dt *= myhalf;
484  continue;
485  }
486  RT f4 = m_ode.rhs_func(u4, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
487 
488  u_new = a_u + dt*(f1+two*f2+two*f3+f4)/amrex::Real(6.0);
489 
490  if (std::isfinite(u_new)) {
491  if (u_new > 0) {
492  step_success = true;
493  break;
494  }
495  }
496  dt *= myhalf;
497  }
498 
499  if (step_success) {
500 
501  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
502  a_u = u_new;
503  cur_time += dt;
504 
505  if (m_verbose) {
506  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
507  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
508  }
509  if (snorm < m_stol) {
510  break;
511  }
512 
513  } else {
514 
515  a_success = false;
516  break;
517 
518  }
519 
520  n_step++;
521  }
522 
523  return;
524  }
525 
526  /*! \brief 1st order implicit backward Euler method */
527  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
528  void be ( RT& a_u, /*!< solution */
529  bool& a_success /*!< success/failure flag */ ) const
530  {
531  RT cur_time = zero;
532  a_success = true;
533 
534  int n_step = 0;
535  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
536 
537  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
538  RT dt = m_cfl / std::sqrt(tau*tau);
539  if ((cur_time + dt) > m_t_final) {
540  dt = m_t_final - cur_time;
541  }
542  if (!std::isfinite(dt)) {
543  a_success = false;
544  break;
545  }
546 
547 #ifdef AMREX_USE_FLOAT
548  RT res_norm_a = std::numeric_limits<float>::max();
549  RT res_norm_r = std::numeric_limits<float>::max();
550 #else
551  RT res_norm_a = std::numeric_limits<double>::max();
552  RT res_norm_r = std::numeric_limits<double>::max();
553 #endif
554  bool converged = false;
555 
556  RT u_new = zero;
557  bool step_success = false;
558  while (!step_success) {
559 
560  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
561  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
562  break;
563  }
564 
565  RT mu = one / dt;
566  RT rhs = mu * a_u;
567  u_new = a_u;
568  m_newton( u_new, rhs, mu,
569  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
570  res_norm_a, res_norm_r, converged );
571 
572  if (converged) {
573  if (std::isfinite(u_new)) {
574  if (u_new > 0) {
575  step_success = true;
576  break;
577  }
578  }
579  }
580  dt *= myhalf;
581  }
582 
583  if (step_success) {
584 
585  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
586  a_u = u_new;
587  cur_time += dt;
588 
589  if (m_verbose) {
590  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
591  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
592  AMREX_DEVICE_PRINTF( " norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
593  res_norm_a, res_norm_r,
594  (converged ? "yes" : "no") );
595  }
596  if (snorm < m_stol) {
597  break;
598  }
599 
600  } else {
601 
602  a_success = false;
603  break;
604 
605  }
606 
607  n_step++;
608  }
609 
610  return;
611  }
612 
613  /*! \brief 2nd-order implicit Crank-Nicolson method */
614  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
615  void cn ( RT& a_u, /*!< solution */
616  bool& a_success /*!< success/failure flag */ ) const
617  {
618  RT cur_time = zero;
619  a_success = true;
620 
621  int n_step = 0;
622  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
623 
624  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
625  RT dt = m_cfl / std::sqrt(tau*tau);
626  if ((cur_time + dt) > m_t_final) {
627  dt = m_t_final - cur_time;
628  }
629  if (!std::isfinite(dt)) {
630  a_success = false;
631  break;
632  }
633 
634  RT res_norm_a = DBL_MAX;
635  RT res_norm_r = DBL_MAX;
636  bool converged = false;
637 
638  RT u_new = zero;
639  bool step_success = false;
640  while (!step_success) {
641 
642  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
643  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
644  break;
645  }
646  RT mu = one / (myhalf*dt);
647 
648  RT u1 = a_u;
649  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
650 
651  RT u2 = u1;
652  RT rhs = mu * (a_u + myhalf*dt*f1);
653  m_newton( u2, rhs, mu,
654  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
655  res_norm_a, res_norm_r, converged );
656  if (u2 <= 0) {
657  dt *= myhalf;
658  continue;
659  }
660  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
661 
662  u_new = a_u + myhalf * dt * (f1 + f2);
663 
664  if (converged) {
665  if (std::isfinite(u_new)) {
666  if (u_new > 0) {
667  step_success = true;
668  break;
669  }
670  }
671  }
672  dt *= myhalf;
673  }
674 
675  if (step_success) {
676 
677  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
678  a_u = u_new;
679  cur_time += dt;
680 
681  if (m_verbose) {
682  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
683  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
684  AMREX_DEVICE_PRINTF( " norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
685  res_norm_a, res_norm_r,
686  (converged ? "yes" : "no") );
687  }
688  if (snorm < m_stol) {
689  break;
690  }
691 
692  } else {
693 
694  a_success = false;
695  break;
696 
697  }
698 
699  n_step++;
700  }
701 
702  return;
703  }
704 
705  /*! \brief 2nd-order, 2-stage diagonally-implicit RK method */
706  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
707  void dirk212 ( RT& a_u, /*!< solution */
708  bool& a_success /*!< success/failure flag */ ) const
709  {
710  RT cur_time = zero;
711  a_success = true;
712 
713  int n_step = 0;
714  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
715 
716  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
717  RT dt = m_cfl / std::sqrt(tau*tau);
718  if ((cur_time + dt) > m_t_final) {
719  dt = m_t_final - cur_time;
720  }
721  if (!std::isfinite(dt)) {
722  a_success = false;
723  break;
724  }
725 
726  bool converged = false;
727  RT res_norm_a = DBL_MAX;
728  RT res_norm_r = DBL_MAX;
729 
730  RT u_new = zero;
731  bool step_success = false;
732  while (!step_success) {
733 
734  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
735  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
736  break;
737  }
738  RT mu = one / dt;
739 
740  converged = true;
741  res_norm_a = zero;
742  res_norm_r = zero;
743 
744  RT u1 = a_u;
745  {
746  RT rhs = mu * a_u;
747  RT res_norm_a_i = DBL_MAX;
748  RT res_norm_r_i = DBL_MAX;
749  bool converged_i = false;
750  m_newton( u1, rhs, mu,
751  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
752  res_norm_a_i, res_norm_r_i, converged_i );
753  converged = converged && converged_i;
754  res_norm_a = std::max(res_norm_a, res_norm_a_i);
755  res_norm_r = std::max(res_norm_r, res_norm_r_i);
756  }
757  if (u1 <= 0) {
758  dt *= myhalf;
759  continue;
760  }
761  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
762 
763  RT u2 = u1;
764  {
765  RT rhs = mu * (a_u - dt*f1);
766  RT res_norm_a_i = DBL_MAX;
767  RT res_norm_r_i = DBL_MAX;
768  bool converged_i = false;
769  m_newton( u2, rhs, mu,
770  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
771  res_norm_a_i, res_norm_r_i, converged_i );
772  converged = converged && converged_i;
773  res_norm_a = std::max(res_norm_a, res_norm_a_i);
774  res_norm_r = std::max(res_norm_r, res_norm_r_i);
775  }
776  if (u2 <= 0) {
777  dt *= myhalf;
778  continue;
779  }
780  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
781 
782  u_new = a_u + myhalf * dt * (f1 + f2);
783 
784  if (converged) {
785  if (std::isfinite(u_new)) {
786  if (u_new > 0) {
787  step_success = true;
788  break;
789  }
790  }
791  }
792  dt *= myhalf;
793  }
794 
795  if (step_success) {
796 
797  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
798  a_u = u_new;
799  cur_time += dt;
800 
801  if (m_verbose) {
802  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
803  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
804  AMREX_DEVICE_PRINTF( " norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
805  res_norm_a, res_norm_r,
806  (converged ? "yes" : "no") );
807  }
808  if (snorm < m_stol) {
809  break;
810  }
811 
812  } else {
813 
814  a_success = false;
815  break;
816 
817  }
818 
819  n_step++;
820  }
821 
822  return;
823  }
824 
825  };
826 }
827 
828 #endif
829 #endif
ERF_Constants.H
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:16
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:89
mp_wsm6::alpha
real(kind=kind_phys), parameter, public alpha
Definition: ERF_module_mp_wsm6.F90:44