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;
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 
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:
259 
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
317 
318  // now we call the base class timings output for Electrics
320  // and then for ionics
322  }
323 
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:
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 
438  sf_mat & mass_i,
439  sf_mat & mass_e,
440  limpet::MULTI_IF *miif,
441  FILE_SPEC logger);
442 
459 void balance_electrode(SF::vector<stimulus> & stimuli, int balance_from, int balance_to);
460 
461 void apply_stim_to_vector(const stimulus & s, sf_vec & vec, bool add);
462 
463 void set_cond_type(MaterialType & m, cond_t type);
464 
467 const char* get_tsav_ext(double time);
468 
471 void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid,
472  SF::vector<mesh_int_t>* & restr, bool async = false);
473 
474 void compute_restr_idx_async(sf_mesh & mesh,
475  SF::vector<mesh_int_t> & inp_idx,
476  SF::vector<mesh_int_t> & idx);
477 
478 void compute_restr_idx(sf_mesh & mesh,
479  SF::vector<mesh_int_t> & inp_idx,
480  SF::vector<mesh_int_t> & idx);
481 
482 } // namespace opencarp
483 
484 #endif
double SF_real
Use the general double as real type.
Definition: SF_globals.h:38
A vector storing arbitrary data.
Definition: SF_vector.h:43
The abstract physics interface we can use to trigger all physics.
Definition: physics_types.h:59
virtual void output_timings()
Definition: physics_types.h:76
const char * name
The name of the physic, each physic should have one.
Definition: physics_types.h:62
SF::vector< stimulus > stimuli
the electrical stimuli
Definition: electrics.h:263
LAT_detector lat
the activation time detector
Definition: electrics.h:275
grid_t
An electrics grid identifier to distinguish between intra and extra grids.
Definition: electrics.h:257
phie_recovery_data phie_rcv
struct holding helper data for phie recovery
Definition: electrics.h:284
generic_timing_stats IO_stats
Definition: electrics.h:286
void destroy()
Currently we only need to close the file logger.
Definition: electrics.cc:378
gvec_data gvec
datastruct holding global IMP state variable output
Definition: electrics.h:278
elliptic_solver ellip_solver
Solver for the elliptic bidomain equation.
Definition: electrics.h:270
MaterialType mtype[2]
the material types of intra_grid and extra_grid grids.
Definition: electrics.h:261
std::string timer_unit(const int timer_id)
figure out units of a signal linked to a given timer
Definition: electrics.cc:772
parabolic_solver parab_solver
Solver for the parabolic bidomain equation.
Definition: electrics.h:272
double timer_val(const int timer_id)
figure out current value of a signal linked to a given timer
Definition: electrics.cc:756
void initialize()
Initialize the Electrics.
Definition: electrics.cc:36
igb_output_manager output_manager
class handling the igb output
Definition: electrics.h:281
Electrics()
Most of the initialization is done with initialize()
Definition: electrics.h:291
SF::index_mapping< mesh_int_t > petsc_to_nodal
Definition: electrics.h:220
int check_quiescence(double tm, double dt)
check for quiescence
Definition: electrics.cc:1684
void output_initial_activations()
output one nodal vector of initial activation time
Definition: electrics.cc:1793
void init(sf_vec &vm, sf_vec &phie, int offset, enum physic_t=elec_phys)
initializes all datastructs after electric solver setup
Definition: electrics.cc:1509
int check_acts(double tm)
check activations at sim time tm
Definition: electrics.cc:1629
SF::vector< Activation > acts
Definition: electrics.h:219
LAT_detector()
constructor, sets up basic datastructs from global_params
Definition: electrics.cc:1381
SF::vector< stimulus > stimuli
the electrical stimuli
Definition: electrics.h:390
elliptic_solver ellip_solver
Solver for the elliptic bidomain equation.
Definition: electrics.h:393
double timer_val(const int timer_id)
figure out current value of a signal linked to a given timer
Definition: electrics.cc:2203
std::string timer_unit(const int timer_id)
figure out units of a signal linked to a given timer
Definition: electrics.cc:2213
MaterialType mtype[2]
the material types of intra_grid and extra_grid grids.
Definition: electrics.h:388
igb_output_manager output_manager
class handling the igb output
Definition: electrics.h:395
manager for dirichlet boundary conditions
Definition: stimulate.h:193
sf_mat * phie_mat
lhs matrix to solve elliptic
Definition: electrics.h:54
void rebuild_stiffness(MaterialType *mtype, SF::vector< stimulus > &stimuli, FILE_SPEC logger)
Definition: electrics.cc:956
void rebuild_matrices(MaterialType *mtype, SF::vector< stimulus > &stimuli, FILE_SPEC logger)
Definition: electrics.cc:945
lin_solver_stats stats
Definition: electrics.h:60
sf_vec * phie_i
phi_e on intracellular grid
Definition: electrics.h:49
void solve(sf_mat &Ki, sf_vec &Vmv, sf_vec &tmp_i)
Definition: electrics.cc:1053
sf_vec * phie
phi_e
Definition: electrics.h:48
sf_sol * lin_solver
petsc or ginkgo lin_solver
Definition: electrics.h:57
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 * currtmp
temp vector for phiesrc
Definition: electrics.h:51
dbc_manager * dbc
dbcs require a dbc manager
Definition: electrics.h:63
int max_it
maximum number of iterations
Definition: electrics.h:68
sf_vec * phiesrc
I_e.
Definition: electrics.h:50
void rebuild_mass(FILE_SPEC logger)
Definition: electrics.cc:1016
sf_vec * Ivol
global Vm vector
Definition: electrics.h:107
double tol
CG stopping tolerance.
Definition: electrics.h:132
sf_vec * Iact
global Vm vector
Definition: electrics.h:108
sf_vec * Diff_term
Diffusion current.
Definition: electrics.h:123
double final_residual
Holds the residual after convergence.
Definition: electrics.h:137
sf_mat * rhs_parab
rhs matrix to solve parabolic
Definition: electrics.h:119
sf_vec * kappa_i
scaling vector for intracellular mass matrix, M
Definition: electrics.h:111
lin_solver_stats stats
Definition: electrics.h:129
void rebuild_matrices(MaterialType *mtype, limpet::MULTI_IF &miif, FILE_SPEC logger)
Definition: electrics.cc:1175
parabolic_t parab_tech
manner in which parabolic equations are solved
Definition: electrics.h:134
void solve(sf_vec &phie_i)
Definition: electrics.cc:1278
int niter
number of iterations
Definition: electrics.h:138
sf_vec * inv_mass_diag
inverse diagonal of mass matrix, for EXPLICIT solving
Definition: electrics.h:115
sf_mat * mass_i
lumped for parabolic problem
Definition: electrics.h:118
sf_vec * Ic
global Vm vector
Definition: electrics.h:106
sf_mat * u_mass_i
unscaled mass matrix, M
Definition: electrics.h:117
sf_vec * tmp_i2
scratch vector for i-grid
Definition: electrics.h:113
int max_it
maximum number of iterations
Definition: electrics.h:133
sf_mat * phie_recov_mat
rhs for phie recovery with pseudo bidomain
Definition: electrics.h:121
sf_vec * tmp_i1
scratch vector for i-grid
Definition: electrics.h:112
sf_mat * lhs_parab
lhs matrix (CN) to solve parabolic
Definition: electrics.h:120
sf_vec * Vmv
global Vm vector
Definition: electrics.h:104
sf_vec * Irhs
weighted transmembrane currents
Definition: electrics.h:114
sf_vec * old_vm
older Vm needed for 2nd order dT
Definition: electrics.h:110
sf_sol * lin_solver
petsc or ginkgo lin_solver
Definition: electrics.h:126
sf_vec * IIon
ionic currents
Definition: electrics.h:103
Electrical ionics functions and LIMPET wrappers.
physic_t
Identifier for the different physics we want to set up.
Definition: physics_types.h:51
int stimidx_from_timeridx(const SF::vector< stimulus > &stimuli, const int timer_id)
determine link between timer and stimulus
Definition: electrics.cc:786
void set_cond_type(MaterialType &m, cond_t type)
Definition: electrics.cc:838
cond_t
description of electrical tissue properties
Definition: fem_types.h:42
void compute_restr_idx_async(sf_mesh &mesh, SF::vector< mesh_int_t > &inp_idx, SF::vector< mesh_int_t > &idx)
Definition: electrics.cc:582
void apply_stim_to_vector(const stimulus &s, sf_vec &vec, bool add)
Definition: electrics.cc:457
void recover_phie_std(sf_vec &vm, phie_recovery_data &rcv)
Definition: electrics.cc:1908
@ ACT_THRESH
Definition: electrics.h:172
bool have_dbc_stims(const SF::vector< stimulus > &stimuli)
return wheter any stimuli require dirichlet boundary conditions
Definition: electrics.cc:861
void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid, SF::vector< mesh_int_t > *&restr, bool async)
Definition: electrics.cc:617
void constant_total_stimulus_current(SF::vector< stimulus > &stimuli, sf_mat &mass_i, sf_mat &mass_e, limpet::MULTI_IF *miif, FILE_SPEC logger)
Scales stimulus current to maintain constant total current across affected regions.
Definition: electrics.cc:2221
int postproc_recover_phie()
Definition: electrics.cc:1974
void balance_electrode(elliptic_solver &ellip, SF::vector< stimulus > &stimuli, int balance_from, int balance_to)
Definition: electrics.cc:393
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 compute_restr_idx(sf_mesh &mesh, SF::vector< mesh_int_t > &inp_idx, SF::vector< mesh_int_t > &idx)
Definition: electrics.cc:549
mesh_t
The enum identifying the different meshes we might want to load.
Definition: sf_interface.h:58
@ intra_elec_msh
Definition: sf_interface.h:59
void setup_phie_recovery_data(phie_recovery_data &data)
Definition: electrics.cc:2057
const char * get_tsav_ext(double time)
Definition: electrics.cc:870
SF::abstract_matrix< SF_int, SF_real > sf_mat
Definition: sf_interface.h:51
file_desc * FILE_SPEC
Definition: basics.h:138
Basic physics types.
Interface to SlimFem.
Simulator-level utility execution control functions.
Electrical stimulation functions.
event detection data structures
Definition: electrics.h:175
char * ID
ID used to name output file.
Definition: electrics.h:190
sf_vec * dvp0
additional vector for derivative
Definition: electrics.h:186
sf_vec * phi
signal
Definition: electrics.h:181
sf_vec * dvp1
additional vector for derivative
Definition: electrics.h:187
char * prv_fname
file name of previous run when restarting
Definition: electrics.h:191
int nacts
number of events detected over time step
Definition: electrics.h:194
int all
determine all or first instants of activation
Definition: electrics.h:179
ActMethod method
method to check whether activation occured
Definition: electrics.h:176
PotType measurand
quantity being monitored
Definition: electrics.h:192
FILE_SPEC fout
output file
Definition: electrics.h:188
sf_vec * tm
activation times
Definition: electrics.h:185
int init
true if intialized
Definition: electrics.h:180
int * ibuf
buffer indices where activation occured
Definition: electrics.h:183
int mode
toggle mode from standard to reverse
Definition: electrics.h:178
int offset
node number offset (used for PS)
Definition: electrics.h:193
float threshold
threshold for detection of activation
Definition: electrics.h:177
double * actbuf
buffer activation times if method==all
Definition: electrics.h:184
char * fname
output file name
Definition: electrics.h:189
sf_vec * phip
previous value of signal
Definition: electrics.h:182
description of materal properties in a mesh
Definition: fem_types.h:110
sentinel for checking activity in the tissue
Definition: electrics.h:198
bool activated
flag sentinel activation
Definition: electrics.h:199
int ID
ID of LAT detector used as sentinel.
Definition: electrics.h:203
double t_start
start of observation window
Definition: electrics.h:200
double t_window
duration of observation window
Definition: electrics.h:201
double t_quiesc
measure current duration of quiescence
Definition: electrics.h:202
File descriptor struct.
Definition: basics.h:133
for analysis of the #iterations to solve CG
Definition: timers.h:69
for analysis of the #iterations to solve CG
Definition: timers.h:13
sf_vec * dphi
Auxiliary vectors.
Definition: electrics.h:244
sf_vec * phie_rec
The phie recovery output vector buffer.
Definition: electrics.h:242
SF::vector< mesh_real_t > pts
The phie recovery locations.
Definition: electrics.h:241
SF_real gBath
Bath conductivity.
Definition: electrics.h:245