nrndae.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <cstdio>
#include "nrndae.h"
#include "nrnoc2iv.h"
 
extern void nrndae_alloc();
extern int nrndae_extra_eqn_count();
extern void nrndae_init();
extern void nrndae_rhs();  // relative to c*dy/dt = -g*y + b
extern void nrndae_lhs();
extern void nrndae_dkmap(double**, double**);
extern void nrndae_dkres(double*, double*, double*);
extern void nrndae_dkpsol(double);
extern void nrndae_update();
extern void nrn_matrix_node_free();
extern int cvode_active_;
extern int nrn_use_daspk_;
extern int secondorder;
extern int nrndae_list_is_empty();
 
static NrnDAEPtrList nrndae_list;
 
int nrndae_list_is_empty() {
    return nrndae_list.empty() ? 1 : 0;
}
 
 
void nrndae_register(NrnDAE* n) {
    nrndae_list.push_back(n);
}
 
void nrndae_deregister(NrnDAE* n) {
    nrndae_list.remove(n);
}
 
int nrndae_extra_eqn_count() {
    int neqn = 0;
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        neqn += (*m)->extra_eqn_count();
    }
    return neqn;
}
 
void nrndae_update() {
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->update();
    }
}
 
void nrndae_alloc() {
    NrnThread* _nt = nrn_threads;
    nrn_thread_error("NrnDAE only one thread allowed");
    int neqn = _nt->end;
    if (_nt->_ecell_memb_list) {
        neqn += _nt->_ecell_memb_list->nodecount * nlayer;
    }
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->alloc(neqn + 1);
        neqn += (*m)->extra_eqn_count();
    }
}
 
 
void nrndae_init() {
    if ((!nrndae_list.empty()) &&
        (secondorder > 0 || ((cvode_active_ > 0) && (nrn_use_daspk_ == 0)))) {
        hoc_execerror("NrnDAEs only work with secondorder==0 or daspk", 0);
    }
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->init();
    }
}
 
void nrndae_rhs() {
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->rhs();
    }
}
 
void nrndae_lhs() {
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->lhs();
    }
}
 
void nrndae_dkmap(double** pv, double** pvdot) {
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->dkmap(pv, pvdot);
    }
}
 
void nrndae_dkres(double* y, double* yprime, double* delta) {
    // c*y' = f(y) so
    // delta = c*y' - f(y)
    for (NrnDAEPtrListIterator m = nrndae_list.begin(); m != nrndae_list.end(); m++) {
        (*m)->dkres(y, yprime, delta);
    }
}
 
inline void NrnDAE::alloc_(int size, int start, int nnode, Node** nodes, int* elayer) {}
 
void NrnDAE::alloc(int start_index) {
    // printf("NrnDAE::alloc %lx\n", (long)this);
    size_ = y_.size();
    if (y0_) {
        assert(y0_->size() == size_);
    }
    assert(c_->nrow() == size_ && c_->ncol() == size_);
    cyp_.resize(size_);
    yptmp_.resize(size_);
    start_ = start_index;
    // printf("start=%d size=%d\n", start_, size_);
    delete[] bmap_;
    bmap_ = new int[size_];
    for (int i = 0; i < size_; ++i) {
        if (i < nnode_) {
            bmap_[i] = nodes_[i]->eqn_index_ + elayer_[i];
            if (elayer_[i] > 0 && !nodes_[i]->extnode) {
                // hoc_execerror(secname(nodes_[i]->sec), "NrnDAE: Referring to an extracellular
                // layer but\nextracellular is not inserted.");
                // instead treat as though connected to ground.
                bmap_[i] = 0;
            }
        } else {
            bmap_[i] = start_ + i - nnode_;
        }
    }
    // printf("c_->alloc start=%d, nnode=%d\n", start_, nnode_);
    c_->alloc(start_, nnode_, nodes_, elayer_);
 
    // allow subclasses to do their own allocations as well
    alloc_(size_, start_, nnode_, nodes_, elayer_);
}
 
NrnDAE::NrnDAE(Matrix* cmat,
               Vect* const yvec,
               Vect* const y0,
               int nnode,
               Node** const nodes,
               Vect* const elayer,
               void (*f_init)(void* data),
               void* const data)
    : y_(*yvec)
    , yptmp_((Object*) NULL)
    , cyp_((Object*) NULL)
    , f_init_(f_init)
    , data_(data) {
    // printf("NrnDAE %lx\n", (long)this);
    if (cmat) {
        assumed_identity_ = NULL;
    } else {
        const int size = y_.size();
        assumed_identity_ = new OcSparseMatrix(size, size);
        // assumed_identity_->setdiag(0, 1);
        for (int i = 0; i < size; i++)
            (*assumed_identity_)(i, i) = 1;
        cmat = assumed_identity_;
    }
    c_ = new MatrixMap(cmat);
    Vect& elay = *elayer;
    nnode_ = nnode;
    nodes_ = nodes;
    if (nnode_ > 0) {
        elayer_ = new int[nnode_];
        if (elayer) {
            for (int i = 0; i < nnode_; ++i) {
                elayer_[i] = int(elay[i]);
            }
        } else {
            for (int i = 0; i < nnode_; ++i) {
                elayer_[i] = 0;
            }
        }
    } else {
        elayer_ = NULL;
    }
    y0_ = y0;
    bmap_ = new int[1];
 
    nrndae_register(this);
    //  use_sparse13 = 1;
    nrn_matrix_node_free();
}
 
 
NrnDAE::~NrnDAE() {
    nrndae_deregister(this);
    delete[] bmap_;
    delete c_;
    delete assumed_identity_;
    if (elayer_) {
        delete[] elayer_;
    }
    //  if (nrndae_list->count() == 0) {
    //      use_sparse13 = 0;
    //  }
    nrn_matrix_node_free();
}
 
 
int NrnDAE::extra_eqn_count() {
    // printf("NrnDAE::extra_eqn_count %lx\n", (long)this);
    // printf("  nnode_=%d g_->nrow()=%d\n", nnode_, g_->nrow());
    return c_->nrow() - nnode_;
}
 
