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  RT res_norm_a = DBL_MAX;
548  RT res_norm_r = DBL_MAX;
549  bool converged = false;
550 
551  RT u_new = zero;
552  bool step_success = false;
553  while (!step_success) {
554 
555  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
556  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
557  break;
558  }
559 
560  RT mu = one / dt;
561  RT rhs = mu * a_u;
562  u_new = a_u;
563  m_newton( u_new, rhs, mu,
564  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
565  res_norm_a, res_norm_r, converged );
566 
567  if (converged) {
568  if (std::isfinite(u_new)) {
569  if (u_new > 0) {
570  step_success = true;
571  break;
572  }
573  }
574  }
575  dt *= myhalf;
576  }
577 
578  if (step_success) {
579 
580  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
581  a_u = u_new;
582  cur_time += dt;
583 
584  if (m_verbose) {
585  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
586  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
587  AMREX_DEVICE_PRINTF( " norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
588  res_norm_a, res_norm_r,
589  (converged ? "yes" : "no") );
590  }
591  if (snorm < m_stol) {
592  break;
593  }
594 
595  } else {
596 
597  a_success = false;
598  break;
599 
600  }
601 
602  n_step++;
603  }
604 
605  return;
606  }
607 
608  /*! \brief 2nd-order implicit Crank-Nicolson method */
609  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
610  void cn ( RT& a_u, /*!< solution */
611  bool& a_success /*!< success/failure flag */ ) const
612  {
613  RT cur_time = zero;
614  a_success = true;
615 
616  int n_step = 0;
617  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
618 
619  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
620  RT dt = m_cfl / std::sqrt(tau*tau);
621  if ((cur_time + dt) > m_t_final) {
622  dt = m_t_final - cur_time;
623  }
624  if (!std::isfinite(dt)) {
625  a_success = false;
626  break;
627  }
628 
629  RT res_norm_a = DBL_MAX;
630  RT res_norm_r = DBL_MAX;
631  bool converged = false;
632 
633  RT u_new = zero;
634  bool step_success = false;
635  while (!step_success) {
636 
637  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
638  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
639  break;
640  }
641  RT mu = one / (myhalf*dt);
642 
643  RT u1 = a_u;
644  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
645 
646  RT u2 = u1;
647  RT rhs = mu * (a_u + myhalf*dt*f1);
648  m_newton( u2, rhs, mu,
649  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
650  res_norm_a, res_norm_r, converged );
651  if (u2 <= 0) {
652  dt *= myhalf;
653  continue;
654  }
655  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
656 
657  u_new = a_u + myhalf * dt * (f1 + f2);
658 
659  if (converged) {
660  if (std::isfinite(u_new)) {
661  if (u_new > 0) {
662  step_success = true;
663  break;
664  }
665  }
666  }
667  dt *= myhalf;
668  }
669 
670  if (step_success) {
671 
672  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
673  a_u = u_new;
674  cur_time += dt;
675 
676  if (m_verbose) {
677  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
678  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
679  AMREX_DEVICE_PRINTF( " norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
680  res_norm_a, res_norm_r,
681  (converged ? "yes" : "no") );
682  }
683  if (snorm < m_stol) {
684  break;
685  }
686 
687  } else {
688 
689  a_success = false;
690  break;
691 
692  }
693 
694  n_step++;
695  }
696 
697  return;
698  }
699 
700  /*! \brief 2nd-order, 2-stage diagonally-implicit RK method */
701  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
702  void dirk212 ( RT& a_u, /*!< solution */
703  bool& a_success /*!< success/failure flag */ ) const
704  {
705  RT cur_time = zero;
706  a_success = true;
707 
708  int n_step = 0;
709  while ((cur_time < m_t_final) && (n_step < m_max_steps)) {
710 
711  RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
712  RT dt = m_cfl / std::sqrt(tau*tau);
713  if ((cur_time + dt) > m_t_final) {
714  dt = m_t_final - cur_time;
715  }
716  if (!std::isfinite(dt)) {
717  a_success = false;
718  break;
719  }
720 
721  bool converged = false;
722  RT res_norm_a = DBL_MAX;
723  RT res_norm_r = DBL_MAX;
724 
725  RT u_new = zero;
726  bool step_success = false;
727  while (!step_success) {
728 
729  if ( (dt < (amrex::Real(1.0e-12)*m_cfl/std::sqrt(tau*tau)))
730  && (dt < (amrex::Real(1.0e-12)*m_t_final)) ) {
731  break;
732  }
733  RT mu = one / dt;
734 
735  converged = true;
736  res_norm_a = zero;
737  res_norm_r = zero;
738 
739  RT u1 = a_u;
740  {
741  RT rhs = mu * a_u;
742  RT res_norm_a_i = DBL_MAX;
743  RT res_norm_r_i = DBL_MAX;
744  bool converged_i = false;
745  m_newton( u1, rhs, mu,
746  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
747  res_norm_a_i, res_norm_r_i, converged_i );
748  converged = converged && converged_i;
749  res_norm_a = std::max(res_norm_a, res_norm_a_i);
750  res_norm_r = std::max(res_norm_r, res_norm_r_i);
751  }
752  if (u1 <= 0) {
753  dt *= myhalf;
754  continue;
755  }
756  RT f1 = m_ode.rhs_func(u1, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
757 
758  RT u2 = u1;
759  {
760  RT rhs = mu * (a_u - dt*f1);
761  RT res_norm_a_i = DBL_MAX;
762  RT res_norm_r_i = DBL_MAX;
763  bool converged_i = false;
764  m_newton( u2, rhs, mu,
765  m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s,
766  res_norm_a_i, res_norm_r_i, converged_i );
767  converged = converged && converged_i;
768  res_norm_a = std::max(res_norm_a, res_norm_a_i);
769  res_norm_r = std::max(res_norm_r, res_norm_r_i);
770  }
771  if (u2 <= 0) {
772  dt *= myhalf;
773  continue;
774  }
775  RT f2 = m_ode.rhs_func(u2, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s);
776 
777  u_new = a_u + myhalf * dt * (f1 + f2);
778 
779  if (converged) {
780  if (std::isfinite(u_new)) {
781  if (u_new > 0) {
782  step_success = true;
783  break;
784  }
785  }
786  }
787  dt *= myhalf;
788  }
789 
790  if (step_success) {
791 
792  RT snorm = std::sqrt((a_u-u_new)*(a_u-u_new)/(a_u*a_u));
793  a_u = u_new;
794  cur_time += dt;
795 
796  if (m_verbose) {
797  AMREX_DEVICE_PRINTF( "Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n",
798  cur_time, dt, dt * std::sqrt(tau*tau), std::sqrt(a_u), snorm);
799  AMREX_DEVICE_PRINTF( " norms = %1.3e (abs), %1.3e (rel), converged = %s\n",
800  res_norm_a, res_norm_r,
801  (converged ? "yes" : "no") );
802  }
803  if (snorm < m_stol) {
804  break;
805  }
806 
807  } else {
808 
809  a_success = false;
810  break;
811 
812  }
813 
814  n_step++;
815  }
816 
817  return;
818  }
819 
820  };
821 }
822 
823 #endif
824 #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