electrics.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#ifndef _ELECTRICS_H
#define _ELECTRICS_H
 
#include "physics_types.h"
#include "sim_utils.h"
#include "stimulate.h"
 
#include "sf_interface.h"
#include "timers.h"
#include "ionics.h"
 
namespace opencarp {
 
class elliptic_solver
{
  public:
  sf_vec* phie       = nullptr; 
  sf_vec* phie_i     = nullptr; 
  sf_vec* phiesrc    = nullptr; 
  sf_vec* currtmp    = nullptr; 
 
  sf_mat* mass_e     = nullptr; 
  sf_mat* phie_mat   = nullptr; 
 
  // the linear solver
  sf_sol* lin_solver = nullptr; 
 
  // linear solver stats
  lin_solver_stats stats;
 
  dbc_manager* dbc                    = nullptr;
  bool         phie_mat_has_nullspace = false;
 
  // solver config
  double       tol            = 1e-8;                 
  int          max_it         = 100;                  
 
  void init();
  void rebuild_matrices(MaterialType* mtype, SF::vector<stimulus> & stimuli, FILE_SPEC logger);
  void rebuild_stiffness(MaterialType* mtype, SF::vector<stimulus> & stimuli, FILE_SPEC logger);
  void rebuild_mass(FILE_SPEC logger);
  void solve(sf_mat & Ki, sf_vec & Vmv, sf_vec & tmp_i);
  void solve_laplace();
 
  ~elliptic_solver()
  {
    if(dbc) delete dbc;
    if(lin_solver) delete lin_solver;
    // matrices
    if(phie_mat) delete phie_mat;
    if(mass_e) delete mass_e;
    // vectors
    if(phie) delete phie;
    if(phie_i) delete phie_i;
    if(phiesrc) delete phiesrc;
    if(currtmp) delete currtmp;
  }
 
  private:
  void enforce_dbc();
  void setup_linear_solver(FILE_SPEC logger);
};
 
class parabolic_solver
{
  public:
 
  // Has to be kept in line with param_globals::parab_solve
  enum parabolic_t { EXPLICIT = 0, CN = 1, O2dT = 2};
 
  sf_vec* IIon           = nullptr; 
  sf_vec* Vmv            = nullptr; 
 
  sf_vec* old_vm         = nullptr; 
  sf_vec* kappa_i        = nullptr; 
  sf_vec* tmp_i1         = nullptr; 
  sf_vec* tmp_i2         = nullptr; 
  sf_vec* Irhs           = nullptr; 
  sf_vec* inv_mass_diag  = nullptr; 
 
  sf_mat* u_mass_i       = nullptr; 
  sf_mat* mass_i         = nullptr; 
  sf_mat* rhs_parab      = nullptr; 
  sf_mat* lhs_parab      = nullptr; 
  sf_mat* phie_recov_mat = nullptr; 
 
  // the linear solver
  sf_sol* lin_solver     = nullptr; 
 
  // linear solver stats
  lin_solver_stats stats;
 
  // solver config
  double       tol               = 1e-8;                 
  int          max_it            = 100;                  
  parabolic_t  parab_tech        = CN;                   
 
  // solver statistics
  double       final_residual    = -1.0;                 
  int          niter             = -1;                   
 
  ~parabolic_solver()
  {
    if (lin_solver) delete lin_solver;
    // matrices
    if (u_mass_i) delete u_mass_i;
    if (mass_i) delete mass_i;
    if (rhs_parab) delete rhs_parab;
    if (lhs_parab) delete lhs_parab;
    if (phie_recov_mat) delete phie_recov_mat;
    // vectors
    // IIon and Vmv are shallow copies of global vectors, do not delete them
    if (old_vm) delete old_vm;
    if (kappa_i) delete kappa_i;
    if (tmp_i1) delete tmp_i1;
    if (tmp_i2) delete tmp_i2;
    if (Irhs) delete Irhs;
    if (inv_mass_diag) delete inv_mass_diag;
  }
 
  void init();
  void rebuild_matrices(MaterialType* mtype, limpet::MULTI_IF & miif, FILE_SPEC logger);
  void solve(sf_vec & phie_i);
 
  private:
  void setup_linear_solver(FILE_SPEC logger);
 