void NrnDAE::dkmap(double** pv, double** pvdot) {
    // printf("NrnDAE::dkmap\n");
    NrnThread* _nt = nrn_threads;
    for (int i = nnode_; i < size_; ++i) {
        // printf("bmap_[%d] = %d\n", i, bmap_[i]);
        pv[bmap_[i] - 1] = y_.data() + i;
        pvdot[bmap_[i] - 1] = _nt->_actual_rhs + bmap_[i];
    }
}
 
void NrnDAE::update() {
    // printf("NrnDAE::update %lx\n", (long)this);
    NrnThread* _nt = nrn_threads;
    // note that the following is correct also for states that refer
    // to the internal potential of a segment. i.e rhs is v + vext[0]
    for (int i = 0; i < size_; ++i) {
        y_[i] += _nt->_actual_rhs[bmap_[i]];
    }
    // for (int i=0; i < size_; ++i) printf(" i=%d bmap_[i]=%d y_[i]=%g\n", i, bmap_[i],
    // y_->elem(i));
}
 
void NrnDAE::init() {
    // printf("NrnDAE::init %lx\n", (long)this);
    // printf("init size_=%d %d %d %d\n", size_, y_->size(), y0_->size(), b_->size());
    Vect& y0 = *y0_;
 
    v2y();
    if (f_init_) {
        f_init_(data_);
    } else {
        if (y0_) {
            for (int i = nnode_; i < size_; ++i) {
                y_[i] = y0[i];
            }
        } else {
            for (int i = nnode_; i < size_; ++i) {
                y_[i] = 0.;
            }
        }
    }
    // for (i=0; i < nnode_; ++i) printf(" i=%d y[i]=%g\n", i, y[i]);
}
 
void NrnDAE::v2y() {
    // vm,vext may be reinitialized between fixed steps and certainly
    // have been adjusted by daspk
 
    for (int i = 0; i < nnode_; ++i) {
        Node* nd = nodes_[i];
        // Note elayer_[0] refers to internal.
        if (elayer_[i] == 0) {
            y_[i] = NODEV(nd);
            if (nd->extnode) {
                y_[i] += nd->extnode->v[0];
            }
        } else if (nd->extnode) {
            y_[i] = nd->extnode->v[elayer_[i] - 1];
        }
    }
}
 
 
void NrnDAE::dkres(double* y, double* yprime, double* delta) {
    // printf("NrnDAE::dkres %lx\n", (long)this);
    // delta is already f(y)
    // now subtract c*y'
    // The problem is the map between y and yprime and the local
    // representation of the y and yprime
 
    for (int i = 0; i < size_; ++i) {
        // printf("%d %d %g %g\n", i, bmap_[i]-1, y[bmap_[i]-1], yprime[bmap_[i]-1]);
        yptmp_[i] = yprime[bmap_[i] - 1];
    }
    Vect* cyp;
    if (assumed_identity_) {
        // if c_ is KNOWN to be the identity matrix, then no multiplication to do
        // for now, this happens only when cmat = NULL
        // note that the user might change cmat, so we can't assume that if it's initially
        // the identity that it will stay that way
        cyp = &yptmp_;
    } else {
        c_->mulv(yptmp_, cyp_);  // mulv cannot multiply in place
        cyp = &cyp_;
    }
    for (int i = 0; i < size_; ++i) {
        delta[bmap_[i] - 1] -= (*cyp)[i];
    }
}
 
 
void NrnDAE::rhs() {
    NrnThread* _nt = nrn_threads;
    v2y();
    f_(y_, yptmp_, size_);
    for (int i = 0; i < size_; ++i) {
        _nt->_actual_rhs[bmap_[i]] += yptmp_[i];
    }
}
 
void NrnDAE::lhs() {
    // printf("NrnDAE::lhs %lx\n", (long)this);
    // printf("  nrn_threads[0].cj = %g\n", nrn_threads[0].cj);
    // left side portion of (c/dt - J)*[dy] =  f(y)
    c_->add(nrn_threads[0].cj);
    v2y();
    jacobian_(y_)->add(jacobian_multiplier_() * -1);
}