src/predictor-corrector.h

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
  
// 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[];
    }
    boundary (output);
  }

  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;
  }

  trace
  void run()
  {
    t = 0., iter = 0;
    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();
  }