grids.h

Go to the documentation of this file.

/******************************************************************
Author: Austin Lachance
Date: 10/28/16
Description: Header File for grids.cpp. Allows access to Grid_node
and Flux_pair structs and their respective functions
******************************************************************/
#include <stdio.h>
#include <assert.h>
#include <nrnmpi.h>

#include <nrnwrap_Python.h>

#define SAFE_FREE(ptr){if((ptr)!=NULL) free(ptr);}
#define IDX(x,y,z)  ((z) + (y) * g->size_z + (x) * g->size_z * g->size_y)
#define INDEX(x,y,z)  ((z) + (y) * grid->size_z + (x) * grid->size_z * grid->size_y)
#define ALPHA(x,y,z) (g->get_alpha(g->alpha,IDX(x,y,z)))
#define VOLFRAC(idx) (g->get_alpha(g->alpha,idx))
#define TORT(idx) (g->get_permeability(g->permeability,idx))
#define PERM(x,y,z) (g->get_permeability(g->permeability,IDX(x,y,z)))
#define SQ(x)       ((x)*(x))
#define CU(x)       ((x)*(x)*(x))
#define TRUE                1
#define FALSE               0
#define TORTUOSITY          2
#define VOLUME_FRACTION     3
#define ICS_ALPHA           4

#define NEUMANN             0
#define DIRICHLET           1

#define MAX(a,b)    ((a)>(b)?(a):(b))
#define MIN(a,b)    ((a)<(b)?(a):(b))

typedef struct {
    PyObject_HEAD
    void* ho_;
    union {
        double x_;
        char* s_;
        void* ho_;
        double* px_;
    }u;
    void* sym_; // for functions and arrays
    void* iteritem_; // enough info to carry out Iterator protocol
    int nindex_; // number indices seen so far (or narg)
    int* indices_; // one fewer than nindex_
    int type_; // 0 HocTopLevelInterpreter, 1 HocObject
        // 2 function (or TEMPLATE)
        // 3 array
        // 4 reference to number
        // 5 reference to string
        // 6 reference to hoc object
        // 7 forall section iterator
        // 8 pointer to a hoc scalar
        // 9 incomplete pointer to a hoc array (similar to 3)
} PyHocObject;

typedef struct Hybrid_data {
    long num_1d_indices;
    long* indices1d;
    long* num_3d_indices_per_1d_seg;
    long* indices3d;
    double* rates;
    double* volumes1d;
    double* volumes3d;
} Hybrid_data;

typedef struct Flux_pair {
    double *flux;           // Value of flux
    int grid_index;         // Location in grid
} Flux;

typedef struct {
    double* destination;    /* memory loc to transfer concentration to */
    long source;            /* index in grid for source */
} Concentration_Pair;

typedef struct {
    long destination;       /* index in grid */
    double* source;         /* memory loc of e.g. ica */
    double scale_factor;
} Current_Triple;

typedef void (*ReactionRate)(double**, double**, double**, double*, double*, double*, double*, double**, double);
typedef void (*ECSReactionRate)(double*, double*, double*, double*);
typedef struct Reaction {
    struct Reaction* next;
    ECSReactionRate reaction;
    unsigned int num_species_involved;
    unsigned int num_params_involved;
    double** species_states;
    unsigned char* subregion;
    unsigned int region_size;
    uint64_t* mc3d_indices_offsets;
    double** mc3d_mults;
} Reaction;

typedef struct {
        double* copyFrom;
        long copyTo;
} AdiLineData;

typedef struct {
    /*TODO: Support for mixed boundaries and by edge (maybe even by voxel)*/
    unsigned char type;
    double value;
} BoundaryConditions; 

class Grid_node {
    public:
    Grid_node *next;

    double *states;         // Array of doubles representing Grid space
    double *states_x;
    double *states_y;
    double *states_z;       //TODO: This is only used by ICS, is it necessary? 
    double *states_cur;
    int size_x;          // Size of X dimension
    int size_y;          // Size of Y dimension
    int size_z;          // Size of Z dimension
    double dc_x;            // X diffusion constant
    double dc_y;            // Y diffusion constant
    double dc_z;            // Z diffusion constant
    double dx;              // ∆X
    double dy;              // ∆Y
    double dz;              // ∆Z
    bool diffusable;
    bool hybrid;
    BoundaryConditions* bc;
    Hybrid_data* hybrid_data;
    Concentration_Pair* concentration_list;
    Current_Triple* current_list;
    ssize_t num_concentrations, num_currents;

    /*used for MPI implementation*/
    int num_all_currents;
    int* proc_offsets;
    int* proc_num_currents;
    int* proc_flux_offsets;
    int* proc_num_fluxes;
    long* current_dest;
    double* all_currents;

