hugogranstrom / nimsundials Goto Github PK

import os

import nimterop/[cimport, paths,git]

const
  version = "v5.0.0-dev.1"
  srcDir = currentSourcePath().parentDir()/"sundials"
  buildDir= srcDir/"build"
  installDir = "lib"
  installDirLib64 = buildDir/installDir/"usr/local/lib64"
  installDirInclude = buildDir/installDir/"usr/local/include"

static:
  # https://github.com/JuliaDiffEq/Sundials.jl
  gitPull("https://github.com/LLNL/sundials", outdir = srcDir, checkout = "tags/v5.0.0-dev.1")
  cDebug()
  cDisableCaching()


  cmake(buildDir,"Makefile","../.") #" -DCMAKE_INSTALL_PREFIX=" & buildDir & " -DCMAKE_INSTALL_LIBDIR=" & installDir )
  make(buildDir,"", "")
  make(buildDir,"prueba.txt", "DESTDIR=\"" & buildDir & "/lib" & "\" install")


#[
https://github.com/LLNL/sundials/blob/master/INSTALL_GUIDE.pdf

In Archlinux: https://git.archlinux.org/svntogit/community.git/tree/trunk/PKGBUILD?h=packages/sundials
cmake ../$pkgname-$pkgver \
  -DCMAKE_INSTALL_PREFIX=/usr \
  -DCMAKE_INSTALL_LIBDIR=lib \
  -DMPI_ENABLE=ON \
  -DPTHREAD_ENABLE=ON	\
  -DOPENMP_ENABLE=ON \
  -DF77_INTERFACE_ENABLE=ON \
  -DKLU_ENABLE=ON \
  -DKLU_LIBRARY_DIR=/usr/lib \
  -DEXAMPLES_INSTALL_PATH=/usr/share/sundials/examples
]#

cIncludeDir(installDirInclude)
#cIncludeDir(buildDir/"lib")
cIncludeDir(installDirLib64)

cPlugin:
    import strutils
    proc onSymbol*(sym: var Symbol) {.exportc, dynlib.} =
        if sym.kind == nskType:
            if sym.name == "_N_VectorContent_Serial":
                sym.name = "N_VectorContent_Serial_Base"
        if sym.name == "ModifiedGS":
            sym.name = "cModifiedGS"
        if sym.name == "ClassicalGS":
            sym.name = "cClassicalGS"
        if sym.name == "_SpgmrMemRec":
            sym.name = "cSpgmrMemRec"
        if sym.name == "iS":
            sym.name = "ciS"
        sym.name = sym.name.strip(chars = {'_'})



const
    libsundials_nvecserial     = installDirLib64/"libsundials_nvecserial.so"
    libsundials_cvode          = installDirLib64/"libsundials_cvode.so"
    #libsundials_cvodes          = installDirLib64/"libsundials_cvodes.so"
    libsundials_sunlinsolspgmr = installDirLib64/"libsundials_sunlinsolspgmr.so"
    libsundials_sunmatrixdense = installDirLib64/"libsundials_sunmatrixdense.so"

cImport(installDirInclude/"sundials/sundials_config.h")
cImport(installDirInclude/"sundials/sundials_types.h")
cImport(installDirInclude/"sundials/sundials_math.h")

cImport(installDirInclude/"sundials/sundials_nvector.h")
#cImport(installDirInclude/"nvector/nvector_serial2.h")
cImport(installDirInclude/"nvector/nvector_serial.h", dynlib = "libsundials_nvecserial")

cImport(installDirInclude/"sundials/sundials_nonlinearsolver.h")
cImport(installDirInclude/"sundials/sundials_direct.h")
cImport(installDirInclude/"sundials/sundials_iterative.h")
cImport(installDirInclude/"sundials/sundials_matrix.h", dynlib="libsundials_cvode")
cImport(installDirInclude/"sundials/sundials_linearsolver.h", dynlib="libsundials_cvode")

