emi.h

Go to the documentation of this file.

//
// 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
// ----------------------------------------------------------------------------
 
#if WITH_EMI_MODEL
 
#ifndef _EMI_H
#define _EMI_H
 
#include "electrics.h"
#include "ionicsOnFace.h"
 
namespace opencarp {
 
class parabolic_solver_emi
{
  public:
 
  // Has to be kept in line with param_globals::parab_solve
  enum parabolic_t {SEMI_IMPLICIT = 0};
  //
  // There are three surface meshes in EMI model:
 
  // - `emi_surface_msh`:
  // Contains face/counter-face information only when the extracellular and intracellular sides are not located on the same MPI rank, i.e., when the membrane interface is split across ranks and a counter-face representation is required.
  // This mesh is associated with the vector `one_face`.
 
  // - `emi_surface_counter_msh`:
  // Contains both faces and counter-faces for all membrane interfaces,  i.e., regardless of whether the extracellular and intracellular parts are located on the same MPI rank or on different ranks.
  // This mesh is associated with the vector `both_face`.
 
  // - `emi_surface_unique_face_msh`:
  // Contains only unique faces, where each physical membrane face appears exactly once, independent of whether its two sides are located on the same MPI rank or on different ranks.
  // This mesh is associated with the vector `unique_face`.
  //
  sf_vec* ui             = nullptr; 
  sf_vec* dui            = nullptr; 
  sf_vec* ui_pre         = nullptr; 
  sf_vec* vb             = nullptr; 
  sf_vec* vb_one_face    = nullptr; 
  sf_vec* vb_both_face   = nullptr; 
  sf_vec* vb_unique_face = nullptr; 
  sf_vec* Ib             = nullptr; 
  sf_vec* Ib_one_face    = nullptr; 
  sf_vec* Ib_both_face   = nullptr; 
  sf_vec* Ib_unique_face = nullptr; 
  sf_vec* Iij_stim       = nullptr; 
  sf_vec* Iij_temp       = nullptr; 
  sf_vec* Irhs           = nullptr; 
 
  sf_mat* mass_emi       = nullptr; 
  sf_mat* mass_surf_emi  = nullptr; 
  sf_mat* lhs_emi        = nullptr; 
  sf_mat* stiffness_emi  = nullptr; 
  sf_mat* B              = nullptr; 
  sf_mat* Bi             = nullptr; 
  sf_mat* BsM            = nullptr; 
  sf_mat* Vm_myocyte_emi = nullptr; 
  
  // Direct operators (unique <-> both, skipping intermediate one-face)
  sf_mat* operator_unique_to_both_faces = nullptr; 
  sf_mat* operator_both_to_unique_face = nullptr;
 
 
  hashmap::unordered_map<SF::tuple<mesh_int_t>, 
                         std::pair<SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>,
                                   SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>>> line_face;
  hashmap::unordered_map<SF::triple<mesh_int_t>, 
                         std::pair<SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>,
                                   SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>>> tri_face;
  hashmap::unordered_map<SF::quadruple<mesh_int_t>, 
                        std::pair<SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>, 
                                  SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>>> quad_face;
 
  hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, mesh_int_t> map_vertex_tag_to_dof;
  hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, std::pair<mesh_int_t,mesh_int_t>> map_vertex_tag_to_dof_petsc;
 
  hashmap::unordered_map<mesh_int_t,mesh_int_t> dof2vertex;  
  hashmap::unordered_map<mesh_int_t,mesh_int_t> dof2petsc;
  hashmap::unordered_map<mesh_int_t,mesh_int_t> petsc2dof; 
 
  hashmap::unordered_set<int> extra_tags;
  hashmap::unordered_set<int> intra_tags;
 
  // here for any element, we set 1 for membrane and 2 for gap-junction
  SF::vector<mesh_int_t> elemTag_emi_mesh;
  SF::vector<mesh_int_t> elemTag_surface_mesh;
  SF::vector<mesh_int_t> elemTag_surface_w_counter_mesh;
 
 
  // map via vector from emi_surfmesh_w_counter_face element indices to emi_surfmesh (both -> one)
  SF::vector<mesh_int_t> vec_both_to_one_face;
 
  hashmap::unordered_map<mesh_int_t, std::pair<SF::emi_index_rank<mesh_int_t>, SF::emi_index_rank<mesh_int_t>>> map_elem_uniqueFace_to_elem_oneface;
 
  
  // Direct mapping (unique -> both face indices, no intermediate)
  hashmap::unordered_map<mesh_int_t, std::pair<SF::emi_index_rank<mesh_int_t>, SF::emi_index_rank<mesh_int_t>>> map_elem_uniqueFace_to_elem_bothface;
 
  hashmap::unordered_map<mesh_int_t, std::pair<mesh_int_t, mesh_int_t>> map_elem_uniqueFace_to_tags;
  // 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;                  
  parabolic_t  parab_tech        = SEMI_IMPLICIT;         
 
