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