ldifus.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include "section.h"
#include "membfunc.h"
#include "neuron.h"
#include "nrniv_mf.h"
#include "parse.hpp"
 
 
#define nt_t  nrn_threads->_t
#define nt_dt nrn_threads->_dt
 
extern "C" int diam_change_cnt;
extern "C" int structure_change_cnt;
 
typedef struct LongDifus {
    int dchange;
    int* mindex; /* index into memb_list[m] */
    int* pindex; /* parent in this struct */
    double** state;
    double* a;
    double* b;
    double* d;
    double* rhs;
    double* af;  /* efficiency benefit is very low, rall/dx^2 */
    double* bf;  /* 1/dx^2 */
    double* vol; /* volatile volume from COMPARTMENT */
    double* dc;  /* volatile diffusion constant * cross sectional
                area from LONGITUDINAL_DIFFUSION */
} LongDifus;
 
typedef struct LongDifusThreadData {
    int nthread;
    LongDifus** ldifus;
    Memb_list** ml;
} LongDifusThreadData;
 
static int ldifusfunccnt;
static ldifusfunc_t* ldifusfunc;
 
static ldifusfunc2_t stagger, ode, matsol, overall_setup;
 
void hoc_register_ldifus1(ldifusfunc_t f) {
    ldifusfunc = (ldifusfunc_t*) erealloc(ldifusfunc, (ldifusfunccnt + 1) * sizeof(ldifusfunc_t));
    ldifusfunc[ldifusfunccnt] = f;
    ++ldifusfunccnt;
}
 
#if MAC
/* this avoids a missing _ptrgl12 in the mac library that was called by
the MrC compiled object
*/
void mac_difusfunc(ldifusfunc2_t* f,
                   int m,
                   ldifusfunc3_t diffunc,
                   void** v,
                   int ai,
                   int sindex,
                   int dindex,
                   NrnThread* nt) {
    (*f)(m, diffunc, v, ai, sindex, dindex, nt);
}
#endif
 
 
extern "C" void nrn_tree_solve(double* a, double* d, double* b, double* rhs, int* pindex, int n) {
    /*
        treesolver
 
        a      - above the diagonal
        d      - diagonal
        b      - below the diagonal
        rhs    - right hand side, which is changed to the result
        pindex - parent indices
        n      - number of states
    */
 
    int i;
 
    /* triang */
    for (i = n - 1; i > 0; --i) {
        int pin = pindex[i];
        if (pin > -1) {
            double p;
            p = a[i] / d[i];
            d[pin] -= p * b[i];
            rhs[pin] -= p * rhs[i];
        }
    }
    /* bksub */
    for (i = 0; i < n; ++i) {
        int pin = pindex[i];
        if (pin > -1) {
            rhs[i] -= b[i] * rhs[pin];
        }
        rhs[i] /= d[i];
    }
}
 
 
void long_difus_solve(int method, NrnThread* nt) {
    ldifusfunc2_t* f = NULL;
    int i;
    if (ldifusfunc) {
        switch (method) {
        case 0: /* normal staggered time step */
            f = stagger;
            break;
        case 1: /* dstate = f(state) */
            f = ode;
            break;
        case 2: /* solve (1 + dt*jacobian)*x = b */
            f = matsol;
            break;
        case 3: /* setup only called by thread 0 */
            f = overall_setup;
            break;
        }
        assert(f);
 
        for (i = 0; i < ldifusfunccnt; ++i) {
            (*ldifusfunc[i])(f, nt);
        }
    }
}
 
static void longdifusfree(LongDifus** ppld) {
    if (*ppld) {
        LongDifus* pld = *ppld;
#if 0
printf("free longdifus structure_change=%d %d\n", pld->schange, structure_change_cnt);
#endif
        free(pld->mindex);
        free(pld->pindex);
        free(pld->state);
        free(pld->a);
        free(pld->b);
        free(pld->d);
        free(pld->rhs);
        free(pld->af);
        free(pld->bf);
        free(pld->vol);
        free(pld->dc);
        free(pld);
    }
    *ppld = (LongDifus*) 0;
}
 