  // solver statistics
  double       final_residual    = -1.0;                 
  int          niter             = -1;                   
 
  ~parabolic_solver_emi()
  {
    if(dbc) delete dbc;
    if (lin_solver) delete lin_solver;
    // matrices
    if (mass_emi) delete mass_emi;
    if (mass_surf_emi) delete mass_surf_emi;
    if (lhs_emi) delete lhs_emi;
    if (stiffness_emi) delete stiffness_emi;
    if (B) delete B;
    if (Bi) delete Bi;
    if (BsM) delete BsM;
    if (Vm_myocyte_emi) delete Vm_myocyte_emi;
 
    if (operator_unique_to_both_faces) delete operator_unique_to_both_faces; 
    if (operator_both_to_unique_face) delete operator_both_to_unique_face; 
 
    // vectors
    if (ui) delete ui;
    if (dui) delete dui;
    if (ui_pre) delete ui_pre;
    if (vb) delete vb;
 
    if (vb_both_face) delete vb_both_face;
    if (vb_unique_face) delete vb_unique_face;
    if (Ib) delete Ib;
 
    if (Ib_both_face) delete Ib_both_face;
    if (Ib_unique_face) delete Ib_unique_face;
    if (Iij_stim) delete Iij_stim;
    if (Iij_temp) delete Iij_temp;
    if (Irhs) delete Irhs;
 
    line_face.clear();
    tri_face.clear();
    quad_face.clear();
  }
 
  void init();
  void rebuild_matrices(MaterialType* mtype, limpet::MULTI_IF & miif, SF::vector<stimulus> & stimuli, FILE_SPEC logger);
  void solve();
 
  private:
  void setup_linear_solver(FILE_SPEC logger);
 
  void solve_semiImplicit();
};
 
class EMI : public Basic_physic
{
  public:
    
  MaterialType mtype_vol[1];
  MaterialType_EMI mtype_face;
  SF::vector<stimulus> stimuli;
 
  IonicsOnFace ion;
 
  parabolic_solver_emi parab_solver;
 
  gvec_data_OnFace gvec;
 
  igb_output_manager output_manager;
 
  phie_recovery_data phie_rcv;
 
  generic_timing_stats IO_stats;
        
  hashmap::unordered_map<mesh_int_t,mesh_int_t> vertex2ptsdata; 
  hashmap::unordered_map<mesh_int_t,mesh_int_t> dof2ptsData; 
                                           
  EMI() : ion(emi_surface_unique_face_msh)
  {
    name = "EMI"; 
  }
 
  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();
  }
 
  ~EMI()
  {
  };
 
  double timer_val(const int timer_id);
 
  std::string timer_unit(const int timer_id);
 
  private:
   void setup_stimuli();
 
   void apply_current_stimulus();
 
   void apply_dbc_stimulus();
 
   void balance_electrodes();
 
   void scale_total_stimulus_current(SF::vector<stimulus>& stimuli,
                                     sf_mat&               mass_vol,
                                     sf_mat&               mass_surf,
                                     FILE_SPEC             logger);
 
   void setup_solvers();
 
   void setup_mappings();
 
   void setup_output();
 
   void dump_matrices();
 
   void checkpointing();
 
   void setup_EMI_mesh();
}; 
 
void log_mesh_local_element_ranges(const sf_mesh& emi_mesh,
                                   const sf_mesh& emi_surfmesh_w_counter_face,
                                   const sf_mesh& emi_surfmesh_unique_face);
void set_elec_tissue_properties_emi_volume(MaterialType* mtype, FILE_SPEC logger);
 
void extract_unique_tag(SF::vector<mesh_int_t>& unique_tags);
 
void compute_tags_per_rank(int num_tags, SF::vector<mesh_int_t>& num_tags_per_rank);
 
void distribute_elements_based_tags(SF::meshdata<mesh_int_t, mesh_real_t>& mesh,
                                    int total_num_tags);
 
void partition_based_tags(int num_tags,
                          SF::vector<mesh_int_t> tag,
                          SF::vector<mesh_int_t> num_tags_per_rank,
                          SF::vector<mesh_int_t>& part);
 
void map_tags_to_rank(int size, const SF::vector<mesh_int_t> & unique_tags, const SF::vector<mesh_int_t> & num_tags_per_rank, hashmap::unordered_map<mesh_int_t, mesh_int_t> &tags_to_rank_map);
 
bool load_partitions_from_file(hashmap::unordered_map<mesh_int_t, mesh_int_t>& tags_to_rank_map,
                               int expected_num_tags,
                               MPI_Comm comm);
 
void permute_mesh_locally_based_on_tag_elemIdx(SF::meshdata<mesh_int_t, mesh_real_t>& mesh);
 
}  // namespace opencarp
 
 
#endif
#endif