#cImport(installDirInclude/"sundials/sundials_spgmr.h")
cImport(installDirInclude/"sunlinsol/sunlinsol_spgmr.h", dynlib="libsundials_sunlinsolspgmr")

cImport(installDirInclude/"sunmatrix/sunmatrix_dense.h", dynlib = "libsundials_sunmatrixdense")

cImport(installDirInclude/"cvode/cvode_ls.h", dynlib="libsundials_cvode")
#cImport(installDirInclude/"cvode/cvode2.h")
cImport(installDirInclude/"cvode/cvode.h", dynlib="libsundials_cvode")

# The following avoids importing twice the same definitions (already defined somewhere else).
#[
static:
   cSkipSymbol @["CV_ADAMS", "CV_BDF", "CV_NORMAL", "CV_ONE_STEP", "CV_SUCCESS", "CV_TSTOP_RETURN", "CV_ROOT_RETURN", "CV_WARNING",
   "CV_TOO_MUCH_WORK", "CV_TOO_MUCH_ACC", "CV_ERR_FAILURE", "CV_CONV_FAILURE", "CV_LINIT_FAIL", "CV_LSETUP_FAIL", "CV_LSOLVE_FAIL",
   "CV_RHSFUNC_FAIL", "CV_FIRST_RHSFUNC_ERR", "CV_REPTD_RHSFUNC_ERR", "CV_UNREC_RHSFUNC_ERR", "CV_RTFUNC_FAIL", "CV_NLS_INIT_FAIL",
   "CV_NLS_SETUP_FAIL", "CV_CONSTR_FAIL", "CV_MEM_FAIL", "CV_MEM_NULL", "CV_ILL_INPUT", "CV_NO_MALLOC", "CV_BAD_K", "CV_BAD_T",
   "CV_BAD_DKY", "CV_TOO_CLOSE", "CV_VECTOROP_ERR", "CVRhsFn", "CVRootFn", "CVEwtFn", "CVErrHandlerFn", "CVodeCreate"]

cImport(installDirInclude/"cvodes/cvodes.h", dynlib="libsundials_cvodes")
]#


template ptr2Array[T](p: pointer): auto = cast[ptr UncheckedArray[T]](p)
template array2Ptr[T](arr: openArray[T]): auto = addr(arr[0])

template `->`(a, b: untyped): untyped =
    a[].b

template NV_CONTENT_S(v: untyped): untyped = cast[N_VectorContent_Serial](v->content)
template NV_LENGTH_S(v: untyped): untyped = NV_CONTENT_S(v) -> length
template NV_OWN_DATA_S(v: untyped): untyped = NV_CONTENT_S(v) -> own_data
template NV_DATA_S(v: untyped): untyped = NV_CONTENT_S(v) -> data
template NV_Ith_S(v: untyped, i: sunindextype): untyped = ptr2Array[realtype](NV_DATA_S(v))[i]



type
    NVectorType* = object
        length*: int
        rawVector*: ref[N_Vector]

proc newNVector*(length: int): NVectorType =
    if length <= 0:
        raise newException(ValueError, "NVector length must be greater than 0")
    result.length = length
    result.rawVector = new N_Vector
    result.rawVector[] = N_VNew_Serial(result.length)

# create newNVector(length) and for over arr, NV_Ith_S(result, i) = arr[i] to keep it in memory.
proc newNVector*(arr: openArray[realtype]): NVectorType =
    if arr.len <= 0:
        raise newException(ValueError, "NVector length must be greater than 0")
    result.length = arr.len
    result.rawVector = new N_Vector
    result.rawVector[] = N_VNew_Serial(result.length)
    for i in 0 ..< result.length:
        NV_Ith_S(result.rawVector[], i) = arr[i]

proc newNVector*(v: N_Vector): NVectorType =
    result = newNVector(NV_LENGTH_S(v).int)
    N_VScale_Serial(1.0, v, result.rawVector[])

proc clone*(v: NVectorType): NVectorType =
    result = newNVector(v.length)
    N_VScale_Serial(1.0, v.rawVector[], result.rawVector[])


