PLUTO
cooling_ode_solver.c
Go to the documentation of this file.
1 #include "pluto.h"
2 #if COOLING == MINEq
3  #include "cooling_defs.h"
4 #endif
5 
6 /* ********************************************************************* */
7 double SolveODE_CK45 (double *v0, double *k1, double *v5th,
8  double dt, double tol, intList *vars)
9 /*
10  *
11  * use explicit Cash-Karp 4-5 integrator
12  *
13  *
14  *
15  *********************************************************************** */
16 {
17  int nv, ksub, kfail;
18 
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;
25 
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;
32 
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;
38 
39  double scrh, err;
40  double t, tend, dt_shrink, dt_grow;
41  double v1[NVAR], v4th[NVAR];
42  double k2[NVAR], k3[NVAR], k4[NVAR], k5[NVAR], k6[NVAR];
43  double vscal[NVAR];
44  double tsub[4096];
45 
46 /* -------------------------------------------
47  copy ALL variables so that density is
48  defined when calling Radiat.
49  ------------------------------------------- */
50 
51  for (nv = 0; nv < NVAR; nv++) v1[nv] = v0[nv];
52 
53  t = tsub[0] = 0.0;
54  tend = dt;
55 
56  ksub = kfail = 0;
57  for (;;){
58 
59  FOR_EACH(nv,0,vars) v1[nv] = v0[nv] + dt*b21*k1[nv];
60 
61  /* -- get K2 -- */
62 
63  Radiat(v1, k2);
64  FOR_EACH(nv,0,vars) v1[nv] = v0[nv] + dt*(b31*k1[nv] + b32*k2[nv]);
65 
66  /* -- get K3 -- */
67 
68  Radiat(v1, k3);
69  FOR_EACH(nv, 0, vars){
70  v1[nv] = v0[nv] + dt*(b41*k1[nv] + b42*k2[nv] + b43*k3[nv]);
71  }
72 
73  /* -- get K4 -- */
74 
75  Radiat(v1, k4);
76  FOR_EACH(nv,0, vars) {
77  v1[nv] = v0[nv] + dt*(b51*k1[nv] + b52*k2[nv] + b53*k3[nv]
78  + b54*k4[nv]);
79  }
80 
81  /* -- get K5 -- */
82 
83  Radiat(v1, k5);
84  FOR_EACH(nv,0,vars){
85  v1[nv] = v0[nv] + dt*(b61*k1[nv] + b62*k2[nv] + b63*k3[nv]
86  + b64*k4[nv] + b65*k5[nv]);
87  }
88 
89  /* -- get K6 -- */
90 
91  Radiat(v1, k6);
92 
93  /* -- 5th order solution -- */
94 
95  FOR_EACH(nv,0,vars){
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]);
98  }
99 
100  /* -- 4th order embedded solution -- */
101 
102  FOR_EACH(nv,0,vars){
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]);
105  }
106 
107  /* -----------------------------------
108  compute error
109  ----------------------------------- */
110 
111  vscal[PRS] = fabs(v0[PRS]) + dt*fabs(k1[PRS]);
112 
113  NIONS_LOOP(nv) vscal[nv] = 1.0; /* fabs(v0[nv]) + dt*fabs(k1[nv]) + 1.e-6; */
114 
115 
116  /* -- do not take error on Fe -- */
117 
118  #if COOLING == MINEq && Fe_IONS > 0
119  vscal[FeI] = vscal[FeII] = vscal[FeIII] = 1.e18;
120  #endif
121 
122  err = 0.0;
123  FOR_EACH(nv,0,vars){
124  scrh = fabs(v5th[nv] - v4th[nv])/fabs(vscal[nv]);
125  err = MAX(err, scrh);
126  }
127  err /= tol;
128 
129  if (err < 1.0){ /* -- ok, accept step -- */
130 
131  ksub++;
132  err = MAX(0.0, 1.e-18);
133 
134  t += dt;
135  tsub[ksub] = t;
136 
137  /* -- provide an estimate for next dt -- */
138 
139  dt_grow = 0.9*dt*pow(err, -0.2);
140  dt_grow = MIN(dt_grow, 5.0*dt); /* -- do not increase more than 5 -- */
141  dt = dt_grow;
142 
143  /* -- exit loop if final time has been reached -- */
144 
145  if (fabs(t/tend - 1.0) < 1.e-9) break;
146 
147  /* -- adjust dt not to exceed exactly tend -- */
148 
149  if (dt > (tend - t)) dt = tend - t;
150 
151  /* -- initialize solution vector before continuing -- */
152 
153  FOR_EACH(nv,0,vars) v0[nv] = v5th[nv];
154  Radiat(v0, k1);
155 
156  if (ksub > 1000) {
157  print ("! SolveODE_CK45: Number of substeps too large (%d)\n",ksub);
158  QUIT_PLUTO(1);
159  }
160 
161  }else{ /* -- shrink dt and redo time step -- */
162 
163  kfail++;
164  dt_shrink = 0.9*dt*pow(err, -0.25);
165  dt = MAX(dt_shrink, 0.05*dt); /* -- do not decrease more than 20 -- */
166  }
167 
168  } /* -- end while -- */
169 
170  if (ksub > 100) {
171  int i;
172  print ("! SolveODE_CK45: number of substeps is %d\n", ksub);
173  }
174 
175  return (dt);
176 }
177 
178 /* ********************************************************************* */
179 double SolveODE_RK4 (double *v0, double *k1, double *v4th,
180  double dt, intList *var_list)
181 /*!
182  * Solve the system of ODE with a standard RK4 integrator (no
183  * adaptive step).
184  *
185  *
186  *
187  *
188  *********************************************************************** */
189 {
190  int nv;
191  double v1[NVAR];
192  double k2[NVAR], k3[NVAR], k4[NVAR];
193 
194 /* -------------------------------------------
195  copy ALL variables so that density is
196  defined when calling Radiat.
197  ------------------------------------------- */
198 
199  for (nv = 0; nv < NVAR; nv++) v1[nv] = v0[nv];
200 
201 /* -- step 1 -- */
202 
203  FOR_EACH(nv,0,var_list) v1[nv] = v0[nv] + 0.5*dt*k1[nv];
204 
205 /* -- step 2 -- */
206 
207  Radiat(v1, k2);
208  FOR_EACH(nv,0,var_list) v1[nv] = v0[nv] + 0.5*dt*k2[nv];
209 
210 /* -- step 3 -- */
211 
212  Radiat(v1, k3);
213  FOR_EACH(nv,0,var_list) v1[nv] = v0[nv] + 0.5*dt*k3[nv];
214 
215 /* -- step 4 -- */
216 
217  Radiat(v1, k4);
218  FOR_EACH(nv,0,var_list) {
219  v4th[nv] = v0[nv] + dt*(k1[nv] + 2.0*(k2[nv] + k3[nv]) + k4[nv])/6.0;
220  }
221 
222  return (0.0);
223 }
224 
225 /* ********************************************************************* */
226 double SolveODE_RKF12 (double *v0, double *k1, double *v2nd,
227  double dt, double tol, intList *vars)
228 /*
229  *
230  *
231  *
232  *********************************************************************** */
233 {
234  int nv;
235  double err, scrh, dt_grow;
236  double v1[NVAR], v1st[NVAR], vscal[NVAR];
237  double k2[NVAR];
238 
239 /* -------------------------------------------
240  copy ALL variables so that density is
241  defined when calling Radiat.
242  ------------------------------------------- */
243 
244  for (nv = 0; nv < NVAR; nv++) v1[nv] = v0[nv];
245 
246 /* -- Get K2 -- */
247 
248  FOR_EACH(nv, 0, vars) v1[nv] = v0[nv] + 0.5*dt*k1[nv];
249 
250  Radiat(v1, k2);
251 
252 /* -- 2nd order solution -- */
253 
254  FOR_EACH(nv, 0, vars) v2nd[nv] = v0[nv] + dt*k2[nv];
255 
256 /* -- 1st order solution -- */
257 
258  FOR_EACH(nv, 0, vars) v1st[nv] = v0[nv] + dt*k1[nv];
259 
260 /* -----------------------------------
261  compute error
262  ----------------------------------- */
263 
264  vscal[PRS] = fabs(v0[PRS]) + dt*fabs(k1[PRS]);
265 
266  NIONS_LOOP(nv) vscal[nv] = 1.0; /* fabs(v0[nv]) + dt*fabs(k1[nv]) + 1.e-6; */
267 
268 
269  err = 0.0;
270  FOR_EACH(nv,0,vars){
271  scrh = fabs(v2nd[nv] - v1st[nv])/fabs(vscal[nv]);
272  err = MAX(err, scrh);
273  }
274  err /= tol;
275 
276  if (err < 1.0) return ( 1.0);
277  else return (-1.0);
278 }
279 
280 /* ********************************************************************* */
281 double SolveODE_RKF23 (double *v0, double *k1, double *v3rd,
282  double dt, double tol, intList *vars)
283 /*
284  *
285  *
286  * y(3rd) = y(0) + dt*(k1 + k2 + 4*k3)/6
287  *
288  * k1 = RHS(y(0))
289  * k2 = RHS(y(0) + dt*k1)
290  * k3 = RHS(y(0) + 0.25*dt*(k1 + k2))
291  *
292  * y(2nd) = y(0) + dt*(k1 + k2)/2
293  *
294  *********************************************************************** */
295 {
296  int nv;
297  double err, scrh, dt_4, dt_6;
298  double v1[NVAR], v2nd[NVAR], vscal[NVAR];
299  double k2[NVAR], k3[NVAR];
300 
301  dt_4 = 0.25*dt;
302  dt_6 = dt/6.0;
303 
304 /* -------------------------------------------
305  copy ALL variables so that density is
306  defined when calling Radiat.
307  ------------------------------------------- */
308 
309  for (nv = 0; nv < NVAR; nv++) v1[nv] = v0[nv];
310 
311 /* -- Get K2 -- */
312 
313  FOR_EACH(nv,0,vars) v1[nv] = v0[nv] + dt*k1[nv];
314 
315  Radiat(v1, k2);
316 
317 /* -- Get K3 -- */
318 
319  FOR_EACH(nv,0,vars) v1[nv] = v0[nv] + dt_4*(k1[nv] + k2[nv]);
320 
321  Radiat(v1, k3);
322 
323 /* -- 3rd order solution -- */
324 
325  FOR_EACH(nv,0,vars) v3rd[nv] = v0[nv] + dt_6*(k1[nv] + k2[nv] + 4.0*k3[nv]);
326 
327 /* -- 2nd order solution -- */
328 
329  FOR_EACH(nv,0,vars) v2nd[nv] = v0[nv] + 0.5*dt*(k1[nv] + k2[nv]);
330 
331 /* -----------------------------------
332  compute error
333  ----------------------------------- */
334 
335  vscal[PRS] = fabs(v0[PRS]) + dt*fabs(k1[PRS]);
336 
337  NIONS_LOOP(nv) vscal[nv] = 1.0; /* fabs(v0[nv]) + dt*fabs(k1[nv]) + 1.e-6; */
338 
339 
340  err = 0.0;
341  FOR_EACH(nv,0,vars){
342  scrh = fabs(v3rd[nv] - v2nd[nv])/fabs(vscal[nv]);
343  err = MAX(err, scrh);
344  }
345  err /= tol;
346 
347  if (err < 1.0) return ( 1.0);
348  else return (-1.0);
349 }
350 
351 #define GAM (1.0/2.0)
352 #define A21 2.0
353 #define A31 (48.0/25.0)
354 #define A32 (6.0/25.0)
355 #define C21 -8.0
356 #define C31 (372.0/25.0)
357 #define C32 (12.0/5.0)
358 #define C41 (-112.0/125.0)
359 #define C42 (-54.0/125.0)
360 #define C43 (-2.0/5.0)
361 #define B1 (19.0/9.0)
362 #define B2 (1.0/2.0)
363 #define B3 (25.0/108.0)
364 #define B4 (125.0/108.0)
365 #define E1 (17.0/54.0)
366 #define E2 (7.0/36.0)
367 #define E3 0.0
368 #define E4 (125.0/108.0)
369 
370 #define C1X (1.0/2.0)
371 #define C2X (-3.0/2.0)
372 #define C3X (121.0/50.0)
373 #define C4X (29.0/250.0)
374 #define A2X 1.0
375 #define A3X (3.0/5.0)
376 
377 /* ********************************************************************* */
378 double SolveODE_ROS34 (double *v0, double *k1, double *v4th,
379  double dt, double tol)
380 /*
381  *
382  * Solve using Semi-Implicit Rosenbrock Method
383  *
384  *********************************************************************** */
385 {
386  int i, j, n, nv, ksub, kfail;
387  static int *indx;
388  double err, scrh, t, tend;
389  double v1[NVAR], vscal[NVAR], k2[NVAR], k3[NVAR];
390  double tsub[4096], dt_grow, dt_shrink;
391  static double **a, **J, **J2;
392  static double *g1, *g2, *g3, *g4;
393 
394 double vbeg[NVAR];
395 
396 /* -------------------------------------------
397  copy ALL variables so that density is
398  defined when calling Radiat.
399  ------------------------------------------- */
400 
401  for (nv = 0; nv < NVAR; nv++) v1[nv] = vbeg[nv] = v0[nv];
402 
403  n = NIONS + 1;
404 
405  if (indx == NULL){
406  indx = ARRAY_1D (n, int);
407  a = ARRAY_2D (n, n, double);
408  J = ARRAY_2D (n, n, double);
409  J2 = ARRAY_2D (n, n, double);
410  g1 = ARRAY_1D (n, double);
411  g2 = ARRAY_1D (n, double);
412  g3 = ARRAY_1D (n, double);
413  g4 = ARRAY_1D (n, double);
414  }
415 
416  for (i = 0; i < n - 1; i++) vscal[i + NFLX] = 1.0;
417 
418  /* -- do not take error on Fe -- */
419 
420  #if COOLING == MINEq && Fe_IONS > 0
421  vscal[FeI] = vscal[FeII] = vscal[FeIII] = 1.e18;
422  #endif
423 
424  vscal[PRS] = fabs(v0[PRS]);
425 
426  Jacobian (v0, k1, J);
427 /*
428  Numerical_Jacobian (v0, J2);
429 {
430  int k,l;
431  double T0;
432 
433  T0 = v0[PRS]/v0[RHO]*KELVIN*MeanMolecularWeight(v0);
434  printf (" T0 = %12.6e\n",T0);
435 
436  for (k = 0; k < n; k++){
437 
438  for (l = 0; l < n; l++){
439  printf ("%4.1f ", 100.*fabs(J2[k][l] - J[k][l])/(fabs(J[k][l]) + 1.e-12));
440  }
441 
442 // for (l = 0; l < 1; l++){
443 // printf ("%12.6e %12.6e", J2[k][l], J[k][l]);
444 // }
445  printf ("\n");
446  }
447  exit(1);
448 }
449 */
450  t = 0.0;
451  tend = dt;
452 
453  ksub = kfail = 0;
454  for (;;){
455 
456  /* -- Compute (I - hJ) -- */
457 
458  for (i = 0; i < n; i++) {
459  for (j = 0; j < n; j++) a[i][j] = -J[i][j];
460  a[i][i] += 1.0/(GAM*dt);
461  }
462  LUDecompose (a, n, indx, &scrh); /* LU decomposition of the matrix. */
463 
464  /* -- set right hand side for g1 -- */
465 
466  for (i = 0; i < n - 1; i++) {
467  g1[i] = k1[i + NFLX];
468  }
469  g1[n - 1] = k1[PRS];
470 
471  /* -- solve for g1 -- */
472 
473  LUBackSubst (a, n, indx, g1);
474  for (i = 0; i < n - 1; i++) {
475  v1[i + NFLX] = v0[i + NFLX] + A21*g1[i];
476  }
477  v1[PRS] = v0[PRS] + A21*g1[n - 1];
478 
479  Radiat (v1, k2);
480 
481  /* -- set right hand side for g2 -- */
482 
483  for (i = 0; i < n - 1; i++) {
484  g2[i] = k2[i + NFLX] + C21*g1[i]/dt;
485  }
486  g2[n - 1] = k2[PRS] + C21*g1[n - 1]/dt;
487 
488  /* -- solve for g2 -- */
489 
490  LUBackSubst (a, n, indx, g2);
491  for (i = 0; i < n - 1; i++) {
492  v1[i + NFLX] = v0[i + NFLX] + A31*g1[i] + A32*g2[i];
493  }
494  v1[PRS] = v0[PRS] + A31*g1[n - 1] + A32*g2[n - 1];
495 
496  Radiat (v1, k3);
497 
498  /* -- set right hand side for g3 -- */
499 
500  for (i = 0; i < n - 1; i++) {
501  g3[i] = k3[i + NFLX] + (C31*g1[i] + C32*g2[i])/dt;
502  }
503  g3[n - 1] = k3[PRS] + (C31*g1[n - 1] + C32*g2[n - 1])/dt;
504 
505  /* -- solve for g3 -- */
506 
507  LUBackSubst (a, n, indx, g3);
508 
509  /* -- set right hand side for g4 -- */
510 
511  for (i = 0; i < n - 1; i++) {
512  g4[i] = k3[i + NFLX] + (C41*g1[i] + C42*g2[i] + C43*g3[i])/dt;
513  }
514  g4[n - 1] = k3[PRS] + (C41*g1[n - 1] + C42*g2[n - 1] + C43*g3[n - 1])/dt;
515 
516  /* -- solve for g4 -- */
517 
518  LUBackSubst (a, n, indx, g4);
519 
520  /* -- 4th order solution & error estimation -- */
521 
522  i = n - 1;
523  v4th[PRS] = v0[PRS] + B1*g1[i] + B2*g2[i] + B3*g3[i] + B4*g4[i];
524  err = fabs(E1*g1[i] + E2*g2[i] + E3*g3[i] + E4*g4[i])/vscal[PRS];
525 
526  for (i = 0; i < n - 1; i++) {
527  v4th[i + NFLX] = v0[i + NFLX] + B1*g1[i] + B2*g2[i] + B3*g3[i] + B4*g4[i];
528  scrh = E1*g1[i] + E2*g2[i] + E3*g3[i] + E4*g4[i];
529  err = MAX(err, fabs(scrh)/vscal[i + NFLX]);
530  }
531 
532  err /= tol;
533 
534  if (err < 1.0) {
535 
536  ksub++;
537  err = MAX(0.0, 1.e-18);
538 
539  t += dt;
540  tsub[ksub] = t;
541 
542  /* -- provide an estimate for next dt -- */
543 
544  dt_grow = 0.9*dt*pow(err, -0.25);
545  dt_grow = MIN(dt_grow, 5.0*dt);
546  dt = dt_grow;
547 
548  /* -- exit loop if final time has been reached -- */
549 
550  if (fabs(t/tend - 1.0) < 1.e-9) break;
551 
552  /* -- adjust dt not to exceed tend -- */
553 
554  if (dt > (tend - t)) dt = tend - t;
555 
556  /* -- initialize solution vector continuing -- */
557 
558  v0[PRS] = v4th[PRS];
559  NIONS_LOOP(nv) v0[nv] = v4th[nv];
560  Radiat (v0, k1);
561  Jacobian (v0, k1, J);
562 
563  if (ksub > 1000){
564  print ("! SolveODE_ROS34: Number of substeps too large (%d)\n",ksub);
565  QUIT_PLUTO(1);
566  }
567 
568  }else{ /* -- shrink dt and redo time step -- */
569 
570  kfail++;
571  dt_shrink = 0.9*dt*pow(err, -1.0/3.0);
572  dt = MAX(dt_shrink, 0.05*dt);
573  }
574  }
575 
576  if (ksub > 2) {
577  int i;
578  print ("! SolveODE_ROS34: number of substeps is %d\n", ksub);
579 /*
580  for (i = 1; i <= ksub; i++){
581  printf ("# %d, dt = %12.6e, t/dt = %f, tend = %12.6e\n",
582  i, tsub[i]-tsub[i-1],tsub[i]/tend, tend);
583  }
584  printf ("kfail = %d\n",kfail);
585  QUIT_PLUTO(1);
586 */
587  }
588 
589  return (dt);
590 }
591 
592 #undef GAM
593 #undef A21
594 #undef A31
595 #undef A32
596 #undef C21
597 #undef C31
598 #undef C32
599 #undef C41
600 #undef C42
601 #undef C43
602 #undef B1
603 #undef B2
604 #undef B3
605 #undef B4
606 #undef E1
607 #undef E2
608 #undef E3
609 #undef E4
610 #undef C1X
611 #undef C2X
612 #undef C3X
613 #undef C4X
614 #undef A2X
615 #undef A3X
616 
617 
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
NIONS_LOOP
#define NIONS_LOOP(n)
Definition: pluto.h:614
E4
#define E4
Definition: cooling_ode_solver.c:368
n
static int n
Definition: analysis.c:3
Radiat
void Radiat(double *, double *)
Definition: radiat.c:94
configure.scrh
tuple scrh
Definition: configure.py:200
E3
#define E3
Definition: cooling_ode_solver.c:367
E1
#define E1
Definition: cooling_ode_solver.c:365
B2
#define B2
Definition: cooling_ode_solver.c:362
FOR_EACH
#define FOR_EACH(nv, beg, list)
Definition: macros.h:89
C41
#define C41
Definition: cooling_ode_solver.c:358
SolveODE_CK45
double SolveODE_CK45(double *v0, double *k1, double *v5th, double dt, double tol, intList *vars)
Definition: cooling_ode_solver.c:7
MIN
#define MIN(a, b)
Definition: macros.h:104
cooling_defs.h
B1
#define B1
Definition: cooling_ode_solver.c:361
SolveODE_RKF12
double SolveODE_RKF12(double *v0, double *k1, double *v2nd, double dt, double tol, intList *vars)
Definition: cooling_ode_solver.c:226
LUBackSubst
void LUBackSubst(double **a, int n, int *indx, double b[])
Definition: math_lu_decomp.c:84
NFLX
#define NFLX
Definition: mod_defs.h:32
C43
#define C43
Definition: cooling_ode_solver.c:360
SolveODE_RKF23
double SolveODE_RKF23(double *v0, double *k1, double *v3rd, double dt, double tol, intList *vars)
Definition: cooling_ode_solver.c:281
NIONS
#define NIONS
Definition: cooling.h:28
j
int j
Definition: analysis.c:2
A32
#define A32
Definition: cooling_ode_solver.c:354
Jacobian
void Jacobian(real *v, real *rhs, real **dfdy)
Definition: jacobian.c:4
E2
#define E2
Definition: cooling_ode_solver.c:366
C32
#define C32
Definition: cooling_ode_solver.c:357
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
A31
#define A31
Definition: cooling_ode_solver.c:353
SolveODE_ROS34
double SolveODE_ROS34(double *v0, double *k1, double *v4th, double dt, double tol)
Definition: cooling_ode_solver.c:378
C31
#define C31
Definition: cooling_ode_solver.c:356
B3
#define B3
Definition: cooling_ode_solver.c:363
INT_LIST
Definition: structs.h:328
A21
#define A21
Definition: cooling_ode_solver.c:352
pluto.h
PLUTO main header file.
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
i
int i
Definition: analysis.c:2
C21
#define C21
Definition: cooling_ode_solver.c:355
GAM
#define GAM
Definition: cooling_ode_solver.c:351
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
LUDecompose
int LUDecompose(double **a, int n, int *indx, double *d)
Definition: math_lu_decomp.c:12
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
B4
#define B4
Definition: cooling_ode_solver.c:364
SolveODE_RK4
double SolveODE_RK4(double *v0, double *k1, double *v4th, double dt, intList *var_list)
Definition: cooling_ode_solver.c:179
C42
#define C42
Definition: cooling_ode_solver.c:359