bench/bench_net_exc.c

   1 /* copyright 2016 Apache 2 sddekit authors */
   2 #include <stdlib.h>
   3 #include <stdio.h>
   4 #include <time.h>
   5
   6 #include "weights.c"
   7 #include "tract_lengths.c"
   8
   9 #include "sddekit.h"
  10
  11 typedef struct { sd_out *out; } outd;
  12
  13 static sd_stat out_igc(void *data, double t,
  14                              uint32_t nx, double * restrict x,
  15                              uint32_t nc, double * restrict c)
  16 {
  17         (void) nc; (void) c;
  18         outd *d = data;
  19         return (d->out)->apply(d->out, t, nx, x, 0, NULL);
  20 }
  21
  22 int main() {
  23         uint32_t i, n=76, nnz, *Or, *Ic;
  24         double *w=weights, *d=tract_lengths, *sw, *sd, *x0;
  25
  26         /* set up outputs
  27          *
  28          * sol -> igc -> tee -> tavg -> out file
  29          *                   \> until t = 1e3
  30          */
  31         sd_out_file *of = sd_out_file_new_from_name("bench_net_exc.dat");
  32         sd_out_tavg *tavg = sd_out_tavg_new(20, SD_AS(of, out));
  33         sd_out *until = sd_out_new_until(1e3);
  34         sd_out_tee *tee = sd_out_tee_new(2);
  35         tee->set_out(tee, 0, until);
  36         tee->set_out(tee, 1, SD_AS(tavg, out));
  37         outd igcd = { .out = SD_AS(tee, out) };
  38         sd_out *out = sd_out_new_cb(&igcd, &out_igc);
  39
  40         /* connectivity, assuming conduction velocity of 1.0 */
  41         /*
  42         sd_util_read_square_matrix("bench/conn76/weights.txt", &n, &w);
  43         sd_util_read_square_matrix("bench/conn76/tract_lengths.txt", &n, &d);
  44         */
  45         sd_sparse_from_dense(n, n, weights, tract_lengths, 0.0, &nnz, &Or, &Ic, &sw, &sd);
  46         fprintf(stdout, "[bench_net_exc] nnz=%d\n", nnz);
  47
  48         /* setup model */
  49         sd_sys_exc *exc = sd_sys_exc_new();
  50         exc->set_a(exc, 1.01);
  51         exc->set_tau(exc, 3.0);
  52         exc->set_D(exc, 0.01);
  53         exc->set_k(exc, 0.001);
  54         sd_net *net = sd_net_new_hom(n, SD_AS(exc, sys), 2, 1, 1, nnz, Or, Ic, w, d);
  55
  56         /* setup scheme & solver */
  57         x0 = sd_malloc(sizeof(double)*2*n);
  58         for (i=0; i<(2*n); i++)
  59                 x0[i] = 0.0;
  60         sd_sch *heun = sd_sch_new_heun(2*n);
  61         sd_hfill *hf = sd_hfill_new_val(0.0);
  62         sd_sol *sol = sd_sol_new_default(SD_AS(net, sys), heun, out, hf,
  63                         42, 2*n, x0, n, nnz, Ic, sd, 0.0, 0.01);
  64
  65         /* solve */
  66         clock_t tic, toc;
  67         tic = clock();
  68         sol->cont(sol);
  69         toc = clock();
  70         double reft = (double) (toc - tic) / CLOCKS_PER_SEC;
  71         sd_log_info("continuation required %.3f s\n", reft);
  72
  73         /* clean up */
  74         sd_free(x0);
  75         /*
  76         sd_free(w);
  77         sd_free(d);
  78         */
  79         sd_free(Or);
  80         sd_free(Ic);
  81         sd_free(sw);
  82         sd_free(sd);
  83
  84         return 0;
  85 }