ION_IF.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#ifndef IONIC_IF_H
#define IONIC_IF_H
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <stdexcept>
#include "LUT.h"
#include "ODEint.h"
#include <math.h>
#include <assert.h>
#include <memory>
#include <functional>
#include <vector>
#include <utility>
#ifdef _OPENMP
#include <omp.h>
#endif
#ifdef HAS_CUDA_MODEL
#include <cuda_runtime.h>
#endif
#ifdef HAS_ROCM_MODEL
#include <hip/hip_runtime.h>
#endif
 
#include "basics.h"
#include "target.h"
 
#ifdef I
#undef I
#endif
 
namespace limpet {
 
#ifndef M_PI  //C99 does not define M_PI
#define M_PI  3.14159265358979323846264338327
#endif
 
//Defines used in Slava's LIMPET model '
#define heav(x) ( (x)<0 ? 0 : 1)
#define sign(x) ( (x)<0 ? -1 : 1 )
#define square(x) ((x)*(x))
#define cube(x) ((x)*(x)*(x))
 
// MPI message tags used in the IMP world
#define MPI_RESTORE_IMP_POLL_BUFSIZE_TAG  10
#define MPI_RESTORE_IMP_SNDRCV_BUFFER_TAG 11
#define MPI_DUMP_IMP_POLL_BUFSIZE_TAG     20
#define MPI_DUMP_IMP_SNDRCV_BUFFER_TAG    21
#define MPI_DUMP_SVS_POLL_BUFSIZE_TAG     30
#define MPI_DUMP_SVS_SNDRCV_BUFFER_TAG    31
 
struct tc_grp {
  float dt;
  int   skp;
  int   rat;
  int   update;
};
 
struct ts {
  int     cnt;  
  int     ng;   
  tc_grp *tcg;  
};
 
 
struct SV_TAB {
  int            svSize;    
  int            numSeg;    
  void          *y;         
};
 
#define NDEF       0
struct cell_geom {
    float      SVratio=NDEF;      
    float      v_cell=NDEF;       
    float      a_cap=NDEF;        
    float      fr_myo=NDEF;       
    float      sl_i2c=NDEF;       
};
 
}
#include "ion_type.h"
namespace limpet {
 
struct IMPinfo {
  char    *name;          
  int      sz;            
  int      nplug;         
  IMPinfo *plug;          
  bool     compatible;    
  int      map;           
  int      offset;        
  uint64_t fingerprint;   
};
 
class LutSet {
public:
  explicit LutSet(Target target) : target(target) {}
  ~LutSet();
 
  LutSet(const LutSet&) = delete;
  LutSet& operator=(const LutSet&) = delete;
 
  std::vector<LUT> host;
  LUT* device = nullptr;
  size_t device_count = 0;
 
private:
  Target target;
};
 
class IonIfBase {
private:
  // Information about the model type.
  const IonType& _type;
  node_count_t _num_node; 
  IonIfBase* _parent; 
  std::vector<IonIfBase*> _plugins; 
  uint32_t _reqdat; 
  uint32_t _moddat; 
public:
  int         miifIdx;              
  Target _target;                   
protected:
  ts          _tstp;                
private:
  cell_geom   _cgeom;               
  float       dt;                   
  std::shared_ptr<LutSet> _tables;  
  LUT* _tables_d = nullptr;          
  size_t _n_tables_d = 0;            
  GlobalData_t **ldata = nullptr; 
  struct sv_field { std::string name; int off; int sz; int type; };
  mutable std::vector<sv_field> _sv_layout; 
  mutable bool _sv_layout_valid = false;    
 
  const std::vector<sv_field>& sv_layout() const;
 
protected:
  bool has_lookup_tables() const { return this->_tables != nullptr; }
 
public:
  IonIfBase(const IonType& type, Target target, node_count_t num_node, const std::vector<std::reference_wrapper<IonType>>& plugins);
 
  virtual ~IonIfBase();
 
  const IonType& get_type() const;
 
  node_count_t get_num_node() const;
 
  std::size_t get_num_threads() const;
 
  IonIfBase* parent() const;
 
  // VERY UNRECOMMENDED TO CALL, THIS IS ONLY USED FOR AN UGLY HACK.
  void set_parent(IonIfBase* parent);
 
  std::vector<IonIfBase*>& plugins();
 
  uint32_t get_reqdat() const;
 
  uint32_t get_moddat() const;
 
  void set_moddat(uint32_t data);
 
#if defined HAS_CUDA_MODEL || defined HAS_ROCM_MODEL
  __device__ __host__
#endif
  cell_geom& cgeom() {
    return this->_cgeom;
  }
 
  float get_dt() const;
 
  void set_dt(float dt);
 
 
  ts& get_tstp();
 
  virtual void *get_sv_address() = 0;
 