static void longdifus_diamchange(LongDifus* pld, int m, int sindex, Memb_list* ml, NrnThread* _nt) {
    int i, n, mi, mpi, j, index, pindex, vnodecount;
    Node *nd, *pnd;
    double rall, dxp, dxc;
 
    if (pld->dchange == diam_change_cnt) {
        return;
    }
    /*printf("longdifus_diamchange %d %d\n", pld->dchange, diam_change_cnt);*/
    vnodecount = _nt->end;
    n = ml->nodecount;
 
    for (i = 0; i < n; ++i) {
        /* For every child with a parent having this mechanism */
        /* Also child may butte end to end with parent or attach to middle */
        mi = pld->mindex[i];
        if (sindex < 0) {
            pld->state[i] = ml->pdata[mi][-sindex - 1].pval;
        } else {
            pld->state[i] = ml->data[mi] + sindex;
        }
        nd = ml->nodelist[mi];
        pindex = pld->pindex[i];
        if (pindex > -1) {
            mpi = pld->mindex[pindex];
            pnd = ml->nodelist[mpi];
            if (nd->sec_node_index_ == 0) {
                rall = nd->sec->prop->dparam[4].val;
            } else {
                rall = 1.;
            }
            dxc = section_length(nd->sec) / ((double) (nd->sec->nnode - 1));
            dxp = section_length(pnd->sec) / ((double) (pnd->sec->nnode - 1));
            pld->af[i] = 2 * rall / dxp / (dxc + dxp);
            pld->bf[i] = 2 / dxc / (dxc + dxp);
        }
    }
    pld->dchange = diam_change_cnt;
}
 
static void longdifusalloc(LongDifus** ppld, int m, int sindex, Memb_list* ml, NrnThread* _nt) {
    LongDifus* pld;
    int i, n, mi, mpi, j, index, pindex, vnodecount;
    int *map, *omap;
    Node *nd, *pnd;
    hoc_Item* qsec;
 
    vnodecount = _nt->end;
    *ppld = pld = (LongDifus*) emalloc(sizeof(LongDifus));
    n = ml->nodecount;
 
    pld->dchange = 0;
    pld->mindex = (int*) ecalloc(n, sizeof(int));
    pld->pindex = (int*) ecalloc(n, sizeof(int));
    pld->state = (double**) ecalloc(n, sizeof(double*));
    pld->a = (double*) ecalloc(n, sizeof(double));
    pld->b = (double*) ecalloc(n, sizeof(double));
    pld->d = (double*) ecalloc(n, sizeof(double));
    pld->rhs = (double*) ecalloc(n, sizeof(double));
    pld->af = (double*) ecalloc(n, sizeof(double));
    pld->bf = (double*) ecalloc(n, sizeof(double));
    pld->vol = (double*) ecalloc(n, sizeof(double));
    pld->dc = (double*) ecalloc(n, sizeof(double));
 
    /* make a map from node_index to memb_list index. -1 means no exist*/
    map = (int*) ecalloc(vnodecount, sizeof(int));
    omap = (int*) ecalloc(n, sizeof(int));
    for (i = 0; i < vnodecount; ++i) {
        map[i] = -1;
    }
    for (i = 0; i < n; ++i) {
        map[ml->nodelist[i]->v_node_index] = i;
    }
#if 0
for (i=0; i < vnodecount; ++i) {
    printf("%d index=%d\n", i, map[i]);
}
#endif
    /* order the indices for efficient gaussian elimination */
    /* But watch out for 0 area nodes. Use the parent of parent */
    /* But if parent of parent does not have diffusion mechanism
       check the parent section */
    /* And watch out for root. Use first node of root section */
    for (i = 0, j = 0; i < vnodecount; ++i) {
        if (map[i] > -1) {
            pld->mindex[j] = map[i];
            omap[map[i]] = j; /* from memb list index to order */
            pnd = _nt->_v_parent[i];
            pindex = map[pnd->v_node_index];
            if (pindex == -1) { /* maybe this was zero area node */
                pnd = _nt->_v_parent[pnd->v_node_index];
                if (pnd) {
                    pindex = map[pnd->v_node_index];
                    if (pindex < 0) { /* but what about the parent section */
                        Section* psec = _nt->_v_node[i]->sec->parentsec;
                        if (psec) {
                            pnd = psec->pnode[0];
                            pindex = map[pnd->v_node_index];
                        }
                    }
                } else { /* maybe this section is not the root */
                    Section* psec = _nt->_v_node[i]->sec->parentsec;
                    if (psec) {
                        pnd = psec->pnode[0];
                        pindex = map[pnd->v_node_index];
                    }
                }
            }
            if (pindex > -1) {
                pld->pindex[j] = omap[pindex];
            } else {
                pld->pindex[j] = -1;
            }
            ++j;
        }
    }
    longdifus_diamchange(pld, m, sindex, ml, _nt);
#if 0
    for (i=0; i < n; ++i) {
printf("i=%d pin=%d mi=%d :%s node %d state[(%i)]=%g\n", i, pld->pindex[i],
    pld->mindex[i], secname(ml->nodelist[pld->mindex[i]]->sec),
    ml->nodelist[pld->mindex[i]]->sec_node_index_
    , sindex, pld->state[i][0]);
    }
#endif
    free(map);
    free(omap);
}
 
