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; 

  sf_vec* Diff_term      = 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, enum physic_t = elec_phys);

  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);

const char* get_tsav_ext(double time);

void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid,
                   SF::vector<mesh_int_t>* & restr, bool async = false);

void compute_restr_idx_async(sf_mesh & mesh,
                             SF::vector<mesh_int_t> & inp_idx,
                             SF::vector<mesh_int_t> & idx);

void compute_restr_idx(sf_mesh & mesh,
                       SF::vector<mesh_int_t> & inp_idx,
                       SF::vector<mesh_int_t> & idx);

}  // namespace opencarp

#endif