1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
| // Generic predictor/corrector time-integration
#include "utils.h"
// Required from solver
// fields updated by time-integration
extern scalar * evolving;
// how to compute updates
double (* update) (scalar * evolving, scalar * updates, double dtmax) = NULL;
// User-provided functions
// gradient
double (* gradient) (double, double, double) = minmod2;
// the timestep
double dt = 0.;
trace
static void advance_generic (scalar * output, scalar * input, scalar * updates,
double dt)
{
if (input != output)
trash (output);
foreach() {
scalar o, i, u;
for (o,i,u in output,input,updates)
o[] = i[] + dt*u[];
}
}
static void (* advance) (scalar * output, scalar * input, scalar * updates,
double dt) = advance_generic;
event defaults (i = 0)
{
// limiting
for (scalar s in all)
s.gradient = gradient;
display ("box();");
}
trace
void run()
{
t = 0., iter = 0; dimensional (dt == DT);
init_grid (N);
// main loop
perf.nc = perf.tnc = 0;
perf.gt = timer_start();
while (events (true)) {
// list of updates
scalar * updates = list_clone (evolving);
dt = dtnext (update (evolving, updates, DT));
if (gradient != zero) {
/* 2nd-order time-integration */
scalar * predictor = list_clone (evolving);
/* predictor */
advance (predictor, evolving, updates, dt/2.);
/* corrector */
update (predictor, updates, dt);
delete (predictor);
free (predictor);
}
advance (evolving, evolving, updates, dt);
delete (updates);
free (updates);
update_perf();
iter = inext, t = tnext;
}
timer_print (perf.gt, iter, perf.tnc);
free_grid();
}
|