  virtual std::size_t get_sv_size() const = 0;
 
//#if defined HAS_CUDA_MODEL || defined HAS_ROCM_MODEL
//  __device__ __host__
//#endif
//  SV_TAB& sv_tab() {
//    return this->_sv_tab;
//  }
 
  std::vector<LUT>& tables();
 
#if defined HAS_CUDA_MODEL || defined HAS_ROCM_MODEL
  __device__ __host__
#endif
  LUT *tables_d() const {
    return this->_tables_d;
  }
 
  size_t get_n_tables_d() const;
 
#if defined HAS_GPU_MODEL
  __device__ __host__
#endif
  Target get_target() const {
    return this->_target;
  };
 
  virtual void set_target(Target target);
 
  void initialize_params();
 
  virtual void initialize(double dt, GlobalData_t **impdat);
 
  virtual void initialize_state(double dt, GlobalData_t **impdat);
 
  void compute(node_index_t start, node_index_t end, GlobalData_t** data);
 
  char* fill_buf(char *buf, int* n, opencarp::Salt_list *l) const;
 
  node_count_t restore(opencarp::FILE_SPEC in, node_count_t n, const node_index_t *pos, IIF_Mask_t *mask,
                       size_t *offset, IMPinfo *impinfo, const global_node_index_t* loc2canon);
 
  node_count_t restore_per_node(opencarp::FILE_SPEC in, node_count_t n, const node_index_t *pos,
                                IIF_Mask_t *mask, size_t *offset, IMPinfo *impinfo,
                                const global_node_index_t* loc2canon);
 
  void get_sv_layout(std::vector<std::pair<int, int>>& fields) const;
 
  template <class Fn>
  void for_each_sv_field(node_index_t node, const std::vector<std::pair<int, int>>& fields, Fn&& fn) {
    const std::size_t vec_size = this->get_type().dlo_vector_size();
    char* block_base = static_cast<char*>(this->get_sv_address())
                     + (static_cast<std::size_t>(node) / vec_size) * this->get_sv_size();
    const std::size_t lane = static_cast<std::size_t>(node) % vec_size;
    for (const auto& f : fields)
      fn(block_base + f.first + lane * static_cast<std::size_t>(f.second),
         static_cast<std::size_t>(f.second));
  }
 
  size_t get_sv_per_node_size() const;
 
  uint64_t sv_fingerprint() const;
 
  int dump_luts(bool zipped);
 
  void destroy_luts();
 
  void rebuild_lookup_tables();
 
  void share_lookup_tables_from(const IonIfBase& other);
 
  void tune(const char *im_par, const char *plugs, const char *plug_par);
 
  int read_svs(FILE* file);
 
  int write_svs(FILE* file, node_index_t node);
 
  virtual void copy_SVs_from(IonIfBase& other, bool alloc) = 0;
 
  virtual void copy_configuration_from(IonIfBase& other) = 0;
 
  void copy_plugins_from(IonIfBase& other);
 
  void for_each(const std::function<void(IonIfBase&)>& consumer);
};
 
template<typename T>
class IonIf : public IonIfBase {
    struct has_rosenbrock_type {}; 
    struct has_rosenbrock_vector_type {}; 
    struct no_rosenbrock_type {}; 
  public:
    template<typename S>
    class LimpetArray {
      static constexpr bool is_void = std::is_void<S>::value; 
      S *_data; 
      std::size_t _size; //<! size of the SV array
      bool _allocated = false; 
      Target _target; 
      public:
 
      LimpetArray() : _size(0), _target(Target::UNKNOWN) {
      }
 
      LimpetArray(Target target, std::size_t size) : _size(size),
      _target(target) {
        this->allocate(target, size);
      }
 
      ~LimpetArray() {
        // Check if memory was allocated, if it's the case, deallocate on the
        // corresponding target
        if (!is_void && _allocated) {
          deallocate_on_target<S>(this->_target, this->_data);
        }
      }
 
      void allocate(Target target, std::size_t size) {
        this->_size = size;
        this->_target = target;
        if (!is_void){
          bool do_zero = false;
          this->_data = allocate_on_target<S>(this->_target, this->_size, do_zero);
          this->_allocated = true;
        }
      }
 
      template<typename Type = S>
      typename std::enable_if<!std::is_void<Type>::value, std::size_t>::type get_element_size() const {
        return sizeof(Type);
      }
 
      template<typename Type = S>
      typename std::enable_if<std::is_void<Type>::value, std::size_t>::type get_element_size() const {
        return 0;
      }
 
#ifdef HAS_GPU_MODEL
  __device__ __host__
#endif
      S *data() const {
        static_assert(!is_void, "The LimpetArray class can't be used with type 'void'");
        return this->_data;
      }
 
#ifdef HAS_GPU_MODEL
  __device__ __host__
#endif
      std::size_t size() const {
        static_assert(!is_void, "The LimpetArray class can't be used with type 'void'");
        return this->_size;
      }
 
