NMODL SympySolver - derivimplicit

This notebook describes the implementation of the derivimplicit part of the SympySolverVisitor, which solves the systems of ODEs defined in DERIVATIVE blocks when these ODEs are nonlinear.

For a higher level overview of the approach to solving ODEs in NMODL, please see the nmodl-odes-overview notebook.

For a more general tutorial on using the NMODL python interface, please see the tutorial notebook.

---

Implementation

The SympySolverVisitor for solver method derivimplicit does the following:

• Construct the implicit Euler equations to convert the sysytem of ODEs to a NONLINEAR block

• This NONLINEAR block is then solved as described in nmodl-nonlinear-solver

---

Implementation Tests

The unit tests may be helpful to understand what these functions are doing - SympySolverVisitor tests are located in test/visitor/sympy_solver.cpp, and tests involving the derivimplicit solver method have the tag “[derivimplicit]”

Examples

%%capture
! pip install nmodl

import nmodl.dsl as nmodl


def run_sympy_solver(mod_string, cse=False):
    # parse NMDOL file (supplied as a string) into AST
    driver = nmodl.NmodlDriver()
    AST = driver.parse_string(mod_string)
    # run SymtabVisitor to generate Symbol Table
    nmodl.symtab.SymtabVisitor().visit_program(AST)
    # constant folding, inlining & local variable renaming passes
    nmodl.visitor.ConstantFolderVisitor().visit_program(AST)
    nmodl.visitor.InlineVisitor().visit_program(AST)
    nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)
    # run SympySolver visitor
    nmodl.visitor.SympySolverVisitor().visit_program(AST)
    # return new DERIVATIVE block
    return nmodl.to_nmodl(
        nmodl.visitor.AstLookupVisitor().lookup(
            AST, nmodl.ast.AstNodeType.DERIVATIVE_BLOCK
        )[0]
    )

Ex. 1

2 coupled linear ODEs

ex1 = """
BREAKPOINT {
    SOLVE states METHOD sparse
}
STATE {
    mc m
}
DERIVATIVE states {
    mc' = -a*mc + b*m
    m' = a*mc - b*m
}
"""
print(run_sympy_solver(ex1, cse=True))

DERIVATIVE states {
    EIGEN_NEWTON_SOLVE[2]{
        LOCAL old_mc, old_m
    }{
        old_mc = mc
        old_m = m
    }{
        nmodl_eigen_x[0] = mc
        nmodl_eigen_x[1] = m
    }{
        nmodl_eigen_f[0] = -nmodl_eigen_x[0]*a*dt-nmodl_eigen_x[0]+nmodl_eigen_x[1]*b*dt+old_mc
        nmodl_eigen_j[0] = -a*dt-1.0
        nmodl_eigen_j[2] = b*dt
        nmodl_eigen_f[1] = nmodl_eigen_x[0]*a*dt-nmodl_eigen_x[1]*b*dt-nmodl_eigen_x[1]+old_m
        nmodl_eigen_j[1] = a*dt
        nmodl_eigen_j[3] = -b*dt-1.0
    }{
        mc = nmodl_eigen_x[0]
        m = nmodl_eigen_x[1]
    }{
    }
}

Ex. 2

5 coupled linear ODEs

ex2 = """
STATE {
    c1 o1 o2 p0 p1
}
BREAKPOINT  {
    SOLVE ihkin METHOD sparse
}
DERIVATIVE ihkin {
    LOCAL alpha, beta, k3p, k4, k1ca, k2
    evaluate_fct(v, cai)
    c1' = (-1*(alpha*c1-beta*o1))
    o1' = (1*(alpha*c1-beta*o1))+(-1*(k3p*o1-k4*o2))
    o2' = (1*(k3p*o1-k4*o2))
    p0' = (-1*(k1ca*p0-k2*p1))
    p1' = (1*(k1ca*p0-k2*p1))
}
"""
print(run_sympy_solver(ex2))

DERIVATIVE ihkin {
    EIGEN_NEWTON_SOLVE[5]{
        LOCAL alpha, beta, k3p, k4, k1ca, k2, old_c1, old_o1, old_o2, old_p0, old_p1
    }{
        evaluate_fct(v, cai)
        old_c1 = c1
        old_o1 = o1
        old_o2 = o2
        old_p0 = p0
        old_p1 = p1
    }{
        nmodl_eigen_x[0] = c1
        nmodl_eigen_x[1] = o1
        nmodl_eigen_x[2] = o2
        nmodl_eigen_x[3] = p0
        nmodl_eigen_x[4] = p1
    }{
        nmodl_eigen_f[0] = -nmodl_eigen_x[0]*alpha*dt-nmodl_eigen_x[0]+nmodl_eigen_x[1]*beta*dt+old_c1
        nmodl_eigen_j[0] = -alpha*dt-1.0
        nmodl_eigen_j[5] = beta*dt
        nmodl_eigen_j[10] = 0
        nmodl_eigen_j[15] = 0
        nmodl_eigen_j[20] = 0
        nmodl_eigen_f[1] = nmodl_eigen_x[0]*alpha*dt-nmodl_eigen_x[1]*beta*dt-nmodl_eigen_x[1]*dt*k3p-nmodl_eigen_x[1]+nmodl_eigen_x[2]*dt*k4+old_o1
        nmodl_eigen_j[1] = alpha*dt
        nmodl_eigen_j[6] = -beta*dt-dt*k3p-1.0
        nmodl_eigen_j[11] = dt*k4
        nmodl_eigen_j[16] = 0
        nmodl_eigen_j[21] = 0
        nmodl_eigen_f[2] = nmodl_eigen_x[1]*dt*k3p-nmodl_eigen_x[2]*dt*k4-nmodl_eigen_x[2]+old_o2
        nmodl_eigen_j[2] = 0
        nmodl_eigen_j[7] = dt*k3p
        nmodl_eigen_j[12] = -dt*k4-1.0
        nmodl_eigen_j[17] = 0
        nmodl_eigen_j[22] = 0
        nmodl_eigen_f[3] = -nmodl_eigen_x[3]*dt*k1ca-nmodl_eigen_x[3]+nmodl_eigen_x[4]*dt*k2+old_p0
        nmodl_eigen_j[3] = 0
        nmodl_eigen_j[8] = 0
        nmodl_eigen_j[13] = 0
        nmodl_eigen_j[18] = -dt*k1ca-1.0
        nmodl_eigen_j[23] = dt*k2
        nmodl_eigen_f[4] = nmodl_eigen_x[3]*dt*k1ca-nmodl_eigen_x[4]*dt*k2-nmodl_eigen_x[4]+old_p1
        nmodl_eigen_j[4] = 0
        nmodl_eigen_j[9] = 0
        nmodl_eigen_j[14] = 0
        nmodl_eigen_j[19] = dt*k1ca
        nmodl_eigen_j[24] = -dt*k2-1.0
    }{
        c1 = nmodl_eigen_x[0]
        o1 = nmodl_eigen_x[1]
        o2 = nmodl_eigen_x[2]
        p0 = nmodl_eigen_x[3]
        p1 = nmodl_eigen_x[4]
    }{
    }
}