openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
bench.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is licensed under the openCARP Academic Public License (APL)
7 // v1.0: You can use and redistribute it and/or modify it in non-commercial
8 // academic environments under the terms of APL as published by the openCARP
9 // project v1.0, or (at your option) any later version. Commercial use requires
10 // a commercial license (info@opencarp.org).
11 //
12 // This program is distributed without any warranty; see the openCARP APL for
13 // more details.
14 //
15 // You should have received a copy of the openCARP APL along with this program
16 // and can find it online: http://www.opencarp.org/license
17 // ----------------------------------------------------------------------------
18 
27 #define PrMGLOBAL
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include "cmdline.h"
31 #include <time.h>
32 #include <sys/time.h>
33 #include <sys/resource.h>
34 #include <unistd.h>
35 #include <limits.h>
36 #include <math.h>
37 #include <stdbool.h>
38 #include <vector>
39 #include <sstream>
40 #include <string>
41 #include <cmath>
42 #include "basics.h"
43 
44 #include "MULTI_ION_IF.h"
45 
46 #include "trace.h"
47 #include "clamp.h"
48 #include "ap_analyzer.h"
49 #include "restitute.h"
50 #include "stretch.h"
51 #include "sv_init.h"
52 #include "bench_utils.h"
53 #include "target.h"
54 
55 #include "libgen.h"
56 #include "build_info.h"
57 
58 #include "petsc_utils.h" //TODO for EXIT, FPRINT*
59 #include "SF_init.h" // for SF::init_xxx()
60 
61 
62 // globals
63 namespace opencarp {
64 namespace user_globals {
68  std::map<mesh_t, sf_mesh> mesh_reg;
70  std::map<SF::quadruple<int>, SF::index_mapping<int> > map_reg;
71  // /// the physics
72  // std::map<physic_t, Basic_physic*> physics_reg;
75  // /// important solution vectors from different physics
76  // std::map<datavec_t, sf_vec*> datavec_reg;
77 } // namespace user_globals
78 } // namespace opencarp
79 
80 using namespace opencarp;
81 using namespace limpet;
82 
83 #undef __FUNCT__
84 #define __FUNCT__ "main"
85 
86 int main(int argc, char *argv[]) {
87 
88  double dt = 10.e-3; // time step needed for the setup of the LUT's (ms)
89 
90  // specify the number of regions and their names
91  const int num_region = 1;
92  enum { R1 }
93  Region;
94 
95  // specifiy the ionic model for each region
96  IonType* RegionDef;
97 
98  // specify the plugins for each region
99  int num_plugins;
100  IonTypeList RegionPlug = {};
101 
102  // specify the region for each node
103  char *Regions;
104 
105 
106  struct gengetopt_args_info params;
107 
108  // we allocate the vectors on the heap since they will be freed by the MIIF
109  sf_vec* Vmv;
110  sf_vec* I_ion;
111 
112  FILE *fhdls[NUM_IMP_DATA_TYPES+1];
113  GVEC_DUMP gvd;
114  IOCtrl io;
115 
116  char *PetscDBfile = COMPAT_PETSC_NULLPTR; // could be used to pass in ~/.petscrc
117  char *help_msg = COMPAT_PETSC_NULLPTR;
118  initialize_PETSc(&argc, argv, PetscDBfile, help_msg);
119  COMPAT_PetscOptionsInsertString("-options_left no");
120 
121  if (cmdline_parser(argc, argv, &params) != 0)
122  exit(1);
123 
124  // performance measurements
125  double t1, t2, t3, t4, t5, t6, t7, t8;
126  struct timeval tt1;
127 
128  event_timing timings[N_TIMINGS];
129  initialize_timings(timings);
130 
131  // start measurementstrace_IA
132  get_time(t1);
133 
134  // assign command line params and depending variables
135  dt = params.dt_arg;
136  double i_stim = params.stim_curr_arg*dt;
137  double bcl = params.bcl_arg;
138  double stim_start = params.stim_start_arg;
139  double stim_dur = params.stim_dur_arg;
140  int nstim = params.numstim_arg ? params.numstim_arg : -1;
141  double dt_out = params.dt_out_arg;
142  double start_out = params.start_out_arg;
143  int numNode = params.num_arg;
144  int stim_assign = params.stim_assign_flag;
145  bool extVmUpdate = params.ext_vm_update_given;
146  bool analyze_AP = params.APstatistics_given || params.restitute_given;
147  bool restitute = params.restitute_given;
148  bool validate = params.validate_given;
149  bool APclamp = params.AP_clamp_file_given;
150  bool do_trace = restitute ? params.res_trace_flag : (!params.no_trace_flag);
151  TrgList stim_lst;
152  float duration = determine_duration(&params, &stim_lst);
153 
154  io.w2file = validate ? 1 : params.fout_given;
155  io.wbin = validate ? 1 : params.bin_flag;
156 
157  if(get_size() > numNode) {
158  numNode = get_size();
159 
160  log_msg(NULL, 3, 0, "Warning: Number of nodes is less than number of processes used!");
161  log_msg(NULL, 3, 0, "Setting number of nodes to %d.", numNode);
162  }
163 
164  // stimulate with trace from file
165  trace stim_trace;
166  bool tr_stim = false;
167  if (params.stim_file_given) {
168  char *fname = dupstr(params.stim_file_arg);
169  read_trace(&stim_trace, fname);
170  resample_trace(&stim_trace, dt);
171  tr_stim = true;
172  stim_dur = stim_trace.dur;
173  if (bcl < stim_trace.dur)
174  log_msg(NULL, 4, 0, "BCL less than stimulus duration!");
175  }
176 
177  // illuminate with trace from file
178  trace light_trace;
179  bool tr_light = false;
180  if( params.light_file_given ) {
181  char* fname = dupstr( params.light_file_arg );
182  read_trace( &light_trace, fname );
183  resample_trace( &light_trace, dt );
184  tr_light = true;
185 
186  int n_clipped_to_0 = 0;
187  for(int i=0; i<light_trace.N; i++) {
188  if( light_trace.s[i] < 0.0 ) {
189  light_trace.s[i] = 0.0;
190  n_clipped_to_0++;
191  }
192  }
193  if( n_clipped_to_0 ) {
194  FPRINTF( WORLD stderr, "WARNING: %d irradiance values in %s were negative (clipped to 0.0)\n", n_clipped_to_0, fname );
195  }
196  }
197 
198  // state variable clamp
199  Clamp *sv_cl;
200  int SVclamp = process_sv_clamps(params.clamp_SVs_arg, params.SV_clamp_files_arg, &sv_cl, dt);
201 
202  // action potential clamp
203  Clamp ap_cl;
204  if (APclamp){
205  const char *Vm_clamp = "Vm";
206  initialize_sv_clamp(&ap_cl, Vm_clamp, params.AP_clamp_file_arg, dt);
207  }
208 
209  // voltage clamp experiment
210  Clamp cl;
211  bool IsVclamp = initialize_clamp(&cl, params.clamp_arg, params.clamp_ini_arg,
212  params.clamp_start_arg, params.clamp_dur_arg,
213  params.clamp_file_arg,
214  !params.clamp_ini_given ? CLAMP_TM_IDX : 0, &duration);
215  if (IsVclamp) {
216  if (APclamp) {
217  log_msg(NULL, 0, 0, "info: cannot clamp voltage and have AP clamp at same time");
218  exit(1);
219  }
220  if (cl.transient) {
221  log_msg(NULL, 0, 0, "info: Vm will be clamped to %f mV over the interval [%f, %f]",
222  params.clamp_arg, params.clamp_start_arg, params.clamp_start_arg + params.clamp_dur_arg);
223  } else {
224  log_msg(NULL, 0, 0, "info: Vm will be clamped from %f to %f mV over the interval [%f, %f]",
225  params.clamp_ini_arg, params.clamp_arg, params.clamp_start_arg,
226  params.clamp_start_arg + params.clamp_dur_arg);
227  }
228  }
229 
230  action_potential AP;
231  if (analyze_AP) initialize_AP_analysis(&AP);
232 
233  Regions = static_cast<char *>(calloc(numNode, sizeof(char) ));
234 
235  SF::init_vector(&Vmv, numNode, -1);
236  SF::init_vector(&I_ion, Vmv);
237 
238  int *stim_list = static_cast<int *>(malloc(numNode*sizeof(int) ));
239  for (int i = 0; i < numNode; i++) {
240  Regions[i] = R1;
241  stim_list[i] = i;
242  }
243 
244  for (unsigned int i = 0; i < params.load_module_given; i++) {
245 #ifdef HAVE_DLOPEN
246  load_ionic_module(params.load_module_arg[i]);
247  if (!params.imp_given) {
248 
249  char *imp = strdup(basename(params.load_module_arg[i]) );
250  char *dot = strrchr(imp, '.');
251  if (dot)
252  *dot = 0;
253  params.imp_arg = imp;
254  }
255 #else // ifdef HAVE_DLOPEN
256  log_msg(NULL, 5, 0, "Compile with USE_DLOPEN to support dynamic module loading.");
257  EXIT(-1);
258 #endif // ifdef HAVE_DLOPEN
259  }
260 
261  if (get_ion_type(std::string(params.imp_arg)) == NULL) {
262  fprintf(stderr, "Illegal IMP specified: %s\n", params.imp_arg);
263  exit(1);
264  } else {
265  RegionDef = get_ion_type(std::string(params.imp_arg));
266  }
267  if (!get_plug_flag(params.plug_in_arg, &num_plugins, RegionPlug)) {
268  fprintf(stderr, "Illegal plugin specified: %s\n", params.plug_in_arg);
269  exit(1);
270  }
271  // Fetch target enum from command line string
272  Target target = get_target_from_string(std::string(params.target_arg));
273  // Automatic target selection
274  if (target == Target::AUTO) {
275  target = RegionDef->select_target(target);
276  }
277  // Check whether target exists and was generated for the given model and plugin
278  if (target == Target::UNKNOWN) {
279  fprintf(stderr, "Unkown target: %s\n", params.target_arg);
280  fprintf(stderr, "Available targets are: %s\n", get_target_list_string().c_str());
281  exit(1);
282  }
283  // Check if the main imp was generated for the requested target
284  if (RegionDef->select_target(target) == Target::UNKNOWN) {
285  fprintf(stderr, "Model %s was not generated for target %s\n",
286  params.imp_arg, get_string_from_target(target).c_str());
287  exit(1);
288  }
289  // Check the same thing if a plugin is specified
290  if (!RegionPlug.empty()) {
291  if (RegionPlug[0].get().select_target(target) == Target::UNKNOWN) {
292  fprintf(stderr, "Plugin %s was not generated for target %s\n",
293  params.plug_in_arg, get_string_from_target(target).c_str());
294  exit(1);
295  }
296  }
297 
298  // data structure to manage multiple IIF's
299  MULTI_IF MIIF;
300  MULTI_IF doppel, *cMIIF = &MIIF;
301 
302  MIIF.N_IIF = num_region;
303  MIIF.iontypes = {*RegionDef};
304  MIIF.targets = {target};
305  MIIF.IIFmask = Regions;
306  MIIF.numplugs = &num_plugins;
307  MIIF.plugtypes = {RegionPlug};
308  MIIF.numNode = Vmv->lsize();
309  MIIF.gdata[Vm] = Vmv;
310  MIIF.gdata[Iion] = I_ion;
311 
312  if (extVmUpdate)
313  MIIF.extUpdateVm = true;
314 
315  double setup_time = timing(t2, t1);
316  update_timing(timings+SETUP_IDX, setup_time);
317 
318  MIIF.initialize_MIIF();
319 
320  if (params.list_imps_flag || params.plugin_outputs_flag) {
321  print_models(params.plugin_outputs_flag);
322  exit(0);
323  }
324 
325  if (params.imp_info_flag) {
326  print_param_help(RegionDef, RegionPlug);
327  exit(0);
328  }
329 
330  fprintf(stderr, "\n*** GIT tag: %s\n", GIT_COMMIT_TAG);
331  fprintf(stderr, "*** GIT hash: %s\n", GIT_COMMIT_HASH);
332  fprintf(stderr, "*** GIT repo: %s\n", GIT_PATH);
333  fprintf(stderr, "*** dependency commits: %s\n\n", SUBREPO_COMMITS);
334 
335  if (params.buildinfo_flag) {
336  fprintf(stderr, "\n*** -buildinfo flag is deprecated and will be removed.\n");
337  fprintf( stderr, "\n*** The build information is printed by default. \n\n");
338  exit(0);
339  }
340 
341  MIIF.IIF[0]->tune(params.imp_par_arg, params.plug_in_arg, params.plug_par_arg);//FIXME
342 
343  // now compute the tables
344  MIIF.initialize_currents(dt, 1);
345 
346  if (analyze_AP) {
347  // initialize AP analyzer
348  double Vm0 = getCellVal(MIIF.gdata[Vm], 0);
349  for (int i = 0; i < VM_HIST_LEN; i++)
350  AP.vm_trc[i] = Vm0;
351 
352  char ap_stats_fname[1024];
353  snprintf(ap_stats_fname, sizeof ap_stats_fname, "%s_AP_stats.dat", params.fout_arg);
354  AP.stats = fopen(ap_stats_fname, "wt");
355  print_AP_stats_header(&AP, AP.stats);
356  }
357 
358  restitution r;
359  int n_dopple;
360  double *t_dopple;
361  if (restitute) {
362  char *res_file = dupstr(params.res_file_arg);
363  restitution_trigger_list(res_file, &r, params.restitute_arg, &n_dopple, &t_dopple);
364  if (r.dur > duration) duration = r.dur;
365 
366  // output restitution data
367  char ap_rstats_fname[1024];
368  snprintf(ap_rstats_fname, sizeof ap_rstats_fname, "%s_APD_restitution.dat", params.fout_arg);
369  AP.rstats = fopen(ap_rstats_fname, "wt");
370  print_AP_stats_header(&AP, AP.rstats);
371  }
372 
373  // read initial state vector from file and/or assignment
374  double t = 0.;
375  // if (params.restore_given)
376  // t = MIIF.restore_state( params.restore_arg, &pl, true);
377 
378  // illumination stuff
379  double light_irrad = params.light_irrad_arg;
380  double light_start = params.light_start_arg;
381  double light_end = params.duration_arg;
382  int light_nstim = params.light_numstim_arg;
383  double light_bcl = params.light_bcl_arg;
384  double light_dur = params.light_dur_arg;
385 
386  if( light_irrad < 0 ) {
387  FPRINTF( WORLD stderr, "WARNING: irradiance value %.3f makes no sense (clipped to 0.0)\n", light_irrad );
388  light_irrad = 0.0;
389  }
390  TrgList light_lst;
391  if(params.light_times_given) //user specified list
392  determine_stim_list(params.light_times_arg, &light_lst, false);
393 
394  if (params.SV_init_given)
395  initial_SVs(&MIIF, params.SV_init_arg, params.imp_arg, params.plug_in_arg, numNode);
396 
397  timer_manager tmo(dt, t, duration);
398  tmo.timers.resize(N_TIMERS);
399 
400  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0., CON_TM_IDX, "CON_TIMER");
401  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0., SVD_TM_IDX, "SVD_TIMER");
402 
403  if (params.save_time_arg)
404  tmo.initialize_singlestep_timer(params.save_time_arg, 0., STA_TM_IDX, "STA_TIMER");
405 
406  if (params.save_ini_time_given)
407  tmo.initialize_singlestep_timer(params.save_ini_time_arg, 0., SSV_TM_IDX, "SSV_TIMER");
408 
409  // use a trigger for stimulation
410  if (params.restitute_given) {
411  std::vector<double> trg;
412  trg.assign(r.trigs.lst, r.trigs.lst + r.trigs.n);
413  tmo.initialize_neq_timer(trg, stim_dur, STM_TM_IDX, "STIM_TIMER");
414 
415  if (params.res_state_vector_given) {
416  trg.assign(r.saveState.lst, r.saveState.lst + r.saveState.n);
417  tmo.initialize_neq_timer(trg, 0, RES_SAVE_TM_IDX, "RES_SAVE_TIMER");
418  }
419 
420  if (!strcmp(params.restitute_arg, "S1S2f") ) {
421  trg.assign(t_dopple, t_dopple + n_dopple);
422  tmo.initialize_neq_timer(trg, 1, DOPPLE_TM_IDX, "DOPPLE_TIMER");
423  }
424  }
425  else if (params.stim_times_given) {
426  std::vector<double> trg;
427  trg.assign(stim_lst.lst, stim_lst.lst + stim_lst.n);
428  tmo.initialize_neq_timer(trg, stim_dur, STM_TM_IDX, "STIM_TIMER");
429  }
430  else {
431  tmo.initialize_eq_timer(stim_start, stim_start+duration, nstim, bcl, stim_dur, STM_TM_IDX, "STIM_TIMER");
432  }
433 
434  // illumination timers
435  if ( MIIF.gdata[illum] != NULL && !tr_light ) {
436  if ( params.light_times_given ) {
437  std::vector<double> trg;
438  trg.assign(light_lst.lst, light_lst.lst + light_lst.n);
439  tmo.initialize_neq_timer(trg, light_dur, LIGHT_TM_IDX, "LIGHT_TIMER", NULL);
440  } else {
441  tmo.initialize_eq_timer(light_start, light_end, light_nstim,
442  light_bcl, light_dur, LIGHT_TM_IDX, "LIGHT_TIMER", NULL);
443  }
444  }
445 
446  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0, TRACE_TM_IDX, "TRACE_TIMER");
447 
448  if (IsVclamp && cl.transient)
449  tmo.initialize_singlestep_timer(params.clamp_start_arg, params.clamp_dur_arg, CLAMP_TM_IDX, "CLAMP_TIMER");
450  if (!restitute && params.doppel_on_given) {
451  tmo.initialize_singlestep_timer(params.doppel_on_arg, params.doppel_dur_arg, DOPPLE_TM_IDX, "DOPPLE_TIMER");
452  }
453 
454  // read initial state from single cell file and spread it out to all cells
455  if (params.read_ini_file_given) {
456  if (read_sv(&MIIF, R1, params.read_ini_file_arg) != 0) {
457  return 1; // exit with error code.
458  }
459  }
460 
461  // stretch experiment
462  stretch s;
463  memset(&s, 0, sizeof(stretch));
464  if (params.strain_given && MIIF.use_stretch()) {
465  initializePulseStretch(params.strain_arg, params.strain_time_arg,
466  params.strain_dur_arg, params.strain_rate_arg,
467  params.strain_rate_arg, &s);
468 
469  sf_vec* lambdavec; SF::init_vector(&lambdavec, MIIF.gdata[Vm]);
470  MIIF.gdata[Lambda] = lambdavec;
471  SF_real *l = lambdavec->ptr();
472 
473  for (int i = 0; i < lambdavec->lsize(); i++)
474  l[i] = s.pulse.sr;
475 
476  lambdavec->release_ptr(l);
477  }
478 
479  // set up assigning stimulus current to ion concentrations
480  std::stringstream stim_species(params.stim_species_arg);
481  std::stringstream stim_ratios(params.stim_ratios_arg);
482  std::string specie, ratio;
483  std::vector<std::string> species_list;
484  std::vector<float> ratios_list;
485  float ratio_sum = 0;
486  if (stim_assign) {
487  // split stim_species by ":" and then iterate over them
488  while(std::getline(stim_species, specie, ':'))
489  {
490  species_list.push_back(specie);
491  }
492  while(std::getline(stim_ratios, ratio, ':'))
493  {
494  ratios_list.push_back(std::stof(ratio));
495  }
496  if (species_list.size() != ratios_list.size()) {
497  fprintf(stderr, "Number of ion species for which stimulus should be assigned does not match number of ratio entries: %s | %s\n", params.stim_species_arg, params.stim_ratios_arg);
498  exit(1);
499  }
500  for (auto i : ratios_list )
501  ratio_sum += i;
502  if (std::abs(ratio_sum - 1) >= 0.001)
503  log_msg(NULL, 3, 0, "Warning: Sum of ion species ratios for stimulus assignement != 1. Normalizing to sum %.1f.", ratio_sum);
504  }
505 
506  // dumping of state variables goes here
507  char sv_dump_names[1024];
508  if (!params.fout_given && !params.validate_flag)
509  snprintf(sv_dump_names, sizeof sv_dump_names, "%s_%s", params.fout_arg, params.imp_arg);
510  else
511  strcpy(sv_dump_names, params.fout_arg);
512 
513  if (params.validate_flag) {
514  dump_all(&MIIF, R1, params.imp_arg, params.plug_in_arg, t, dt, sv_dump_names);
515  } else {
516  MIIF.sv_dump_add_by_name_list(R1, params.imp_arg, sv_dump_names,
517  params.imp_sv_dump_arg, params.plug_in_arg,
518  params.plug_sv_dump_arg, t, dt_out);
519  }
520 
521  // dumping of LUTs
522  if (params.dump_lut_flag)
523  MIIF.dump_luts_MIIF(true);
524 
525  // initialize tracing
526  int trace_nodes[] = {0};
527  if (do_trace)
528  open_trace(&MIIF, sizeof(trace_nodes)/sizeof(trace_nodes[0]), trace_nodes, &params.trace_no_arg, NULL);
529 
530  // dumping global vectors to file?
531  open_globalvec_dump(fhdls, &gvd, &MIIF, sv_dump_names, &io);
532 
533  // initialization time
534  double ini_time = timing(t3, t2);
535  update_timing(timings + INIT_IDX, ini_time);
536  doppel_MIIF(&MIIF, &doppel);
537 
538  int gpu_Vm_update = 1;
539  double counter = 0.0;
540  double epsilon = 0.0000001f;
541  /***********************************************/
542  /***************** MAIN LOOP *******************/
543  /***********************************************/
544 
545  // Print simulation target
546  fprintf(stderr, "Running simulation on target %s\n", get_string_from_target(target).c_str());
547  while (!tmo.elapsed())
548  {
549  get_time(t4);
550  t = tmo.time;
551 
552  if (tmo.triggered_now(DOPPLE_TM_IDX) ) {
553  printf("[%f] Switching to new ion structure\n", tmo.time);
554  doppel_update(&MIIF, &doppel);
555  cMIIF = &doppel;
556  }
557 
558  if (IsVclamp && cl.transient && tmo.trigger(CLAMP_TM_IDX) )
559  cMIIF->extUpdateVm = true;
560  else
561  cMIIF->extUpdateVm = extVmUpdate;
562 
563  if (IsVclamp)
564  clamp_signal(cMIIF, &cl, &tmo);
565 
566  if (APclamp)
567  AP_clamp(&ap_cl, &tmo, cMIIF->gdata[Vm], tmo.triggered_now(STM_TM_IDX) );
568 
569  for (int i = 0; i < SVclamp; i++) {
570  bool trigger;
571  if (params.SV_I_trigger_flag)
572  trigger = tmo.triggered_now(STM_TM_IDX);
573  else
574  trigger = !tmo.d_time;
575 
576  sv_clamp(sv_cl+i, &tmo, cMIIF, trigger);
577  }
578 
579  // output selected state variables
580  cMIIF->dump_svs(tmo.timers[SVD_TM_IDX]);
581  globalvec_dump(fhdls, &gvd, cMIIF, &tmo, &io, numNode);
582 
583  // save state of all cells to binary file
584  if (tmo.trigger(STA_TM_IDX) )
585  cMIIF->dump_state(params.save_file_arg, tmo.time, mesh_t::unset_msh, false, 0);
586 
587  // save state of cell 0 to text file
588  if (tmo.trigger(SSV_TM_IDX) )
589  save_sv(cMIIF, R1, params.save_ini_file_arg);
590 
591  // save state within an restitution protocol
592  if (tmo.trigger(RES_SAVE_TM_IDX) )
593  restitution_save_sv(cMIIF, R1, &r, &AP);
594 
595  // add stim current
596  if (tmo.trigger(STM_TM_IDX) ) {
597  gpu_Vm_update = 0;
598  if (tr_stim) {
599  i_stim = stim_trace.s[tmo.trigger_elapse(STM_TM_IDX)]*dt;
600  }
601  else if (params.stim_volt_given) {
602  SF_real *d = cMIIF->gdata[Vm]->ptr();
603  // peak_vm = last_vm = d[0];
604  i_stim = (params.stim_volt_arg-d[0])/params.resistance_arg;
605  cMIIF->gdata[Vm]->release_ptr(d);
606  }
607 
608 #ifndef ALGEBRAIC
609  *(cMIIF->gdata[Vm]) += i_stim;
610 #endif // ifndef ALGEBRAIC
611 
612  if (stim_assign) {
613  int iIon = 0;
614  for (auto ion : species_list)
615  {
616  cMIIF->transmem_stim_species( -i_stim*(ratios_list[iIon]/ratio_sum), ion.c_str(), params.surface_to_volume_arg, stim_list, numNode);
617  FPRINTF( WORLD stderr, "[stim_assign] @ t = %.2f ms; magnitude: %fpA/pF; type: %s\n", tmo.time, i_stim, ion.c_str() );
618  iIon++;
619  }
620  }
621  }
622 
623  if ( MIIF.gdata[illum] != NULL ) {
624  MIIF.gdata[illum]->set(tmo.trigger(LIGHT_TM_IDX) ? light_irrad : 0.0);
625 
626  if(tmo.triggered_now(LIGHT_TM_IDX) && light_irrad) {
627  FPRINTF( WORLD stderr, "[stim] illum ON @ t = %.2f ms; E_e = %.3f mW/mm^2\n", tmo.time, light_irrad );
628  }
629  if(tmo.trigger_end(LIGHT_TM_IDX) && light_irrad) {
630  FPRINTF( WORLD stderr, "[stim] illum OFF @ t = %.2f ms\n", tmo.time );
631  }
632  }
633 
634  // strain
635  if (params.strain_given && cMIIF->use_stretch())
636  apply_stretch(cMIIF, &s, &tmo);
637 
638  if (analyze_AP) {
639  if (restitute) AP.pmat = r.trigs.pmat[AP.beat];
640  check_events(getCellVal(cMIIF->gdata[Vm], 0), &AP, &tmo);
641  }
642 
643  get_time(t5);
644 
645 
646  if (is_gpu(target)) {
647  if (gpu_Vm_update == 0) {
648  cMIIF->compute_ionic_current(1, 1);
649  } else if (gpu_Vm_update == 1) {
650  cMIIF->compute_ionic_current(fabs(counter - 0.0) < epsilon || fabs(counter - 1.0) < epsilon , fabs((counter+dt) - 1.0) < epsilon);
651  }
652  }
653  else {
654  cMIIF->compute_ionic_current();
655  }
656 
657  // update Vm from Iion
658  if (extVmUpdate)
659  cMIIF->gdata[Vm]->add_scaled(*cMIIF->gdata[Iion], -cMIIF->dt);
660 
661  double ode_time = timing(t6, t5);
662  update_timing(timings+ODE_IDX, ode_time);
663 
664  if (do_trace && tmo.trigger(TRACE_TM_IDX)) {
665  dump_trace(cMIIF, t);
666  }
667 
668  if (!restitute && tmo.trigger_end(DOPPLE_TM_IDX) ) {
669  log_msg(NULL, 0, 0, "[%f] Switching back to original ion structure\n", tmo.time);
670  cMIIF = &MIIF;
671  }
672 
673  tmo.update_timers();
674 
675  double loop_time = timing(t7, t4);
676  update_timing(timings+LOOP_IDX, loop_time);
677 
678  gpu_Vm_update = 1;
679  counter = counter + dt;
680  if ( fabs(counter - 1.0) < epsilon ) {
681  counter = 0.0;
682  }
683  }//loop end
684 
685  // loop_current time
686  double main_time = timing(t8, t3);
687 
688  // write header so that we know which files belong together
689  const char *ExpID = params.fout_arg;
690  write_dump_header(&gvd, &MIIF.svd, ExpID);
691 
692  // clean up everything
693  if (tr_stim) free_trace(&stim_trace);
694  if (analyze_AP) cleanup_AP_analysis(&AP);
695  if (do_trace) close_trace(&MIIF);
696  MIIF.close_svs_dumps();
697  free_doppel(&doppel);
698  MIIF.free_MIIF();
699 
700  close_globalvec_dump(fhdls, &gvd, &io);
701 
702  log_msg(NULL, 0, 0, "\n\n\nAll done!\n\n");
703 
704  // print timing data
705  event_timing *t_setup = timings+SETUP_IDX;
706  log_msg(NULL, 0, 0, "setup time %.6f s", (double)t_setup->tot);
707  event_timing *t_init = timings+INIT_IDX;
708  log_msg(NULL, 0, 0, "initialization time %.6f s", (double)t_init->tot);
709  event_timing *t_loop = timings+LOOP_IDX;
710  log_msg(NULL, 0, 0, "main loop time %.6f s", (double)t_loop->tot);
711  event_timing *t_ode = timings+ODE_IDX;
712  log_msg(NULL, 0, 0, "total ode time %.6f s\n", (double)t_ode->tot);
713 
714  log_msg(NULL, 0, 0, "mn/avg/mx loop time %.6f %.6f %.6f ms",
715  (double)t_loop->mn*1e3, (double)t_loop->avg*1e3, (double)t_loop->mx*1e3);
716  log_msg(NULL, 0, 0, "mn/avg/mx ODE time %.6f %.6f %.6f ms",
717  (double)t_ode->mn*1e3, (double)t_ode->avg*1e3, (double)t_ode->mx*1e3);
718 
719  log_msg(NULL, 0, 0, "real time factor %.6f\n\n", duration/(t_ode->tot*1e3) );
720 
721 
722  CHKERRQ(PetscFinalize());
723  return 0;
724 } // main
int trigger_elapse(int ID) const
Definition: timer_utils.h:186
timing for setup phase
Definition: bench_utils.h:75
int main(int argc, char *argv[])
Definition: bench.cc:86
int read_sv(MULTI_IF *, int, const char *)
void update_timing(event_timing *t, double event_duration)
Definition: bench_utils.cc:376
void free_trace(trace *tr)
Definition: trace.cc:76
bool pmat
indiacte premature AP
Definition: ap_analyzer.h:91
void restitution_trigger_list(char *r_file, restitution *r, char *protocol, int *n_dop, double **t_dop)
Definition: restitute.cc:45
std::string get_string_from_target(Target const target)
Get a string representation of a given target.
Definition: target.cc:46
TrgList trigs
trigger list for defining stim sequence
Definition: restitute.h:50
virtual void release_ptr(S *&p)=0
void open_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, char *base_name, IOCtrl *io)
Definition: bench_utils.cc:176
int transient
Definition: clamp.h:37
char * dupstr(const char *old_str)
Definition: basics.cc:44
int beat
beat counter
Definition: ap_analyzer.h:89
virtual void add_scaled(const abstract_vector< T, S > &vec, S k)=0
void close_trace(MULTI_IF *MIIF)
void dump_all(MULTI_IF *MIIF, int reg, char *imp, char *plugs, double t, double ddt, char *fout)
Definition: bench_utils.cc:344
double getCellVal(sf_vec *v, int ind)
Definition: bench_utils.cc:392
int read_trace(trace *tr, const char *name)
Definition: trace.cc:43
SV_DUMP svd
state variable dump
Definition: MULTI_ION_IF.h:203
void resample_trace(trace *tr, double dt)
Definition: trace.cc:115
std::vector< std::reference_wrapper< IonType > > IonTypeList
Definition: ion_type.h:289
void sv_dump_add_by_name_list(int, char *, char *, char *, char *, char *, double, double)
void initialize_timings(event_timing *t)
Definition: bench_utils.cc:366
FILE * rstats
output restitution statistics only
Definition: ap_analyzer.h:104
int * numplugs
number of plugins for each region
Definition: MULTI_ION_IF.h:209
float determine_duration(struct gengetopt_args_info *p, TrgList *stim_lst)
determine time of last stimulus
Definition: bench_utils.cc:503
double avg
average duration of event
Definition: bench_utils.h:69
double dur
duration
Definition: trace.h:29
pulseStretch pulse
Definition: stretch.h:52
void sv_clamp(Clamp *cl, timer_manager *tm, MULTI_IF *miif, bool trigger)
Definition: clamp.cc:138
void compute_ionic_current(bool flag_send=1, bool flag_receive=1)
GPU kernel to emulate the add_scaled call made to adjust the Vm values when the update to Vm is not m...
timing for initialization
Definition: bench_utils.h:76
void doppel_MIIF(MULTI_IF *orig, MULTI_IF *miif_doppel)
bool initialize_clamp(Clamp *cl, double cl_val, double ini_val, double start, double dur, const char *f, int trans, float *duration)
Definition: clamp.cc:45
std::map< mesh_t, sf_mesh > mesh_reg
Registry for the different meshes used in a multi-physics simulation.
Definition: main.cc:46
IIF_Mask_t * IIFmask
region for each node
Definition: MULTI_ION_IF.h:214
void initialize_eq_timer(double istart, double iend, int ntrig, double iintv, double idur, int ID, const char *iname, const char *poolname=nullptr)
Definition: timer_utils.cc:48
void dump_luts_MIIF(bool)
void print_AP_stats_header(action_potential *AP, FILE *outbuf)
Definition: ap_analyzer.cc:526
virtual T lsize() const =0
timer_manager * tm_manager
a manager for the various physics timers
Definition: main.cc:52
bool trigger_end(int ID) const
Definition: timer_utils.h:178
void save_sv(MULTI_IF *, int, const char *)
void print_param_help(IonType *im, IonTypeList &plugs)
Definition: bench_utils.cc:468
timing for main loop (including IO and ODE solve)
Definition: bench_utils.h:77
IonTypeList iontypes
type for each region
Definition: MULTI_ION_IF.h:212
char wbin
write to file in binary format
Definition: bench_utils.h:50
bool check_events(double vm, action_potential *AP, timer_manager *tm)
Definition: ap_analyzer.cc:263
int get_plug_flag(char *plgstr, int *out_num_plugins, IonTypeList &out_plugins)
std::vector< IonTypeList > plugtypes
plugins types for each region
Definition: MULTI_ION_IF.h:215
bool trigger(int ID) const
Definition: timer_utils.h:166
int N_IIF
how many different IIF&#39;s
Definition: MULTI_ION_IF.h:211
void initialize_currents(double, int)
bool is_gpu(Target const target)
Checks if this is a GPU target.
Definition: target.cc:64
TrgList saveState
instants at wich we save state vectors
Definition: restitute.h:51
void init_vector(SF::abstract_vector< T, S > **vec)
Definition: SF_init.h:99
double mx
maximum duration of event
Definition: bench_utils.h:68
std::vector< IonIfBase * > IIF
array of IIF&#39;s
Definition: MULTI_ION_IF.h:202
int n
number of pulses required for protocol
Definition: restitute.h:29
void apply_stretch(MULTI_IF *miif, stretch *s, timer_manager *tm)
Definition: stretch.cc:55
void print_models(bool)
virtual Target select_target(Target target) const =0
Gets a supported target from the given target.
int initialize_AP_analysis(action_potential *AP)
Definition: ap_analyzer.cc:71
special value to handle unknown targets
Definition: target.h:47
double tot
total duration of all events
Definition: bench_utils.h:70
void doppel_update(MULTI_IF *orig, MULTI_IF *miif_doppel)
virtual S * ptr()=0
void dump_state(char *, float, opencarp::mesh_t gid, bool, unsigned int)
void initialize_singlestep_timer(double tg, double idur, int ID, const char *iname, const char *poolname=nullptr)
Definition: timer_utils.h:156
Define multiple ionic models to be used in different regions.
void initializePulseStretch(float strain, float onset, float duration, float rise, float fall, stretch *s)
Definition: stretch.cc:38
double time
current time
Definition: timer_utils.h:76
IonType * get_ion_type(const std::string &name)
std::map< SF::quadruple< int >, SF::index_mapping< int > > map_reg
Registriy for the inter domain mappings.
Definition: main.cc:48
double * lst
store instants of pulse delivery
Definition: restitute.h:30
void cleanup_AP_analysis(action_potential *AP)
Definition: ap_analyzer.cc:169
virtual void set(const vector< T > &idx, const vector< S > &vals, const bool additive=false)=0
double vm_trc[VM_HIST_LEN]
Definition: ap_analyzer.h:100
void free_doppel(MULTI_IF *m)
float * s
samples
Definition: trace.h:28
void determine_stim_list(char *stl, TrgList *trg, bool DIAs)
Definition: bench_utils.cc:101
long d_time
current time instance index
Definition: timer_utils.h:77
std::string get_target_list_string()
Returns a string containing the list of available targets.
Definition: target.cc:55
void clamp_signal(MULTI_IF *pMIIF, Clamp *cl, timer_manager *tm)
Definition: clamp.cc:69
void restitution_save_sv(MULTI_IF *miif, int R1, restitution *r, action_potential *AP)
Definition: restitute.cc:275
#define Vm_clamp
Definition: stimulate.h:47
char w2file
write to file
Definition: bench_utils.h:49
void globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, timer_manager *tmo, IOCtrl *io, int numNode)
Definition: bench_utils.cc:272
void get_time(double &tm)
Definition: basics.h:436
FILE * stats
output statistics for each AP
Definition: ap_analyzer.h:103
void open_trace(MULTI_IF *MIIF, int n_traceNodes, int *traceNodes, int *label, opencarp::sf_mesh *imesh)
Set up ionic model traces at some global node numbers.
number of benchmark timings we use
Definition: bench_utils.h:79
V timing(V &t2, const V &t1)
Definition: basics.h:448
Index mapping class. This is a bijective mapping.
Definition: SF_container.h:207
int write_dump_header(GVEC_DUMP *gvd, SV_DUMP *svd, const char *ExpID)
Definition: bench_utils.cc:135
V dot(const vec3< V > &p1, const vec3< V > &p2)
Definition: vect.h:125
int dump_svs(opencarp::base_timer *)
void transmem_stim_species(float, const char *, float, int *, int)
opencarp::sf_vec * gdata[NUM_IMP_DATA_TYPES]
data used by all IMPs
Definition: MULTI_ION_IF.h:216
Defines valid targets for an ionic model to run on and an allocator for allocating memory on a specif...
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
Definition: basics.cc:72
#define VM_HIST_LEN
Definition: ap_analyzer.h:87
timing for ODE solve
Definition: bench_utils.h:78
Target get_target_from_string(std::string const str)
Returns a value from the Target enum from a given string.
Definition: target.cc:36
int N
number of samples
Definition: trace.h:26
void dump_trace(MULTI_IF *MIIF, limpet::Real time)
std::vector< Target > targets
target for each region
Definition: MULTI_ION_IF.h:213
int numNode
local number of nodes
Definition: MULTI_ION_IF.h:210
Target
enum that represents different targets to run ionic models on.
Definition: target.h:45
std::vector< base_timer * > timers
vector containing individual timers
Definition: timer_utils.h:84
double dur
total duration of protocol
Definition: restitute.h:52
bool triggered_now(int ID) const
Definition: timer_utils.h:171
void AP_clamp(Clamp *cl, timer_manager *tm, sf_vec *v, bool trigger)
Definition: clamp.cc:178
void initialize_neq_timer(const std::vector< double > &itrig, double idur, int ID, const char *iname, const char *poolname=nullptr)
Definition: timer_utils.cc:63
int load_ionic_module(const char *)
double dt
time step (ms)
Definition: MULTI_ION_IF.h:219
void close_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, IOCtrl *io)
Definition: bench_utils.cc:315
double mn
minimum duration of event
Definition: bench_utils.h:67
Abstract class representing an ionic model type.
Definition: ion_type.h:59
double SF_real
Use the general double as real type.
Definition: SF_globals.h:38
void initial_SVs(MULTI_IF *miif, char *SVs, char *imp, char *plgins, int num)
Definition: bench_utils.cc:412
void initialize_sv_clamp(Clamp *cl, const char *sv, char *file, double dt)
Definition: clamp.cc:117
bool extUpdateVm
flag indicating update function for Vm
Definition: MULTI_ION_IF.h:208
Basic utility structs and functions, mostly IO related.
int process_sv_clamps(char *SVs, char *files, Clamp **clamps, double dt)
Definition: clamp.cc:197
bool * pmat
store flag to indicate prematurity
Definition: restitute.h:31
int get_size(MPI_Comm comm=PETSC_COMM_WORLD)
Definition: basics.h:290
manage input, output, resampling of traces
Definition: trace.h:25
centralize time managment and output triggering
Definition: timer_utils.h:73
SF::scatter_registry scatter_reg
Registry for the different scatter objects.
Definition: main.cc:44
The scatterer registry class.