19 const double c1 = 37.0/378.0;
20 const double c2 = 0.0;
21 const double c3 = 250.0/621.0;
22 const double c4 = 125.0/594.0;
23 const double c5 = 0.0;
24 const double c6 = 512.0/1771.0;
26 const double cs1 = 2825.0/27648.0;
27 const double cs2 = 0.0;
28 const double cs3 = 18575.0/48384.0;
29 const double cs4 = 13525.0/55296.0;
30 const double cs5 = 277.0/14336.0;
31 const double cs6 = 0.25;
33 const double b21 = 0.2;
34 const double b31 = 3.0/40.0 , b32 = 9.0/40.0;
35 const double b41 = 0.3 , b42 = -0.9 , b43 = 1.2;
36 const double b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 , b54 = 35.0/27.0;
37 const double b61 = 1631.0/55296.0, b62 = 175.0/512.0, b63 = 575.0/13824.0, b64 = 44275.0/110592.0, b65 = 253.0/4096.0;
40 double t, tend, dt_shrink, dt_grow;
51 for (nv = 0; nv <
NVAR; nv++) v1[nv] = v0[nv];
59 FOR_EACH(nv,0,vars) v1[nv] = v0[nv] + dt*b21*k1[nv];
64 FOR_EACH(nv,0,vars) v1[nv] = v0[nv] + dt*(b31*k1[nv] + b32*k2[nv]);
70 v1[nv] = v0[nv] + dt*(b41*k1[nv] + b42*k2[nv] + b43*k3[nv]);
77 v1[nv] = v0[nv] + dt*(b51*k1[nv] + b52*k2[nv] + b53*k3[nv]
85 v1[nv] = v0[nv] + dt*(b61*k1[nv] + b62*k2[nv] + b63*k3[nv]
86 + b64*k4[nv] + b65*k5[nv]);
96 v5th[nv] = v0[nv] + dt*(c1*k1[nv] + c2*k2[nv] + c3*k3[nv]
97 + c4*k4[nv] + c5*k5[nv] + c6*k6[nv]);
103 v4th[nv] = v0[nv] + dt*(cs1*k1[nv] + cs2*k2[nv] + cs3*k3[nv]
104 + cs4*k4[nv] + cs5*k5[nv] + cs6*k6[nv]);
111 vscal[PRS] = fabs(v0[PRS]) + dt*fabs(k1[PRS]);
118 #if COOLING == MINEq && Fe_IONS > 0
119 vscal[FeI] = vscal[FeII] = vscal[FeIII] = 1.e18;
124 scrh = fabs(v5th[nv] - v4th[nv])/fabs(vscal[nv]);
125 err =
MAX(err, scrh);
132 err =
MAX(0.0, 1.e-18);
139 dt_grow = 0.9*dt*pow(err, -0.2);
140 dt_grow =
MIN(dt_grow, 5.0*dt);
145 if (fabs(t/tend - 1.0) < 1.e-9)
break;
149 if (dt > (tend - t)) dt = tend - t;
153 FOR_EACH(nv,0,vars) v0[nv] = v5th[nv];
157 print (
"! SolveODE_CK45: Number of substeps too large (%d)\n",ksub);
164 dt_shrink = 0.9*dt*pow(err, -0.25);
165 dt =
MAX(dt_shrink, 0.05*dt);
172 print (
"! SolveODE_CK45: number of substeps is %d\n", ksub);
void Radiat(double *, double *)
#define FOR_EACH(nv, beg, list)
void print(const char *fmt,...)
#define QUIT_PLUTO(e_code)