proc `[]`*(v: NVectorType, i: int): realtype =
    if v.length <= i:
        raise newException(ValueError, "index i is out of range. `[]`")
    NV_Ith_S(v.rawVector[], i)

proc `[]=`*(v: NVectorType, i: int, c: realtype) =
    if v.length <= i:
        raise newException(ValueError, "index i is out of range. `[]`")
    NV_Ith_S(v.rawVector[], i) = c

proc `$`*(v: NVectorType): string =
    result = "NVector("
    for i in 0 ..< v.length:
        result = result & $v[i] & ", "
    result = result & ")"

proc `@`*(v: NVectorType): seq[realtype] =
    for i in 0 ..< v.length:
        result.add(v[i])

proc `==`*(a, b: NVectorType): bool =
    let length = a.length
    if length != b.length:
        raise newException(ValueError, "NVectors must have same lengths for `==`")
    let aData = ptr2Array[realtype](NV_DATA_S(a.rawVector[]))
    let bData = ptr2Array[realtype](NV_DATA_S(b.rawVector[]))
    for i in 0 ..< length:
        if aData[i] != bData[i]:
            return false
    return true

proc `+`*(a, b: NVectorType): NVectorType =
    if a.length != b.length:
        raise newException(ValueError, "NVector must be of same length for addition")
    result = newNVector(a.length)
    N_VLinearSum_Serial(1.0, a.rawVector[], 1.0, b.rawVector[], result.rawVector[])

proc `+`*(v: NVectorType, c: realtype): NVectorType =
    result = newNVector(v.length)
    N_VAddConst_Serial(v.rawVector[], c, result.rawVector[])

proc `+`*(c: realtype, v: NVectorType): NVectorType =
    result = newNVector(v.length)
    N_VAddConst_Serial(v.rawVector[], c, result.rawVector[])

proc `+=`*(a: var NVectorType, b: NVectorType) =
    if a.length != b.length:
        raise newException(ValueError, "NVector must be of same length for addition")
    N_VLinearSum_Serial(1.0, a.rawVector[], 1.0, b.rawVector[], a.rawVector[])

proc `+=`*(v: var NVectorType, c: realtype) =
    N_VAddConst_Serial(v.rawVector[], c, v.rawVector[])


proc `-`*(a, b: NVectorType): NVectorType =
    if a.length != b.length:
        raise newException(ValueError, "NVector must be of same length for addition")
    result = newNVector(a.length)
    N_VLinearSum_Serial(1.0, a.rawVector[], -1.0, b.rawVector[], result.rawVector[])

proc `-`*(v: NVectorType, c: realtype): NVectorType =
    result = newNVector(v.length)
    N_VAddConst_Serial(v.rawVector[], -c, result.rawVector[])

# TODO
#[
proc `-`*(c: realtype, v: NVectorType): NVectorType =
    result = newNVector(v.length)
    N_VAddConst_Serial(v.rawVector[], c, result.rawVector[])
]#

proc `-=`*(a: var NVectorType, b: NVectorType) =
    if a.length != b.length:
        raise newException(ValueError, "NVector must be of same length for addition")
    N_VLinearSum_Serial(1.0, a.rawVector[], -1.0, b.rawVector[], a.rawVector[])

proc `-=`*(v: var NVectorType, c: realtype) =
    N_VAddConst_Serial(v.rawVector[], -c, v.rawVector[])

proc `*`*(v: NVectorType, c: realtype): NVectorType =
    result = newNVector(v.length)
    N_VScale_Serial(c, v.rawVector[], result.rawVector[])
proc `*`*(c: realtype, v: NVectorType): NVectorType =
    result = newNVector(v.length)
    N_VScale_Serial(c, v.rawVector[], result.rawVector[])


