electrics_eikonal.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 _EIKONAL_H
#define _EIKONAL_H

#include "SF_globals.h"
#include "electrics.h"
#include "physics_types.h"
#include "sim_utils.h"
#include "stimulate.h"

#include "sf_interface.h"
#include "timers.h"
#include "ionics.h"
#include "basics.h"

namespace opencarp
{
struct node_stats {
  enum status { in  = 1,
                out = 0 };
  enum reason { none = 0,
                nbn  = 1,
                conv = 2,
                stim = 3 };
  mesh_int_t  idX;                                                   
  mesh_int_t  cycle     = 0;                                         
  mesh_int_t  idXNB     = std::numeric_limits<Int>::min();           
  SF_real     T_A       = std::numeric_limits<double>::quiet_NaN();  
  SF_real     nbn_T_A   = std::numeric_limits<double>::quiet_NaN();  
  SF_real     T_A_      = std::numeric_limits<double>::quiet_NaN();  
  SF_real     T_R       = std::numeric_limits<double>::quiet_NaN();  
  SF_real     D_I       = std::numeric_limits<double>::quiet_NaN();  
  const char* status    = "out";                                     
  const char* reasonIn  = "-";                                       
  const char* reasonOut = "-";                                       

  FILE_SPEC logger = NULL;

  ~node_stats()
  {
    f_close(logger);
  }

  void init_logger(const char* filename);
  void log_stats(double tm, bool cflg);
  void update_status(enum status s, enum reason r);
};

class eikonal_solver
{
 public:
  eikonal_solver() : stimuliRef(nullptr) {}

  size_t num_pts;

  mesh_int_t numPtsxElem;        // number of points per element
  SF_real    time2stop_eikonal;  // time before the fim stops iterating and assess
  SF_real    CV_l;               // Variable to store the conduction velocity in each update
  mesh_int_t Index_currStim;     // index of the row in the table of stimuli and points

  const double inf = std::numeric_limits<double>::infinity();  // infinity

  // Variables for solver for triangles
  SF_real A, B, C;
  SF_real lam;
  SF_real ae;
  SF_real be;
  SF_real ce;
  SF_real det;
  SF_real Num, Frac;
  SF_real TA, TA_ori;
  SF_real p1, p2;

  // data to dump
  node_stats nodeData;

  // Variables for solver for tetrahedra
  SF_real c1, c2, c3, c4, c5, c6;
  SF_real a1, s1, s2, s3, a2, s4, b1, b2;
  SF_real R1, R2, R3, R4, R5;
  SF_real asol, bsol;
  SF_real a, b;
  SF_real derivat_a, derivat_b;

  double DI;  // diastolic interval

  // Has to be kept in line with param_globals::eik_solve
  enum Idiff_t { FOOT  = 0,
                 GAUSS = 1 };
  Idiff_t Idiff_tech = FOOT;

  struct diffusion_current {
    Idiff_t model   = FOOT;
    SF_real A_F     = 1;
    SF_real tau_F   = 1;
    SF_real V_th    = 1;
    SF_real alpha_1 = 1;
    SF_real beta_1  = 1;
    SF_real gamma_1 = 1;
    SF_real alpha_2 = 1;
    SF_real beta_2  = 1;
    SF_real gamma_2 = 1;
    SF_real alpha_3 = 1;
    SF_real beta_3  = 1;
    SF_real gamma_3 = 1;
  };

  // Vectors
  SF::vector<diffusion_current> diff_cur;
  SF::vector<SF_real>           T_A, T_R, TA_old, D_I;                          // Time of activation, time of repolarization, time of activation previously calculated, diastolic interval, final conduction velocity
  SF::vector<mesh_int_t>        List, num_changes, nReadded2List, stim_status;  // List of active nodes (L), number of changes per node to enter iterate in the list, number of times reentering the list per node, status per node after stimulus
  SF::vector<mesh_int_t>        StimulusPoints;                                 // Nodes that have initial stimulus
  SF::vector<SF_real>           StimulusTimes;                                  // Times when the initial nodes are stimulated

  SF::vector<mesh_int_t>  e2n_con;       // Connectivity vector with nodes that belong to each element
  SF::vector<mesh_int_t>  elem_start;    // For the i_th element elem_start[j] stores its starting position in the "connect" vector
  bool                    twoFib;        // true if there is sheet information, false if there is only fiber orientation

