ION_IF.h

Go to the documentation of this file.

#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 <functional>
#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;        
};

class IonIfBase {
private:
  // Information about the model type.
  const IonType& _type;
  int _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::vector<LUT> _tables;         
  LUT * _tables_d = nullptr;        
  size_t _n_tables_d = 0;           
  GlobalData_t **ldata = nullptr; 
public:
  IonIfBase(const IonType& type, Target target, int num_node, const std::vector<std::reference_wrapper<IonType>>& plugins);

  virtual ~IonIfBase();

  const IonType& get_type() const;

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

  void compute(int start, int end, GlobalData_t** data);

  char* fill_buf(char *buf, int* n, opencarp::Salt_list *l) const;

  int restore(opencarp::FILE_SPEC in, int n, const int *pos, IIF_Mask_t *mask,
              size_t *offset, IMPinfo *impinfo, const int* loc2canon);

  int dump_luts(bool zipped);

  void destroy_luts();

  void tune(const char *im_par, const char *plugs, const char *plug_par);

  int read_svs(FILE* file);

  int write_svs(FILE* file, int node);

  virtual void copy_SVs_from(IonIfBase& other, bool alloc) = 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, int 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();
      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));
    }

    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 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());
      }
      else {
        memcpy(this->ion_private().data(), other.ion_private().data(), sizeof(typename T::private_type) * other.ion_private().size());
      }
    }

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

    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