/* called at end of v_setup_vectors only for thread 0 */
/* only v makes sense and the purpose is to free the old, allocate space */
/* for the new, and setup the tml field */
/* the args used are m and v */
static void overall_setup(int m,
                          ldifusfunc3_t diffunc,
                          void** v,
                          int ai,
                          int sindex,
                          int dindex,
                          NrnThread* _nt) {
    int i;
    LongDifusThreadData** ppldtd = (LongDifusThreadData**) v;
    LongDifusThreadData* ldtd = *ppldtd;
    if (ldtd) { /* free the whole thing */
        free(ldtd->ml);
        for (i = 0; i < ldtd->nthread; ++i) {
            if (ldtd->ldifus[i]) {
                longdifusfree(ldtd->ldifus + i);
            }
        }
        free(ldtd->ldifus);
        free(ldtd);
        *ppldtd = (LongDifusThreadData*) 0;
    }
    /* new overall space */
    *ppldtd = (LongDifusThreadData*) emalloc(sizeof(LongDifusThreadData));
    ldtd = *ppldtd;
    ldtd->nthread = nrn_nthread;
    ldtd->ldifus = (LongDifus**) ecalloc(nrn_nthread, sizeof(LongDifus*));
    ldtd->ml = (Memb_list**) ecalloc(nrn_nthread, sizeof(Memb_list*));
    /* which have memb_list pointers */
    for (i = 0; i < nrn_nthread; ++i) {
        NrnThreadMembList* tml;
        for (tml = nrn_threads[i].tml; tml; tml = tml->next) {
            if (tml->index == m) {
                ldtd->ml[i] = tml->ml;
                longdifusalloc(ldtd->ldifus + i, m, sindex, tml->ml, nrn_threads + i);
                break;
            }
        }
    }
}
 
static LongDifus* v2ld(void** v, int tid) {
    LongDifusThreadData** ppldtd = (LongDifusThreadData**) v;
    return (*ppldtd)->ldifus[tid];
}
 
static Memb_list* v2ml(void** v, int tid) {
    LongDifusThreadData** ppldtd = (LongDifusThreadData**) v;
    return (*ppldtd)->ml[tid];
}
 
static void
stagger(int m, ldifusfunc3_t diffunc, void** v, int ai, int sindex, int dindex, NrnThread* _nt) {
    LongDifus* pld;
    int i, n, di;
    double dc, vol, dfdi, dx;
    double** data;
    Datum** pdata;
    Datum* thread;
    Memb_list* ml;
 
    di = dindex + ai;
 
    pld = v2ld(v, _nt->id);
    if (!pld)
        return;
    ml = v2ml(v, _nt->id);
 
    n = ml->nodecount;
    data = ml->data;
    pdata = ml->pdata;
    thread = ml->_thread;
 
    longdifus_diamchange(pld, m, sindex, ml, _nt);
    /*flux and volume coefficients (if dc is constant this is too often)*/
    for (i = 0; i < n; ++i) {
        int pin = pld->pindex[i];
        int mi = pld->mindex[i];
        pld->dc[i] = (*diffunc)(ai, data[mi], pdata[mi], pld->vol + i, &dfdi, thread, _nt);
        pld->d[i] = 0.;
#if 0
        if (dfdi) {
            pld->d[i] += fabs(dfdi)/pld->vol[i]/pld->state[i][ai];
        }
#endif
        if (pin > -1) {
            /* D * area between compartments */
            dc = (pld->dc[i] + pld->dc[pin]) / 2.;
 
            pld->a[i] = -pld->af[i] * dc / pld->vol[pin];
            pld->b[i] = -pld->bf[i] * dc / pld->vol[i];
        }
    }
    /* setup matrix */
    for (i = 0; i < n; ++i) {
        int pin = pld->pindex[i];
        int mi = pld->mindex[i];
        pld->d[i] += 1. / nt_dt;
        pld->rhs[i] = pld->state[i][ai] / nt_dt;
        if (pin > -1) {
            pld->d[i] -= pld->b[i];
            pld->d[pin] -= pld->a[i];
        }
    }
#if 0
for (i=0; i < n; ++i) { double a,b;
    if (pld->pindex[i] > -1) {
        a = pld->a[i];
        b = pld->b[i];
    }else{ a=b=0.;}
    printf("i=%d a=%g b=%g d=%g rhs=%g state=%g\n",
     i, a, b, pld->d[i], pld->rhs[i], pld->state[i][ai]);
}
#endif
 
    /* we've set up the matrix; now solve it */
    nrn_tree_solve(pld->a, pld->d, pld->b, pld->rhs, pld->pindex, n);
 
    /* update answer */
    for (i = 0; i < n; ++i) {
        pld->state[i][ai] = pld->rhs[i];
    }
}
 