  SF::vector<SF_real> fiber_node, sheet_node;                                                          // Fiber and sheet orientation per node
  SF::vector<SF_real> x_coord, y_coord, z_coord;                                                       // x,y,z coordinates of each node
  SF::vector<SF_real> CVnodes, CVnodes_mod, p_cvrest, k_cvrest, j_cvrest, l_cvrest, cv_ar_t, cv_ar_n;  // CV restitution parameters
  sf_vec*             Idiff = nullptr;
  sf_vec*             AT    = nullptr;
  sf_vec*             RT    = nullptr;

  std::vector<mesh_int_t> n2n_connect, n2n_dsp;
  SF::vector<mesh_int_t>  e2n_cnt, n2e_cnt, n2e_con, n2e_dsp;

  SF::dmat<double> D, iD;

  SF_real actMIN = 0, actMAX = 0;

  // eikonal solver stats
  eikonal_solver_stats stats;

  bool condign1, condign2, condign3, condign4;

  ~eikonal_solver()
  {
    if (Idiff) delete Idiff;
    if (AT) delete AT;
    if (RT) delete RT;
    if (stimuliRef) delete stimuliRef;
  }

  void init();

  bool determine_model_to_run(double& time);

  void set_initial_cv();

  void set_stimuli(SF::vector<stimulus>& stimuli) { stimuliRef = &stimuli; }

  void clean_list();

  void save_eikonal_state(const char* tsav_ext);

  void update_repolarization_times(const Ionics& ion);

  void cycFIM();

  void FIM();

  void update_repolarization_times_from_rd(sf_vec& Vmv, sf_vec& Vmv_old, double time);

  void compute_diffusion_current(const double& time, sf_vec& vm);

  void compute_bc();

 private:
  // stimulus pointer
  SF::vector<stimulus>* stimuliRef = nullptr;

  void bubble_sort(SF::vector<SF_real>& Array, SF::vector<mesh_int_t>& IndexArray, int& length);

  void update_Ta_in_active_list();

  SF_real compute_H(mesh_int_t& indX, SF_real& T_A_indX);

  SF_real compute_coherence(mesh_int_t& indX);

  SF_real compute_F(mesh_int_t& indX, SF_real& T_A_indX);

  SF_real solve_tet(mesh_int_t& indx3, std::vector<mesh_int_t>& Elem_cur);

  SF_real solve_tri(mesh_int_t& indx, std::vector<mesh_int_t>& Elem_cur);

  SF_real solve_edges(mesh_int_t& indX, mesh_int_t& indY);

  SF_real mult_VtMV(SF::vector<SF_real>& X, SF::dmat<double>& iD, SF::vector<SF_real>& Y);

  SF_real compute_AT_at_node(SF_real& pp, mesh_int_t& indxf1, mesh_int_t& indyf1, mesh_int_t& indzf1);

  SF_real compute_AT_at_node_in_tet(SF_real& pf2, SF_real& qf2, mesh_int_t indxf2, mesh_int_t indyf2, mesh_int_t indzf2, mesh_int_t indwf2);

  bool add_node_neighbor_to_list(SF_real& RT, SF_real& oldTA, SF_real& newTA);

  void translate_stim_to_eikonal();

  void create_node_to_node_graph();

  void load_state_file();

  void init_diffusion_current();
};

class Eikonal : 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;

  sf_vec* phie_dummy = nullptr;

  Ionics ion;

  parabolic_solver parab_solver;

  eikonal_solver eik_solver;
  // Has to be kept in line with param_globals::eik_solve
  enum eikonal_t { EIKONAL = 0,
                   REp     = 1,
                   REm     = 2,
                   DREAM   = 3 };
  eikonal_t eik_tech = EIKONAL;

  LAT_detector lat;

  gvec_data gvec;

  igb_output_manager output_manager_time;
  igb_output_manager output_manager_cycle;

  bool do_output_eikonal = false;

  generic_timing_stats IO_stats;

  Eikonal() : ion(intra_elec_msh)
  {
    name = "Eikonal";
  }

  void initialize();

  void destroy();
  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();
  }

  ~Eikonal();

  double timer_val(const int timer_id);

  std::string timer_unit(const int timer_id);

 private:
  void solve_EIKONAL();
  void solve_RE();
  void solve_DREAM();

  void solve_RD();

  void setup_stimuli();

  void stimulate_intracellular();
  void clamp_Vm();

  void setup_solvers();

  void setup_mappings();

  void setup_output();

  void dump_matrices();

  void checkpointing();

  void prepace();

  void set_elec_tissue_properties(MaterialType* mtype, Eikonal::grid_t g, FILE_SPEC logger);
};

}  // namespace opencarp

#endif