sk_solv.c

   1 /* Apache 2.0 INS-AMU 2015 */
   2
   3 #include <string.h>
   4 #include <stdlib.h>
   5
   6 #include "sk_solv.h"
   7 #include "sk_util.h"
   8
   9
  10 int sk_solv_init(
  11         struct sk_solv *s,
  12         sk_sys sys, void *sys_data,
  13         sk_sch scheme, void *scheme_data,
  14         sk_out out, void *out_data,
  15         sk_hist_filler hf,
  16         int seed,
  17         int nx, double *x0,
  18         int nc, int *vi, double *vd,
  19         double t0, double dt
  20         )
  21 {
  22         s->nx = nx;
  23         s->nc = nc;
  24         s->cont = 1;
  25         s->sys = sys;
  26         s->sysd = sys_data;
  27         s->sch = scheme;
  28         s->schd = scheme_data;
  29         s->out = out;
  30         s->outd = out_data;
  31         s->hf = hf;
  32         SK_MALLOCHECK(s->x = malloc (sizeof(double) * nx));
  33         if (nc > 0 && vi!=NULL && vd!=NULL) {
  34                 sk_hist_init(&s->hist, nc, vi, vd, t0, dt);
  35                 sk_hist_fill(&s->hist, hf);
  36                 SK_MALLOCHECK(s->c = malloc(sizeof(double) * s->nc));
  37         } else {
  38                 s->nc = 0;
  39                 s->c = NULL;
  40         }
  41         rk_seed(seed, &(s->rng));
  42         memcpy(s->x, x0, sizeof(double) * s->nx);
  43         s->t = t0;
  44         s->dt = dt;
  45         return 0;
  46 }
  47
  48 void sk_solv_free(struct sk_solv *s)
  49 {
  50         free(s->x);
  51         if (s->c != NULL) {
  52                 free(s->c);
  53                 sk_hist_free(&s->hist);
  54         }
  55 }
  56
  57 int sk_solv_cont(struct sk_solv *s)
  58 {
  59         s->cont = 1;
  60         if (s->nc)
  61                 do {
  62                         s->t += s->dt;
  63                         sk_hist_get(&s->hist, s->t, s->c);
  64                         s->sch(s->schd, &(s->rng), s->sys, s->sysd, s->t, s->dt,
  65                                s->nx, s->x, s->nc, s->c);
  66                         sk_hist_set(&s->hist, s->t, s->x);
  67                         s->cont = s->out(s->outd, s->t, s->nx, s->x);
  68                 } while (s->cont);
  69         else
  70                 do {
  71                         s->t += s->dt;
  72                         s->sch(s->schd, &(s->rng), s->sys, s->sysd, s->t, s->dt,
  73                                s->nx, s->x, s->nc, s->c);
  74                         s->cont = s->out(s->outd, s->t, s->nx, s->x);
  75                 } while (s->cont);
  76
  77         return 0;
  78 }