static void
ode(int m, ldifusfunc3_t diffunc, void** v, int ai, int sindex, int dindex, NrnThread* _nt) {
    LongDifus* pld;
    int i, n, di;
    double dc, vol, dfdi;
    double** data;
    Datum** pdata;
    Datum* thread;
    Memb_list* ml;
 
    di = dindex + ai;
 
    pld = v2ld(v, _nt->id);
    if (!pld)
        return;
    ml = v2ml(v, _nt->id);
 
    n = ml->nodecount;
    data = ml->data;
    pdata = ml->pdata;
    thread = ml->_thread;
 
    longdifus_diamchange(pld, m, sindex, ml, _nt);
    /*flux and volume coefficients (if dc is constant this is too often)*/
    for (i = 0; i < n; ++i) {
        int pin = pld->pindex[i];
        int mi = pld->mindex[i];
        pld->dc[i] = (*diffunc)(ai, data[mi], pdata[mi], pld->vol + i, &dfdi, thread, _nt);
        if (pin > -1) {
            /* D * area between compartments */
            dc = (pld->dc[i] + pld->dc[pin]) / 2.;
 
            pld->a[i] = pld->af[i] * dc / pld->vol[pin];
            pld->b[i] = pld->bf[i] * dc / pld->vol[i];
        }
    }
    /* add terms to diffeq */
    for (i = 0; i < n; ++i) {
        double dif;
        int pin = pld->pindex[i];
        int mi = pld->mindex[i];
#if 0
        pld->d[i] = data[mi][di];
#endif
        if (pin > -1) {
            dif = (pld->state[pin][ai] - pld->state[i][ai]);
            data[mi][di] += dif * pld->b[i];
            data[pld->mindex[pin]][di] -= dif * pld->a[i];
        }
    }
#if 0
    for (i=0; i < n; ++i) {
        int mi = pld->mindex[i];
        printf("%d olddstate=%g new=%g\n", i, pld->d[i], data[mi][di]);
    }
#endif
}
 
 
static void
matsol(int m, ldifusfunc3_t diffunc, void** v, int ai, int sindex, int dindex, NrnThread* _nt) {
    LongDifus* pld;
    int i, n, di;
    double dc, vol, dfdi;
    double** data;
    Datum** pdata;
    Datum* thread;
    Memb_list* ml;
 
    di = dindex + ai;
 
    pld = v2ld(v, _nt->id);
    if (!pld)
        return;
    ml = v2ml(v, _nt->id);
 
    n = ml->nodecount;
    data = ml->data;
    pdata = ml->pdata;
    thread = ml->_thread;
 
    /*flux and volume coefficients (if dc is constant this is too often)*/
    for (i = 0; i < n; ++i) {
        int pin = pld->pindex[i];
        int mi = pld->mindex[i];
        pld->dc[i] = (*diffunc)(ai, data[mi], pdata[mi], pld->vol + i, &dfdi, thread, _nt);
        pld->d[i] = 0.;
        if (dfdi) {
            pld->d[i] += fabs(dfdi) / pld->vol[i] / pld->state[i][ai];
#if 0
printf("i=%d state=%g vol=%g dfdc=%g\n", i, pld->state[i][ai],pld->vol[i], pld->d[i]);
#endif
        }
        if (pin > -1) {
            /* D * area between compartments */
            dc = (pld->dc[i] + pld->dc[pin]) / 2.;
 
            pld->a[i] = -pld->af[i] * dc / pld->vol[pin];
            pld->b[i] = -pld->bf[i] * dc / pld->vol[i];
        }
    }
    /* setup matrix */
    for (i = 0; i < n; ++i) {
        int pin = pld->pindex[i];
        int mi = pld->mindex[i];
        pld->d[i] += 1. / nt_dt;
        pld->rhs[i] = data[mi][di] / nt_dt;
        if (pin > -1) {
            pld->d[i] -= pld->b[i];
            pld->d[pin] -= pld->a[i];
        }
    }
#if 0
for (i=0; i < n; ++i) { double a,b;
    int mi = pld->mindex[i];
    if (pld->pindex[i] > -1) {
        a = pld->a[i];
        b = pld->b[i];
    }else{ a=b=0.;}
    printf("i=%d a=%g b=%g d=%g rhs=%g dstate=%g\n",
     i, a, b, pld->d[i], pld->rhs[i], data[mi][di]);
}
#endif
    /* triang */
    for (i = n - 1; i > 0; --i) {
        int pin = pld->pindex[i];
        if (pin > -1) {
            double p;
            p = pld->a[i] / pld->d[i];
            pld->d[pin] -= p * pld->b[i];
            pld->rhs[pin] -= p * pld->rhs[i];
        }
    }
    /* bksub */
    for (i = 0; i < n; ++i) {
        int pin = pld->pindex[i];
        if (pin > -1) {
            pld->rhs[i] -= pld->b[i] * pld->rhs[pin];
        }
        pld->rhs[i] /= pld->d[i];
    }
    /* update answer */
    for (i = 0; i < n; ++i) {
        int mi = pld->mindex[i];
        data[mi][di] = pld->rhs[i];
    }
}