      bool is_allocated() const {
        return this->_allocated;
      }
 
      LimpetArray& operator=(LimpetArray&& other) {
        if (this == &other) {
          return *this;
        }
        this->_target = other._target;
        this->_data = other._data;
        this->_size = other._size;
        if (other._allocated) {
          this->_allocated = true;
          other._allocated = false;
        }
        return *this;
      }
    };
  private:
    using SvTab = LimpetArray<typename T::state_type>; 
    using PrivateTab = LimpetArray<typename T::private_type>; 
    using PrivateVectorTab = LimpetArray<typename T::private_type_vector>; 
 
    using rosenbrock_usage = typename std::conditional<std::is_void<typename T::private_type>::value, no_rosenbrock_type, typename std::conditional<std::is_void<typename T::private_type_vector>::value, has_rosenbrock_type, has_rosenbrock_vector_type>::type>::type;
 
    // These fields are C-style arrays because they need to be accessible from device code
    SvTab _sv_tabs[Target::N_TARGETS]; 
    typename T::params_type* _params[Target::N_TARGETS] = {nullptr}; 
    PrivateTab _ion_private[Target::N_TARGETS]; 
    PrivateVectorTab _ion_private_vector; 
    size_t private_sz; 
  public:
 
    IonIf(const IonType &type, Target target, node_count_t num_node, const
        std::vector<std::reference_wrapper<IonType>>& plugins) : IonIfBase(type,
          target, num_node, plugins) {
          this->allocate_model_data();
        }
 
    ~IonIf() {
      for (int i = 0; i < Target::N_TARGETS; ++i) {
        if (this->_params[i] != nullptr) {
          deallocate_on_target((Target) i, this->_params[i]);
        }
      }
    }
 
#ifdef HAS_GPU_MODEL
    __device__ __host__
#endif
      typename T::params_type *params() const {
        return this->_params[this->get_target()];
      }
 
#ifdef HAS_GPU_MODEL
    __device__ __host__
#endif
      SvTab& sv_tab() {
        return this->_sv_tabs[this->get_target()];
      }
 
#ifdef HAS_GPU_MODEL
    __device__ __host__
#endif
    PrivateTab& ion_private() {
      return this->_ion_private[this->get_target()];
    }
 
    PrivateVectorTab& ion_private_vector() {
      if (this->get_target() != Target::MLIR_CPU) {
        throw std::logic_error("the vectorized ion private structure can only be used with the MLIR_CPU (vectorized CPU) target");
      }
      return this->_ion_private_vector;
    }
 
    void * get_sv_address() override {
      return (void *) this->_sv_tabs[this->get_target()].data();
    }
 
    std::size_t get_sv_size() const override {
      return sizeof(typename T::state_type);
    }
 
    void set_target(Target target) override {
      Target old_target = this->get_target();
      bool had_lookup_tables = this->has_lookup_tables();
      IonIfBase::set_target(target);
      this->allocate_model_data();
      memcpy(this->sv_tab().data(), this->_sv_tabs[old_target].data(), this->sv_tab().size());
      // TODO: Copy ion private data for Rosenbrock, making transformations as necessary between
      // vectorized and non-vectorized structures
      //    memcpy(this->ion_private().data(), this->_ion_private[old_target].data(), this->ion_private().size());
      memcpy(this->params(), this->_params[old_target], sizeof(typename T::params_type));
      if (had_lookup_tables && old_target != this->get_target())
        this->rebuild_lookup_tables();
    }
 
    void allocate_model_data() {
      std::size_t vec_size = this->get_type().dlo_vector_size();
      if (!this->_sv_tabs[this->get_target()].is_allocated())
        this->_sv_tabs[this->get_target()] =
          SvTab(this->get_target(), (this->get_num_node() + vec_size - 1) / vec_size);
      if (!this->_ion_private[this->get_target()].is_allocated())
        this->_ion_private[this->get_target()] = PrivateTab(this->get_target(), this->get_num_threads());
      // The private vector structure is always made for the vectorized CPU target
      if (!this->_ion_private_vector.is_allocated() && this->get_target() == Target::MLIR_CPU)
        this->_ion_private_vector = PrivateVectorTab(this->get_target(), this->get_num_threads());
      if (this->_params[this->get_target()] == nullptr)
        this->_params[this->get_target()] = allocate_on_target<typename T::params_type>(this->get_target(), 1, true);
    }
 
    void initialize(double dt, GlobalData_t **impdat) override {
      IonIfBase::initialize(dt, impdat);
      this->init_ion_private(rosenbrock_usage {});
    }
 
