src/sk_scheme.c

   1 /* Apache 2.0 INS-AMU 2015 */
   2
   3 #include <math.h>
   4
   5 #include <stdlib.h>
   6
   7 #include "sk_scheme.h"
   8 #include "sk_util.h"
   9
  10 int sk_sch_id_init(sk_sch_id_data *d, int nx)
  11 {
  12         SK_MALLOCHECK(d->f=malloc(sizeof(double)*nx));
  13         SK_MALLOCHECK(d->g=malloc(sizeof(double)*nx));
  14         SK_MALLOCHECK(d->z=malloc(sizeof(double)*nx));
  15         return 0;
  16 }
  17
  18 void sk_sch_id_free(sk_sch_id_data *d)
  19 {
  20         free(d->f);
  21         free(d->g);
  22         free(d->z);
  23 }
  24
  25 SK_DEFSCH(sk_sch_id) {
  26         int i;
  27         sk_sch_id_data *d;
  28         /* unused */ (void) dt;
  29         d = data;
  30         sk_hist_get(hist, t, c);              /*  Jf    Jg           Jc */
  31         (*sys)(sysd, hist, t, 0, nx, x, d->f, d->g, NULL, NULL, nc, c, NULL);
  32         sk_util_fill_gauss(rng, nx, d->z);
  33         for (i=0; i<nx; i++)
  34                 x[i] = d->f[i] + d->g[i] * d->z[i];
  35         sk_hist_set(hist, t, c);
  36         return 0;
  37 }
  38
  39 int sk_sch_em_init(sk_sch_em_data *d, int nx)
  40 {
  41         SK_MALLOCHECK(d->f=malloc(sizeof(double)*nx));
  42         SK_MALLOCHECK(d->g=malloc(sizeof(double)*nx));
  43         SK_MALLOCHECK(d->z=malloc(sizeof(double)*nx));
  44         return 0;
  45 }
  46
  47 void sk_sch_em_free(sk_sch_em_data *d)
  48 {
  49         free(d->f);
  50         free(d->g);
  51         free(d->z);
  52 }
  53
  54 SK_DEFSCH(sk_sch_em)
  55 {
  56         int i;
  57         double sqrt_dt;
  58         sk_sch_em_data *d = data;
  59         sk_hist_get(hist, t, c);             /*  Jf    Jg           Jc */
  60         (*sys)(sysd, hist, t, 0, nx, x, d->f, d->g, NULL, NULL, nc, c, NULL);
  61         sk_util_fill_gauss(rng, nx, d->z);
  62         sqrt_dt = sqrt(dt);
  63         for (i=0; i<nx; i++)
  64                 x[i] += dt * d->f[i] + sqrt_dt * d->g[i] * d->z[i];
  65         sk_hist_set(hist, t, c);
  66         return 0;
  67 }
  68
  69 int sk_sch_emcolor_init(sk_sch_emcolor_data *d, int nx, double lam)
  70 {
  71         SK_MALLOCHECK(d->f=malloc(sizeof(double)*nx));
  72         SK_MALLOCHECK(d->g=malloc(sizeof(double)*nx));
  73         SK_MALLOCHECK(d->z=malloc(sizeof(double)*nx));
  74         SK_MALLOCHECK(d->eps=malloc(sizeof(double)*nx));
  75         d->first_call = 1;
  76         d->lam = lam;
  77         return 0;
  78 }
  79
  80 void sk_sch_emcolor_free(sk_sch_emcolor_data *d)
  81 {
  82         free(d->f);
  83         free(d->g);
  84         free(d->z);
  85         free(d->eps);
  86 }
  87
  88 SK_DEFSCH(sk_sch_emcolor)
  89 {
  90         int i;
  91         double E; /* not stored so can be chaned while running */
  92         sk_sch_emcolor_data *d = data;
  93         if (d->first_call) {
  94                 sk_util_fill_gauss(rng, nx, d->z);      /*  Jf    Jg           Jc */
  95                 (*sys)(sysd, hist, t-dt, 0, nx, x, d->f, d->g, NULL, NULL, nc, c, NULL);
  96                 for (i=0; i<nx; i++)
  97                         d->eps[i] = sqrt(d->g[i] * d->lam) * d->z[i];
  98                 d->first_call = 0;
  99         }
 100         E = exp(-d->lam * dt);
 101         sk_util_fill_gauss(rng, nx, d->z);
 102         sk_hist_get(hist, t, c);             /*  Jf    Jg           Jc */
 103         (*sys)(sysd, hist, t, 0, nx, x, d->f, d->g, NULL, NULL, nc, c, NULL);
 104         for (i=0; i<nx; i++) {
 105                 x[i] += dt * (d->f[i] + d->eps[i]);
 106                 d->eps[i] *= E;
 107                 d->eps[i] += sqrt(d->g[i] * d->lam * (1 - E*E)) * d->z[i];
 108         }
 109         sk_hist_set(hist, t, c);
 110         return 0;
 111 }
 112
 113
 114 int sk_sch_heun_init(sk_sch_heun_data *d, int nx)
 115 {
 116         SK_MALLOCHECK(d->fl=malloc(sizeof(double)*nx));
 117         SK_MALLOCHECK(d->fr=malloc(sizeof(double)*nx));
 118         SK_MALLOCHECK(d->gl=malloc(sizeof(double)*nx));
 119         SK_MALLOCHECK(d->gr=malloc(sizeof(double)*nx));
 120         SK_MALLOCHECK(d->z=malloc(sizeof(double)*nx));
 121         SK_MALLOCHECK(d->xr=malloc(sizeof(double)*nx));
 122         return 0;
 123 }
 124
 125 void sk_sch_heun_free(sk_sch_heun_data *d)
 126 {
 127         free(d->fl);
 128         free(d->fr);
 129         free(d->gl);
 130         free(d->gr);
 131         free(d->z);
 132         free(d->xr);
 133 }
 134
 135 SK_DEFSCH(sk_sch_heun)
 136 {
 137         int i;
 138         double sqrt_dt;
 139         sk_sch_heun_data *d = data;
 140         /* predictor */
 141         sk_hist_get(hist, t, c);             /*    Jf    Jg           Jc */
 142         (*sys)(sysd, hist, t, 0, nx, x, d->fl, d->gl, NULL, NULL, nc, c, NULL);
 143         for (i=0; i<nx; i++)
 144                 d->xr[i] = x[i] + dt * d->fl[i];
 145         sk_hist_set(hist, t, c);
 146         /* corrector */
 147         sk_hist_get(hist, t+dt, c);                 /*    Jf    Jg           Jc */
 148         (*sys)(sysd, hist, t+dt, 0, nx, d->xr, d->fr, d->gr, NULL, NULL, nc, c, NULL);
 149         sk_util_fill_gauss(rng, nx, d->z);
 150         sqrt_dt = sqrt(dt);
 151         for (i=0; i<nx; i++)
 152                 x[i] += 0.5 * (dt*(d->fl[i] + d->fr[i])
 153                                 + sqrt_dt*(d->gl[i] + d->gr[i])*d->z[i]);
 154         /* TODO double check this has an effect */
 155         sk_hist_set(hist, t+dt, c);
 156         return 0;
 157 }