openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
electrics.h
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 #ifndef _ELECTRICS_H
28 #define _ELECTRICS_H
29 
30 #include "physics_types.h"
31 #include "sim_utils.h"
32 #include "stimulate.h"
33 
34 #include "sf_interface.h"
35 #include "timers.h"
36 #include "ionics.h"
37 
38 #define DUMP_NONE 0
39 #define DUMP_IC 1
40 #define DUMP_IVOL 2
41 #define DUMP_IACT 4
42 
43 namespace opencarp {
44 
46 {
47  public:
48  sf_vec* phie = nullptr;
49  sf_vec* phie_i = nullptr;
50  sf_vec* phiesrc = nullptr;
51  sf_vec* currtmp = nullptr;
52 
53  sf_mat* mass_e = nullptr;
54  sf_mat* phie_mat = nullptr;
55 
56  // the linear solver
57  sf_sol* lin_solver = nullptr;
58 
59  // linear solver stats
61 
63  dbc_manager* dbc = nullptr;
64  bool phie_mat_has_nullspace = false;
65 
66  // solver config
67  double tol = 1e-8;
68  int max_it = 100;
69 
70  void init();
71  void rebuild_matrices(MaterialType* mtype, SF::vector<stimulus> & stimuli, FILE_SPEC logger);
72  void rebuild_stiffness(MaterialType* mtype, SF::vector<stimulus> & stimuli, FILE_SPEC logger);
73  void rebuild_mass(FILE_SPEC logger);
74  void solve(sf_mat & Ki, sf_vec & Vmv, sf_vec & tmp_i);
75  void solve_laplace();
76 
78  {
79  if(dbc) delete dbc;
80  if(lin_solver) delete lin_solver;
81  // matrices
82  if(phie_mat) delete phie_mat;
83  if(mass_e) delete mass_e;
84  // vectors
85  if(phie) delete phie;
86  if(phie_i) delete phie_i;
87  if(phiesrc) delete phiesrc;
88  if(currtmp) delete currtmp;
89  }
90 
91  private:
92  void enforce_dbc();
93  void setup_linear_solver(FILE_SPEC logger);
94 };
95 
97 {
98  public:
99 
100  // Has to be kept in line with param_globals::parab_solve
101  enum parabolic_t { EXPLICIT = 0, CN = 1, O2dT = 2};
102 
103  sf_vec* IIon = nullptr;
104  sf_vec* Vmv = nullptr;
105 
106  sf_vec* Ic = nullptr;
107  sf_vec* Ivol = nullptr;
108  sf_vec* Iact = nullptr;
109 
110  sf_vec* old_vm = nullptr;
111  sf_vec* kappa_i = nullptr;
112  sf_vec* tmp_i1 = nullptr;
113  sf_vec* tmp_i2 = nullptr;
114  sf_vec* Irhs = nullptr;
115  sf_vec* inv_mass_diag = nullptr;
116 
117  sf_mat* u_mass_i = nullptr;
118  sf_mat* mass_i = nullptr;
119  sf_mat* rhs_parab = nullptr;
120  sf_mat* lhs_parab = nullptr;
121  sf_mat* phie_recov_mat = nullptr;
122 
123  sf_vec* Diff_term = nullptr;
124 
125  // the linear solver
126  sf_sol* lin_solver = nullptr;
127 
128  // linear solver stats
130 
131  // solver config
132  double tol = 1e-8;
133  int max_it = 100;
134  parabolic_t parab_tech = CN;
135 
136  // solver statistics
137  double final_residual = -1.0;
138  int niter = -1;
139 
141  {
142  if (lin_solver) delete lin_solver;
143  // matrices
144  if (u_mass_i) delete u_mass_i;
145  if (mass_i) delete mass_i;
146  if (rhs_parab) delete rhs_parab;
147  if (lhs_parab) delete lhs_parab;
148  if (phie_recov_mat) delete phie_recov_mat;
149  // vectors
150  // IIon and Vmv are shallow copies of global vectors, do not delete them
151  if (old_vm) delete old_vm;
152  if (kappa_i) delete kappa_i;
153  if (tmp_i1) delete tmp_i1;
154  if (tmp_i2) delete tmp_i2;
155  if (Irhs) delete Irhs;
156  if (inv_mass_diag) delete inv_mass_diag;
157  }
158 
159  void init();
160  void rebuild_matrices(MaterialType* mtype, limpet::MULTI_IF & miif, FILE_SPEC logger);
161  void solve(sf_vec & phie_i);
162 
163  private:
164  void setup_linear_solver(FILE_SPEC logger);
165 
166  void solve_CN(sf_vec & phie_i);
167  void solve_O2dT(sf_vec & phie_i);
168  void solve_EF(sf_vec & phie_i);
169 };
170 
171 enum PotType {VM, PHIE};
173 
175 struct Activation {
177  float threshold;
178  int mode;
179  int all;
180  int init;
181  sf_vec* phi = nullptr;
182  sf_vec* phip = nullptr;
183  int *ibuf = nullptr;
184  double *actbuf = nullptr;
185  sf_vec* tm = nullptr;
186  sf_vec* dvp0 = nullptr;
187  sf_vec* dvp1 = nullptr;
189  char *fname = nullptr;
190  char *ID = nullptr;
191  char *prv_fname = nullptr;
193  int offset;
194  int nacts;
195 };
196 
198 struct Sentinel {
199  bool activated = false;
200  double t_start = -1.0;
201  double t_window = -1.0;
202  double t_quiesc = -1.0;
203  int ID = -1;
204 };
205 
207 {
208  private:
210  int check_cross_threshold(sf_vec & vm, sf_vec & vmp, double tm,
211  int *ibuf, double *actbuf, float threshold, int mode);
212 
214  int check_mx_derivative(sf_vec & vm, sf_vec & vmp, double tm,
215  int *ibuf, double *actbuf, sf_vec & dvp0, sf_vec & dvp1,
216  float threshold, int mode);
217 
218  public:
222 
224  LAT_detector();
225 
227  void init(sf_vec & vm, sf_vec & phie, int offset, enum physic_t = elec_phys);
228 
230  int check_acts(double tm);
231 
233  int check_quiescence(double tm, double dt);
234 
236  void output_initial_activations();
237 };
238 
240 {
242  sf_vec* phie_rec = nullptr;
243  sf_vec* Im = nullptr;
244  sf_vec* dphi = nullptr;
246 };
247 
248 class Electrics : public Basic_physic
249 {
250  public:
257  enum grid_t {intra_grid = 0, extra_grid};
258  enum linsys_t {elliptic_sys, parabolic_sys};
259 
261  MaterialType mtype[2];
264 
268 
273 
276 
279 
282 
285 
287 
292  {
293  name = "Electrics";
294  }
295 
304  void initialize();
305 
306  void destroy();
307 
308  // This funcs from the Basic_physic interface are currently empty
309  void compute_step();
310  void output_step();
311 
312  inline void output_timings()
313  {
314  // since ionics are included in electrics, we substract ionic timings from electric timings
315  compute_time -= ion.compute_time;
316  initialize_time -= ion.initialize_time;
317 
318  // now we call the base class timings output for Electrics
320  // and then for ionics
321  ion.output_timings();
322  }
323 
324  ~Electrics();
325 
327  double timer_val(const int timer_id);
328 
330  std::string timer_unit(const int timer_id);
331 
332 
333  private:
334 
336  void setup_stimuli();
337 
339  void stimulate_intracellular();
342  void clamp_Vm();
343 
345  void stimulate_extracellular();
346 
348  void setup_solvers();
349 
351  void setup_mappings();
352 
354  void setup_output();
355 
357  void dump_matrices();
358 
359  void checkpointing();
360 
361  void balance_electrodes();
362 
363  void prepace();
364 };
365 
373 
375 bool have_dbc_stims(const SF::vector<stimulus> & stimuli);
376 
378 int stimidx_from_timeridx(const SF::vector<stimulus> & stimuli, const int timer_id);
379 
381 void recover_phie_std(sf_vec & vm, phie_recovery_data & rcv);
383 
384 class Laplace : public Basic_physic
385 {
386  public:
388  MaterialType mtype[2];
391 
396 
398  {
399  name = "Laplace solver";
400  }
401 
402  void initialize();
403  void destroy();
404  void compute_step();
405  void output_step();
407  double timer_val(const int timer_id);
409  std::string timer_unit(const int timer_id);
410 };
411 
423  sf_mat & mass_i,
424  sf_mat & mass_e,
425  limpet::MULTI_IF *miif,
426  FILE_SPEC logger);
427 
428 void balance_electrode(SF::vector<stimulus> & stimuli, int balance_from, int balance_to);
429 
430 void apply_stim_to_vector(const stimulus & s, sf_vec & vec, bool add);
431 
432 void set_cond_type(MaterialType & m, cond_t type);
433 
436 const char* get_tsav_ext(double time);
437 
440 void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid,
441  SF::vector<mesh_int_t>* & restr, bool async = false);
442 
443 void compute_restr_idx_async(sf_mesh & mesh,
444  SF::vector<mesh_int_t> & inp_idx,
445  SF::vector<mesh_int_t> & idx);
446 
447 void compute_restr_idx(sf_mesh & mesh,
448  SF::vector<mesh_int_t> & inp_idx,
449  SF::vector<mesh_int_t> & idx);
450 
451 } // namespace opencarp
452 
453 #endif
int nacts
number of events detected over time step
Definition: electrics.h:194
SF_real gBath
Bath conductivity.
Definition: electrics.h:245
bool have_dbc_stims(const SF::vector< stimulus > &stimuli)
return wheter any stimuli require dirichlet boundary conditions
Definition: electrics.cc:856
dbc_manager * dbc
dbcs require a dbc manager
Definition: electrics.h:63
igb_output_manager output_manager
class handling the igb output
Definition: electrics.h:281
gvec_data gvec
datastruct holding global IMP state variable output
Definition: electrics.h:278
description of materal properties in a mesh
Definition: fem_types.h:110
void setup_phie_recovery_data(phie_recovery_data &data)
Definition: electrics.cc:2052
void recover_phie_std(sf_vec &vm, phie_recovery_data &rcv)
Definition: electrics.cc:1903
SF::vector< mesh_real_t > pts
The phie recovery locations.
Definition: electrics.h:241
sf_vec * currtmp
temp vector for phiesrc
Definition: electrics.h:51
cond_t
description of electrical tissue properties
Definition: fem_types.h:42
void solve(sf_mat &Ki, sf_vec &Vmv, sf_vec &tmp_i)
Definition: electrics.cc:1048
SF::vector< stimulus > stimuli
the electrical stimuli
Definition: electrics.h:390
elliptic_solver ellip_solver
Solver for the elliptic bidomain equation.
Definition: electrics.h:270
lin_solver_stats stats
Definition: electrics.h:60
SF::index_mapping< mesh_int_t > petsc_to_nodal
Definition: electrics.h:220
int all
determine all or first instants of activation
Definition: electrics.h:179
SF::vector< Activation > acts
Definition: electrics.h:219
void constant_total_stimulus_current(SF::vector< stimulus > &stimuli, sf_mat &mass_i, sf_mat &mass_e, limpet::MULTI_IF *miif, FILE_SPEC logger)
Definition: electrics.cc:2216
void balance_electrode(elliptic_solver &ellip, SF::vector< stimulus > &stimuli, int balance_from, int balance_to)
Definition: electrics.cc:393
mesh_t
The enum identifying the different meshes we might want to load.
Definition: sf_interface.h:58
physic_t
Identifier for the different physics we want to set up.
Definition: physics_types.h:51
PotType measurand
quantity being monitored
Definition: electrics.h:192
generic_timing_stats IO_stats
Definition: electrics.h:286
void rebuild_mass(FILE_SPEC logger)
Definition: electrics.cc:1011
int mode
toggle mode from standard to reverse
Definition: electrics.h:178
sf_mat * phie_mat
lhs matrix to solve elliptic
Definition: electrics.h:54
igb_output_manager output_manager
class handling the igb output
Definition: electrics.h:395
parabolic_solver parab_solver
Solver for the parabolic bidomain equation.
Definition: electrics.h:272
Electrics()
Most of the initialization is done with initialize()
Definition: electrics.h:291
event detection data structures
Definition: electrics.h:175
for analysis of the #iterations to solve CG
Definition: timers.h:69
The abstract physics interface we can use to trigger all physics.
Definition: physics_types.h:59
int init
true if intialized
Definition: electrics.h:180
ActMethod method
method to check whether activation occured
Definition: electrics.h:176
sf_sol * lin_solver
petsc or ginkgo lin_solver
Definition: electrics.h:57
float threshold
threshold for detection of activation
Definition: electrics.h:177
sentinel for checking activity in the tissue
Definition: electrics.h:198
const char * get_tsav_ext(double time)
Definition: electrics.cc:865
int max_it
maximum number of iterations
Definition: electrics.h:68
manager for dirichlet boundary conditions
Definition: stimulate.h:192
Electrical ionics functions and LIMPET wrappers.
sf_vec * phie
phi_e
Definition: electrics.h:48
void set_elec_tissue_properties(MaterialType *mtype, Electrics::grid_t g, FILE_SPEC logger)
Fill the RegionSpec of an electrics grid with the associated inputs from the param parameters...
Definition: electrics.cc:100
void set_cond_type(MaterialType &m, cond_t type)
Definition: electrics.cc:833
elliptic_solver ellip_solver
Solver for the elliptic bidomain equation.
Definition: electrics.h:393
Interface to SlimFem.
int offset
node number offset (used for PS)
Definition: electrics.h:193
void compute_restr_idx(sf_mesh &mesh, SF::vector< mesh_int_t > &inp_idx, SF::vector< mesh_int_t > &idx)
Definition: electrics.cc:544
SF::vector< stimulus > stimuli
the electrical stimuli
Definition: electrics.h:263
File descriptor struct.
Definition: basics.h:133
virtual void output_timings()
Definition: physics_types.h:76
sf_mat * mass_e
mass matrix for RHS elliptic calc
Definition: electrics.h:53
double tol
CG stopping tolerance.
Definition: electrics.h:67
sf_vec * phiesrc
I_e.
Definition: electrics.h:50
void rebuild_stiffness(MaterialType *mtype, SF::vector< stimulus > &stimuli, FILE_SPEC logger)
Definition: electrics.cc:951
Electrical stimulation functions.
lin_solver_stats stats
Definition: electrics.h:129
A vector storing arbitrary data.
Definition: SF_vector.h:42
grid_t
An electrics grid identifier to distinguish between intra and extra grids.
Definition: electrics.h:257
sf_vec * phie_i
phi_e on intracellular grid
Definition: electrics.h:49
for analysis of the #iterations to solve CG
Definition: timers.h:13
LAT_detector lat
the activation time detector
Definition: electrics.h:275
FILE_SPEC fout
output file
Definition: electrics.h:188
void rebuild_matrices(MaterialType *mtype, SF::vector< stimulus > &stimuli, FILE_SPEC logger)
Definition: electrics.cc:940
double SF_real
Use the general double as real type.
Definition: SF_globals.h:38
void apply_stim_to_vector(const stimulus &s, sf_vec &vec, bool add)
Definition: electrics.cc:452
void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid, SF::vector< mesh_int_t > *&restr, bool async)
Definition: electrics.cc:612
int stimidx_from_timeridx(const SF::vector< stimulus > &stimuli, const int timer_id)
determine link between timer and stimulus
Definition: electrics.cc:781
Simulator-level utility execution control functions.
phie_recovery_data phie_rcv
struct holding helper data for phie recovery
Definition: electrics.h:284
void compute_restr_idx_async(sf_mesh &mesh, SF::vector< mesh_int_t > &inp_idx, SF::vector< mesh_int_t > &idx)
Definition: electrics.cc:577
int postproc_recover_phie()
Definition: electrics.cc:1969
Basic physics types.