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