template CVodeProc*(name, body: untyped): untyped {.dirty.} =
    proc `name`(t: realtype, y_raw: N_Vector, ydot_raw: N_Vector, user_data: pointer): cint {.cdecl.} =
        let y = newNVector(y_raw)
        var ydot = newNVector(ydot_raw)
        body
        NV_DATA_S(ydot_raw) = NV_DATA_S(ydot.rawVector[])
    proc `name`(t: realtype, y: NVectorType): NVectorType =
        var ydot = newNVector(y.length)
        body
        result = ydot


echo "\n\n\n\nLet the real testing begin:"
echo "Addition:"
var v1 = newNVector([1.0, 2.0, 3.0])
var v2 = newNVector([4.0, 5.0, 6.0])
var vsum = v1 + v2
N_VPrint_Serial(vsum.rawVector[])
N_VPrint_Serial((1.0 + v1 + 1.0).rawVector[])
echo "Inplace Addition:"
v2 += v1
N_VPrint_Serial(v2.rawVector[])
v2 += 2.2
N_VPrint_Serial(v2.rawVector[])
echo v2 == v1

echo v2
echo v2[2]
echo -1.0 * v2


proc f(t: realtype, y_raw: N_Vector, ydot_raw: N_Vector, user_data: pointer): cint {.cdecl.} =
    let y = newNVector(y_raw)
    var ydot = newNVector(ydot_raw)
    ydot = clone(y)
    NV_DATA_S(ydot_raw) = NV_DATA_S(ydot.rawVector[])
    #NV_Ith_S(ydot_raw, 0) = NV_Ith_S(y_raw, 0)
    #NV_Ith_S(ydot_raw, 1) = NV_Ith_S(y_raw, 1)
    #NV_Ith_S(ydot_raw, 2) = NV_Ith_S(y_raw, 2)

CVodeProc f2:
    #ydot = 1/(t+1) * y
    ydot = y

var cvode_mem: pointer = nil
cvode_mem = CVodeCreate(CV_ADAMS)
var t0 = 0.0
var y0 = newNVector([1.0, 1.0, 1.0])
var A: SUNMatrix = nil#SUNDenseMatrix(3, 3)
let reltol = 1e-8
let abstol = 1e-8
var flag = CVodeInit(cvode_mem, f2, t0, y0.rawVector[])
flag = CVodeSStolerances(cvode_mem, reltol, abstol)
var LS = SUNLinSol_SPGMR(y0.rawVector[], 0, 0)
flag = CVodeSetLinearSolver(cvode_mem, LS, A)
var t: realtype = 0.0
var tout = 1.0
flag = CVode(cvode_mem, tout, y0.rawVector[], addr(t), CV_NORMAL)

import math
var correct = newNVector(@[exp(1.0), exp(1.0), exp(1.0)])
echo "Answer", y0
echo "Error: ", correct - y0
echo "f2: ", f2(10.0, y0)
echo "Seq: ", @y0
CVodeFree(addr(cvode_mem))
flag = SUNLinSolFree(LS)
SUNMatDestroy(A)

# tout should be openArray[realtype]
proc CVodeSolve*(f: proc(t: realtype, y_raw: N_Vector, ydot_raw: N_Vector, user_data: pointer): cint {.cdecl.}, y0: NVectorType, t0, tout: realtype, abstol = 1e-8, reltol = 1e-8): NVectorType =
    var cvode_mem = CVodeCreate(CV_ADAMS)
    var A: SUNMatrix = nil #SUNDenseMatrix(neqs, neqs)
    var flag = CVodeInit(cvode_mem, f, t0, y0.rawVector[])
    flag = CVodeSStolerances(cvode_mem, reltol, abstol)
    var LS = SUNLinSol_SPGMR(y0.rawVector[], 0, 0)
    flag = CVodeSetLinearSolver(cvode_mem, LS, A)
    var t = 0.0
    var y = newNVector(y0.length)
    flag = CVode(cvode_mem, tout, y.rawVector[], addr(t), CV_NORMAL)
    CVodeFree(addr(cvode_mem))
    flag = SUNLinSolFree(LS)
    SUNMATDestroy(A)
    result = y