    void initialize_state(double dt, GlobalData_t **impdat) override {
      IonIfBase::initialize_state(dt, impdat);
      this->init_ion_private(rosenbrock_usage {});
    }
 
    void copy_configuration_from(IonIfBase& other_base) override {
      if (this->get_type() != other_base.get_type()) {
        throw std::logic_error("cannot copy IMP configuration between different model types");
      }
 
      IonIf<T>& other = static_cast<IonIf<T>&>(other_base);
      memcpy(this->params(), other.params(), sizeof(typename T::params_type));
      this->cgeom() = other.cgeom();
    }
 
    void copy_SVs_from(IonIfBase& other_base, bool alloc) override {
      IonIf<T>& other = static_cast<IonIf<T>&>(other_base);
      int tcg_size = other.get_tstp().ng * sizeof(tc_grp);
 
      if (this->get_num_node() != other.get_num_node()) {
        throw std::logic_error("cannot copy SVs if both IMPs don't handle the same amount of cells");
      }
 
      memcpy(this->sv_tab().data(), other.sv_tab().data(), sizeof(typename T::state_type) * other.sv_tab().size());
      // Only copy memory when the model actually defines a private type
      this->copy_ion_private(other, rosenbrock_usage {});
      memcpy(this->get_tstp().tcg, other.get_tstp().tcg, tcg_size);
    }
 
    void copy_ion_private(IonIf<T>& other, has_rosenbrock_vector_type) {
      if (this->get_target() == Target::MLIR_CPU) {
        memcpy(this->ion_private_vector().data(), other.ion_private_vector().data(), sizeof(typename T::private_type_vector) * other.ion_private_vector().size());
        // Preserve the copied workspace, but keep its self-reference clone-local.
        for (std::size_t i = 0; i < this->ion_private_vector().size(); ++i)
          this->ion_private_vector().data()[i].IF = this;
      }
      else {
        memcpy(this->ion_private().data(), other.ion_private().data(), sizeof(typename T::private_type) * other.ion_private().size());
        // Preserve the copied workspace, but keep its self-reference clone-local.
        for (std::size_t i = 0; i < this->ion_private().size(); ++i)
          this->ion_private().data()[i].IF = this;
      }
    }
 
    void copy_ion_private(IonIf<T>& other, has_rosenbrock_type) {
      memcpy(this->ion_private().data(), other.ion_private().data(), sizeof(typename T::private_type) * other.ion_private().size());
      // Preserve the copied workspace, but keep its self-reference clone-local.
      for (std::size_t i = 0; i < this->ion_private().size(); ++i)
        this->ion_private().data()[i].IF = this;
    }
 
    void copy_ion_private(IonIf<T>& other, no_rosenbrock_type) { }
 
    void init_ion_private(has_rosenbrock_vector_type) {
      // We need to fill the vector version of the private data for vectorized CPU
      if (this->get_target() == Target::MLIR_CPU) {
        for (std::size_t i = 0; i < this->ion_private_vector().size(); ++i) {
          this->ion_private_vector().data()[i].node_number = 0;
          this->ion_private_vector().data()[i].IF = this;
        }
      }
      else {
        for (std::size_t i = 0; i < this->ion_private().size(); ++i) {
          this->ion_private().data()[i].node_number = 0;
          this->ion_private().data()[i].IF = this;
        }
      }
    }
 
    void init_ion_private(has_rosenbrock_type) {
        for (std::size_t i = 0; i < this->ion_private().size(); ++i) {
          this->ion_private().data()[i].node_number = 0;
          this->ion_private().data()[i].IF = this;
        }
    }
 
    void init_ion_private(no_rosenbrock_type) { }
};
 
/*
   function prototypes
 */
 
void     initialize_ts(ts *tstp, int ng, int *skp, double dt);
void     update_ts(ts *tstp);
void     SV_alloc( SV_TAB *psv, int numSeg, int struct_size );
void     SV_free( SV_TAB *psv );
void     free_sv_table( void * );
void     print_IMPs(void);
bool     flag_set( const char *flags, const char *target );
void     print_models(bool );
float    modify_param( float a, char *expr );
int      process_param_mod( char *pstr, char *par, char *mod );
char*    get_typename(int type);
bool     verify_flags( const char *flags, const char* given );
int      load_ionic_module(const char*);
char    *get_next_list( char *lst, char delimiter );
 
// Functions for getting SVs in sv_init.c
SVputfcn getPutSV( SVgetfcn );
 
#include "ION_IF_sv.h"
 
#define CHANGE_PARAM( T, P, V, F )  do { \
  ((T##_Params *)P)->V = modify_param( ((T##_Params *)P)->V, F ); \
  log_msg( _nc_logf,0, 0, " %-20s modifier: %-15s value: %g",\
  #V,F,(float)((T##_Params *)P)->V);} while (0)
 
}  // namespace limpet
 
#endif