  void solve_CN(sf_vec & phie_i);
  void solve_O2dT(sf_vec & phie_i);
  void solve_EF(sf_vec & phie_i);
};
 
enum PotType {VM, PHIE};
enum ActMethod {NONE, ACT_THRESH, ACT_DT};
 
struct Activation {
  ActMethod  method;              
  float      threshold;           
  int        mode;                
  int        all;                 
  int        init;                
  sf_vec*    phi       = nullptr; 
  sf_vec*    phip      = nullptr; 
  int       *ibuf      = nullptr; 
  double    *actbuf    = nullptr; 
  sf_vec*    tm        = nullptr; 
  sf_vec*    dvp0      = nullptr; 
  sf_vec*    dvp1      = nullptr; 
  FILE_SPEC  fout;                
  char      *fname     = nullptr; 
  char      *ID        = nullptr; 
  char      *prv_fname = nullptr; 
  PotType    measurand;           
  int        offset;              
  int        nacts;               
};
 
struct Sentinel {
  bool           activated = false;  
  double         t_start  = -1.0;    
  double         t_window = -1.0;    
  double         t_quiesc = -1.0;    
  int            ID = -1;            
};
 
class LAT_detector
{
  private:
  int check_cross_threshold(sf_vec & vm, sf_vec & vmp, double tm,
                            int *ibuf, double *actbuf, float threshold, int mode);
 
  int check_mx_derivative(sf_vec & vm, sf_vec & vmp, double tm,
                          int *ibuf, double *actbuf, sf_vec & dvp0, sf_vec & dvp1,
                          float threshold, int mode);
 
  public:
  SF::vector<Activation>        acts;
  SF::index_mapping<mesh_int_t> petsc_to_nodal;
  Sentinel sntl;
 
  LAT_detector();
 
  void init(sf_vec & vm, sf_vec & phie, int offset);
 
  int check_acts(double tm);
 
  int check_quiescence(double tm, double dt);
 
  void output_initial_activations();
};
 
struct phie_recovery_data
{
  SF::vector<mesh_real_t> pts; 
  sf_vec* phie_rec = nullptr;  
  sf_vec* Im = nullptr;
  sf_vec* dphi = nullptr;      
  SF_real gBath;               
};
 
class Electrics : public Basic_physic
{
  public:
  enum grid_t {intra_grid = 0, extra_grid};
  enum linsys_t {elliptic_sys, parabolic_sys};
 
  MaterialType mtype[2];
  SF::vector<stimulus> stimuli;
 
  Ionics ion;
 
  elliptic_solver  ellip_solver;
  parabolic_solver parab_solver;
 
  LAT_detector lat;
 
  gvec_data gvec;
 
  igb_output_manager output_manager;
 
  phie_recovery_data phie_rcv;
 
  generic_timing_stats IO_stats;
 
  Electrics() : ion(intra_elec_msh)
  {
    name = "Electrics";
  }
 
  void initialize();
 
  void destroy();
 
  // This funcs from the Basic_physic interface are currently empty
  void compute_step();
  void output_step();
 
  inline void output_timings()
  {
    // since ionics are included in electrics, we substract ionic timings from electric timings
    compute_time    -= ion.compute_time;
    initialize_time -= ion.initialize_time;
 
    // now we call the base class timings output for Electrics
    Basic_physic::output_timings();
    // and then for ionics
    ion.output_timings();
  }
 
  ~Electrics();
 
  double timer_val(const int timer_id);
 
  std::string timer_unit(const int timer_id);
 
 
  private:
 
  void setup_stimuli();
 
  void stimulate_intracellular();
  void clamp_Vm();
 
  void stimulate_extracellular();
 
  void setup_solvers();
 
  void setup_mappings();
 
  void setup_output();
 
  void dump_matrices();
 
  void checkpointing();
 
  void balance_electrodes();
 
  void prepace();
};
 
void set_elec_tissue_properties(MaterialType* mtype, Electrics::grid_t g, FILE_SPEC logger);
 
bool have_dbc_stims(const SF::vector<stimulus> & stimuli);
 
int stimidx_from_timeridx(const SF::vector<stimulus> & stimuli, const int timer_id);
 
void setup_phie_recovery_data(phie_recovery_data & data);
void recover_phie_std(sf_vec & vm, phie_recovery_data & rcv);
int  postproc_recover_phie();
 
class Laplace : public Basic_physic
{
  public:
  MaterialType mtype[2];
  SF::vector<stimulus> stimuli;
 
  elliptic_solver  ellip_solver;
  igb_output_manager output_manager;
 
  Laplace()
  {
    name = "Laplace solver";
  }
 
  void initialize();
  void destroy();
  void compute_step();
  void output_step();
  double timer_val(const int timer_id);
  std::string timer_unit(const int timer_id);
};
 
void constant_total_stimulus_current(SF::vector<stimulus> & simuli,
                                     sf_mat & mass_i,
                                     sf_mat & mass_e,
                                     limpet::MULTI_IF *miif,
                                     FILE_SPEC logger);
 
void balance_electrode(SF::vector<stimulus> & stimuli, int balance_from, int balance_to);
 
void apply_stim_to_vector(const stimulus & s, sf_vec & vec, bool add);
 
void set_cond_type(MaterialType & m, cond_t type);
}  // namespace opencarp
 
#endif