    /*Extension to handle a variable diffusion characteristics of a grid*/
    unsigned char   VARIABLE_ECS_VOLUME;    /*FLAG which variable volume fraction
                                            methods should be used*/
    /*diffusion characteristics are arrays of a single value or
    * the number of voxels (size_x*size_y*size_z)*/
    double *permeability;       /* 1/tortuosities^2 squared 
                                   D_eff = D_free*permeability */
    double *alpha;      /* volume fractions */
    /*Function that will be assigned when the grid is created to either return
    * the single value or the value at a given index*/ 
    double (*get_alpha)(double*,int);
    double (*get_permeability)(double*,int);
    double atolscale;

    int64_t* ics_surface_nodes_per_seg;
    int64_t* ics_surface_nodes_per_seg_start_indices;
    double** ics_concentration_seg_ptrs;
    double** ics_current_seg_ptrs;
    double* ics_scale_factors;
    int ics_num_segs;

    int insert(int grid_list_index);
    int node_flux_count;
    long * node_flux_idx;
    double * node_flux_scale;
    PyObject ** node_flux_src;

    
    virtual ~Grid_node(){};
    virtual void set_diffusion(double*, int) = 0;
    virtual void set_num_threads(const int n) = 0;
    virtual void do_grid_currents(double*, double, int) = 0;
    virtual void apply_node_flux3D(double dt, double* states) = 0;
    virtual void volume_setup() = 0;
    virtual int dg_adi() = 0;
    virtual void variable_step_diffusion(const double* states, double* ydot) = 0;
    virtual void variable_step_ode_solve(double* RHS, double dt) = 0;
    virtual void scatter_grid_concentrations() = 0;
    virtual void hybrid_connections() = 0;
    virtual void variable_step_hybrid_connections(const double* cvode_states_3d, double* const ydot_3d, const double* cvode_states_1d, double *const  ydot_1d) = 0;
};

class ECS_Grid_node : public Grid_node{
    public:
        //Data for DG-ADI
        ECS_Grid_node();
        ECS_Grid_node(PyHocObject*, int, int, int, double, double, double, 
                      double, double, double, PyHocObject*, PyHocObject*, int,
                      double, double);
        ~ECS_Grid_node();
        struct ECSAdiGridData* ecs_tasks;
        struct ECSAdiDirection* ecs_adi_dir_x;
        struct ECSAdiDirection* ecs_adi_dir_y;
        struct ECSAdiDirection* ecs_adi_dir_z;
        
        //Data for multicompartment reactions
        int induced_idx;
        int* react_offsets;
        int react_offset_count;
        int* reaction_indices;
        int* all_reaction_indices;
        int* proc_num_reactions;
        int* proc_num_reaction_states;
        int total_reaction_states;
        unsigned char multicompartment_inititalized;
        int* induced_currents_index;
        int induced_current_count;
        int* proc_induced_current_count;
        int* proc_induced_current_offset;
        double* all_reaction_states;
        double* induced_currents;
        double* local_induced_currents;
        double* induced_currents_scale;

        void set_num_threads(const int n);
        void do_grid_currents(double *, double, int);
        void apply_node_flux3D(double dt, double* states);
        void volume_setup();
        int dg_adi();
        void variable_step_diffusion(const double* states, double* ydot);
        void variable_step_ode_solve(double* RHS, double dt);
        void variable_step_hybrid_connections(const double* cvode_states_3d, double* const ydot_3d, const double* cvode_states_1d, double *const  ydot_1d);
        void scatter_grid_concentrations();
        void hybrid_connections();
        void set_diffusion(double*, int);
        void set_tortuosity(PyHocObject*);
        void set_volume_fraction(PyHocObject*);
        void do_multicompartment_reactions(double*);
        void initialize_multicompartment_reaction();
        void clear_multicompartment_reaction();
        int add_multicompartment_reaction(int, int*, int);
        double* set_rxd_currents(int, int*, PyHocObject**);
};

typedef struct ECSAdiDirection{
    void (*ecs_dg_adi_dir)(ECS_Grid_node*, const double, const int, const int, double const * const, double* const, double* const);
    double* states_in;
    double* states_out;
    int line_size;
} ECSAdiDirection;

typedef struct ECSAdiGridData{
    int start, stop;
    double* state;
    ECS_Grid_node* g;
    int sizej;
    ECSAdiDirection* ecs_adi_dir;
    double* scratchpad;
} ECSAdiGridData;

class ICS_Grid_node : public Grid_node{
    public:
        ICS_Grid_node();
        ~ICS_Grid_node();
        //fractional volumes
        double* _ics_alphas;
        //stores the positive x,y, and z neighbors for each node. [node0_x, node0_y, node0_z, node1_x ...]
        long* _neighbors;

        /*Line definitions from Python. In pattern of [line_start_node, line_length, ...]
        Array is sorted from longest to shortest line */
        long* _sorted_x_lines;
        long* _sorted_y_lines;
        long* _sorted_z_lines;

        //Lengths of _sorted_lines arrays. Used to find thread start and stop indices
        long _x_lines_length;
        long _y_lines_length;
        long _z_lines_length;

        //maximum line length for scratchpad memory allocation
        long _line_length_max;

        //total number of nodes for this grid
        long _num_nodes;

        //indices for thread start and stop positions 
        long* _nodes_per_thread;
        
