electrics_eikonal.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#pragma once
 
#ifndef _EIKONAL_H
#define _EIKONAL_H
 
#include <array>
 
#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) {}
 
  const double inf = std::numeric_limits<double>::infinity();
 
  // data to dump
  node_stats nodeData;
 
  // persistent handle for the detailed per-node debug logger (dream.output.detailed_debugNode),
  // opened once in init() and closed in the destructor instead of per debug call
  FILE_SPEC deepDebugLogger = NULL;
 
  // cached once in init() so LogEikonalCandidate call sites can skip building their
  // arguments entirely when no debug node is configured, rather than allocating them
  // just to have LogEikonalCandidate itself immediately discard them
  bool detailed_debug_enabled = false;
 
  // 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;
  };
 
  // enum classes used to define status related to stimulus of each node
  enum class StimStatus : mesh_int_t {
    None         = 0,  // no stimulus attempted
    Success      = 1,  // successfully stimulated
    Failed       = 2,  // stimulus attempted but unsuccessful: node not recovered
    Subthreshold = 3   // node recovered, but its cell is predicted not to fire
  };
 
  // A pulse captures only if the predicted upstroke starts within this time after its onset
  static constexpr double STIM_CAPTURE_WINDOW = 30.0;  // ms
 
  // Stimulus table rows identify a node's pulse (onset, strength, duration). With NO_PULSE as its
  // row, predict_upstroke() integrates the node's cell without a pulse. This happens after a
  // checkpoint load for a node whose non-capturing pulse lies before the save time: the
  // checkpoint stores only the pulse onset, and the pulse's effect is already in the loaded
  // cell state.
  static constexpr size_t NO_PULSE = static_cast<size_t>(-1);
 
  // enum classes used to define the type of repolarization times to find
  enum class RepolMode {
    Normal,
    FailedStim
  };
 
  // classification of a freshly-converged node after post-sweep FIM verification
  enum class VerifyResult {
    Verified,       // sweep T_A is consistent with converged-only evaluation
    Overestimated,  // sweep T_A was inflated; re-enter node with corrected value
    Unresolved      // converged-only evaluation yielded no activation time; a genuine
                    // conduction block and not-yet-settled neighbors are indistinguishable
                    // at this point, so the node is kept active either way
  };
 
  // record produced by post-sweep Phase 1 verification of a freshly-converged node
  struct VerifyEntry {
    mesh_int_t   indX;         // global node index
    SF_real      ta_verified;  // T_A from converged-only evaluation; inf if none could be formed
    VerifyResult result;       // convergence quality classification
  };
 
  // 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;                                                                                                                             // List of active nodes (L), number of changes per node to enter iterate in the list, number of times reentering the list per node
  SF::vector<mesh_int_t>        StimulusPoints;                                                                                                                                               // Nodes that have initial stimulus
  std::vector<bool>             is_stimulus_target;                                                                                                                                           // per node: does it appear anywhere in StimulusPoints?
  SF::vector<StimStatus>        stim_status;                                                                                                                                                  // status per node after stimulus
  SF::vector<SF_real>           StimulusTimes, StimulusDurations, StimulusStrengths, stim_strength, stim_duration, stim_start;                                                                // Times when the initial nodes are stimulated
  std::vector<long>             capture_candidate;                                                                                                                                            // per node: index into the stimulus rows whose capture the current compute_bc() call predicts (-1: none)
  std::vector<size_t>           last_noncapturing_row;                                                                                                                                        // per node: stimulus table row of that pulse (NO_PULSE: none, or it lies before the loaded checkpoint)
  SF::vector<SF_real>           last_noncapturing_stim;                                                                                                                                       // per node: onset of the last stimulus pulse that failed or was subthreshold at the node (-inf: none)
  SF::vector<SF_real>           rejected_TA, rejected_at;                                                                                                                                     // per node: the T_A reject_non_firing_activations() last turned down and the RD time it did so (-inf: none). Within one RD step the verdict cannot change, so Phase 2A does not propose that T_A again
  SF::vector<SF_real>           captured_onset, captured_TA;                                                                                                                                  // onset and T_A of the stimulus pulse that last captured the node (-inf: none); Idiff is suppressed there from the onset to T_A + 5 ms (solve_RD)
  SF::vector<int>               StimulusOrigIndex;                                                                                                                                            // logical stimulus index (into *stimuliRef) for each row of the tables above
  std::vector<bool>             warned_failed_capture;                                                                                                                                        // per logical stimulus: already warned that it failed to capture a target node?
  std::vector<bool>             warned_zero_strength;                                                                                                                                         // per logical stimulus: already warned that its pulse.strength is 0?
  std::vector<bool>             stim_failed_now;                                                                                                                                              // per node: did the current compute_bc() call fail to capture it?
  std::vector<mesh_int_t>       stim_failed_nodes;                                                                                                                                            // the nodes stim_failed_now is set for, so clearing it costs one pass over those
  SF::vector<SF_real>           failed_stim_time, T_R_before_failed_stim;                                                                                                                     // per node: earliest failed stimulus still ahead of the clock (-inf: none) and T_R before it, see save_eikonal_state
  bool                          warned_cycle_stall   = false;                                                                                                                                 // already warned that solve_DREAM's cycle loop made no progress?
  bool                          warned_sweep_stall   = false;                                                                                                                                 // already warned that cycFIM's list sweep made no progress?
  long long                     num_past_activations = 0;                                                                                                                                     // activation times set before the current RD time (see check_past_activation)
  size_t                        Index_currStim;                                                                                                                                               // index of the row in the table of stimuli and points
  mesh_int_t                    num_active_nodes = 0;                                                                                                                                         // number of nodes in the active List
  SF_real                       actMIN = std::numeric_limits<double>::quiet_NaN(), actMAX = std::numeric_limits<double>::quiet_NaN(), maxadvance = std::numeric_limits<double>::quiet_NaN();  // min/max T_A in active list (wavefront bounds); maxadvance = global FIM upper bound
  SF_real                       loaded_save_time = -1.0;                                                                                                                                      // Time at which checkpoint was saved (-1 if no checkpoint loaded)
 
  // mesh related data
  mesh_int_t              MESH_SIZE;  // number of points per element
  bool                    twoFib;     // true if there is sheet information, false if there is only fiber orientation
  size_t                  num_pts;
  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
  std::vector<mesh_int_t> n2n_connect, n2n_dsp;
  SF::vector<mesh_int_t>  e2n_cnt, n2e_cnt, n2e_con, n2e_dsp;
 
  // Element data storage for anisotropy and CV
  std::vector<SF::dmat<double>> S;
  std::vector<double>           CV_L;
 
  // CV restitution parameters
  SF::vector<SF_real> rho_cvrest, kappa_cvrest, theta_cvrest, denom_cvrest;
 
  // export vectors
  sf_vec* Idiff = nullptr;
  sf_vec* AT    = nullptr;
  sf_vec* RT    = nullptr;
 
  // eikonal solver statistics
  eikonal_solver_stats stats;
 
  // enum classes used in the debugger file
  enum class DebugCase {
    ListIter,
    ListReenter,
    NBStimulusEnter,
    FailedReenter,
    Stimulated,
    UnsuccessfulStimulated,
    RejectedNonFiring,
    ExitsListStimulated,
    ExitsListInvalidUpdate,
    MaskedUntrusted,
    UnmaskedUntrusted,
    Unknown
  };
 
  // Which element a candidate activation time was formed on. SELECTED instead marks the
  // local solver's chosen candidate, whose element is read from the contributor count.
  enum class ElementCandidate {
    SELECTED,
    TETRAHEDRON,
    TRIANGLE,
    EDGE,
    UNKNOWN
  };
 
  ~eikonal_solver()
  {
    if (Idiff) delete Idiff;
    if (AT) delete AT;
    if (RT) delete RT;
    if (stimuliRef) delete stimuliRef;
    f_close(deepDebugLogger);
  }
 
  void init();
 
  bool should_solve_reaction_diffusion(double& time);
 
  void set_stimuli(SF::vector<stimulus>& stimuli) { stimuliRef = &stimuli; }
 
  void clean_list();
 
  void save_eikonal_state(const char* roe_fnm);
 
  void update_repolarization_times(const Ionics& ion, RepolMode mode);
 
  void cycFIM(const Ionics& ion);
 
  void FIM();
 
  void deep_debug_log(DebugCase debug_case, mesh_int_t index_ref = 0, mesh_int_t ind_target = 0);
 
  void LogEikonalCandidate(ElementCandidate type, mesh_int_t ind_target, const int* contributors, size_t num_contributors, mesh_real_t candidate_TA, mesh_real_t candidate_scaling_CV, int counter);
 
  void LogPhaseOneVerify(mesh_int_t indX, SF_real ta_sweep, SF_real ta_verified, VerifyResult result, bool forceConverge);
 
  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 init_imp_region_properties(FILE_SPEC logger);
 
 private:
  // stimulus pointer
  SF::vector<stimulus>* stimuliRef = nullptr;
 
  void set_successful_stimulus(mesh_int_t node, SF_real TA);
 
  void predict_upstroke(const Ionics& ion, const std::vector<size_t>& rows, std::vector<SF_real>& t_up, const std::vector<SF_real>* idiff_TA = nullptr, std::vector<SF_real>* vm_end = nullptr);
 
  void reject_non_firing_activations(const Ionics& ion, std::vector<mesh_int_t>& converged);
 
  SF_real idiff_template(const diffusion_current& node, SF_real dT, SF_real vm) const;
 
  void check_past_activation(mesh_int_t indX);
 
  void set_active(mesh_int_t i);
 
  void clear_active(mesh_int_t i);
 
  void mask_untrusted(mesh_int_t i);
 
  void unmask_untrusted(mesh_int_t i);
 
  void set_unsuccessful_stimulus(mesh_int_t node, size_t stimIndex, SF_real stimTime);
 
  void warn_zero_strength_stimulus(int stim_idx, SF_real stimTime);
 
  void activate_neighbors(mesh_int_t& node);
 
  void compute_bc(const Ionics& ion, bool EmptyList, SF_real time_upper_bound);
 
  void precompute_squared_anisotropy_metric();
 
  // --- Small helpers for active list management ---
  inline void add_to_active(std::vector<mesh_int_t>& active,
                            std::vector<char>&       in_active,
                            mesh_int_t               node)
  {
    if (!in_active[node]) {
      in_active[node] = 1;
      active.push_back(node);
    }
  }
 
  inline void remove_from_active(std::vector<char>& in_active, mesh_int_t node)
  {
    in_active[node] = 0;
  }
 
  void update_Ta_in_active_list();
 
  SF_real compute_H(mesh_int_t& indX, SF_real& tmpTA);
 
  void record_D_I(mesh_int_t indX);
 
  SF_real compute_coherence(mesh_int_t& indX, bool convergedOnly = false);
 
  SF_real update(mesh_int_t& indX, SF_real CVrest_factor = 1.0, bool isDREAM = false, bool convergedOnly = false);
 
  template <int N>
  SF_real update_impl(mesh_int_t& indX, SF_real CVrest_factor = 1.0, bool CheckValidity = false, bool convergedOnly = false);
 
  template <int N>
  bool is_not_valid_update(const std::array<int, N>& nodeIDs, double time, int e_i, double cv, bool convergedOnly = false);
 
  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();
};
 
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 set_elec_tissue_properties(MaterialType* mtype, Eikonal::grid_t g, FILE_SPEC logger);
};
 
template <int MESH_SIZE>
class LocalSolver
{
 private:
  SF::dmat<double>&                 D;
  std::array<SF::Point, MESH_SIZE>& points;
  std::array<double, MESH_SIZE>&    values;
 
  double tsitsiklis_update_line(const std::array<SF::Point, 2>& base,
                                const SF::dmat<double>&         D,
                                const double&                   value);
  double tsitsiklis_update_triangle(const std::array<SF::Point, 3>& base,
                                    const SF::dmat<double>&         D,
                                    const std::array<double, 2>&    values);
  double tsitsiklis_update_tetra(const std::array<SF::Point, 4>& base,
                                 const SF::dmat<double>&         D,
                                 const std::array<double, 3>&    values);
 
 public:
  LocalSolver(SF::dmat<double>& D, std::array<SF::Point, MESH_SIZE>& points, std::array<double, MESH_SIZE>& values)
      : D(D), points(points), values(values) {}
 
  // This function performs the local solver, returning the solution value for the last point with reference to
  // point order of constructor input
  double solve();
};
 
template class LocalSolver<2>;
 
template class LocalSolver<3>;
 
template class LocalSolver<4>;
 
}  // namespace opencarp
 
#endif