        //Data for DG-ADI
        struct ICSAdiGridData* ics_tasks;
        struct ICSAdiDirection* ics_adi_dir_x;
        struct ICSAdiDirection* ics_adi_dir_y;
        struct ICSAdiDirection* ics_adi_dir_z;
 
        ICS_Grid_node(PyHocObject*, long, long*, long*, long, long*, long,
                      long*, long, double*, double*, double, bool, double,
                      double*); 
        void divide_x_work(const int nthreads);
        void divide_y_work(const int nthreads);
        void divide_z_work(const int nthreads);
        void set_num_threads(const int n);
        void do_grid_currents(double*, double, int);
        void apply_node_flux3D(double dt, double* states); 
        void volume_setup();
        int dg_adi();
        void variable_step_diffusion(const double* states, double* ydot);
        void variable_step_ode_solve(double* RHS, double dt);
        void hybrid_connections();
        void variable_step_hybrid_connections(const double* cvode_states_3d, double* const ydot_3d, const double* cvode_states_1d, double *const  ydot_1d);
        void scatter_grid_concentrations();
        void run_threaded_ics_dg_adi(struct ICSAdiDirection*);
        void set_diffusion(double*, int);
        void do_multicompartment_reactions(double*);
        void initialize_multicompartment_reaction();
        void clear_multicompartment_reaction();
        int add_multicompartment_reaction(int, int*, int);
        double* set_rxd_currents(int, int*, PyHocObject**);
};

typedef struct ICSAdiDirection{
    void (*ics_dg_adi_dir)(ICS_Grid_node* g, int, int, int, double, double*, double*, double*, double*, double*, double*);
    double* states_in;
    double* states_out;
    double* deltas;
    long* ordered_line_defs;
    long* ordered_nodes;
    long* ordered_start_stop_indices; 
    long* line_start_stop_indices;
    double dc;
    double* dcgrid;
    double d;
}ICSAdiDirection;


typedef struct ICSAdiGridData{
    //Start and stop node indices for lines
    int line_start, line_stop;
    //Where to start in the ordered_nodes array
    int ordered_start;
    double* state;
    ICS_Grid_node* g;
    ICSAdiDirection* ics_adi_dir;
    double* scratchpad;
    double* RHS;
    double* l_diag;
    double* diag;
    double* u_diag;
    //double* deltas;
 }ICSAdiGridData;
 
/***** GLOBALS *******************************************************************/
extern double *dt_ptr;              // Universal ∆t
extern double *t_ptr;               // Universal t

// static int N = 100;                 // Number of grid_lists (size of Parallel_grids)
extern Grid_node *Parallel_grids[100];// Array of Grid_node * lists
/*********************************************************************************/


// Set the global ∆t
void make_dt_ptr(PyHocObject* my_dt_ptr);


// Create a single Grid_node 
/* Parameters:  Python object that includes array of double pointers,
                size of x, y, and z dimensions
                x, y, and z diffusion constants
                delta x, delta y, and delta z*/

// Free a single Grid_node "grid"
//void free_Grid(Grid_node *grid);

// Insert a Grid_node "new_Grid" into the list located at grid_list_index in Parallel_grids
extern "C" int ECS_insert(int grid_list_index, PyHocObject* my_states, int my_num_states_x,
    int my_num_states_y, int my_num_states_z, double my_dc_x, double my_dc_y,
    double my_dc_z, double my_dx, double my_dy, double my_dz, 
    PyHocObject* my_alpha, PyHocObject* my_permeability, int, double, double);

Grid_node *ICS_make_Grid(PyHocObject* my_states, long num_nodes, long* neighbors, 
                long* x_line_defs, long x_lines_length, long* y_line_defs, long y_lines_length, long* z_line_defs,
                long z_lines_length, double* dcs, double dx, bool is_diffusable, double atolscale, double* ics_alphas);

// Insert an  ICS_Grid_node "new_Grid" into the list located at grid_list_index in Parallel_grids
extern "C" int ICS_insert(int grid_list_index, PyHocObject* my_states, long num_nodes, long* neighbors,
                long* x_line_defs, long x_lines_length, long* y_line_defs, long y_lines_length, long* z_line_defs,
                long z_lines_length, double* dcs, double dx, bool is_diffusable, double atolscale, double* ics_alphas);

extern "C" int ICS_insert_inhom(int grid_list_index, PyHocObject* my_states, long num_nodes, long* neighbors,
                long* x_line_defs, long x_lines_length, long* y_line_defs, long y_lines_length, long* z_line_defs,
                long z_lines_length, double* dcs, double dx, bool is_diffusable, double atolscale, double* ics_alphas);



// Set the diffusion coefficients for a given grid_id 
extern "C" int set_diffusion(int, int, double*,  int);

// Delete a specific Grid_node "find" from the list "head"
int remove(Grid_node **head, Grid_node *find);

// Destroy the list located at list_index and free all memory
void empty_list(int list_index);

void apply_node_flux(int, long*, double*, PyObject**, double, double*);