Invalid sparse method

[alkino]

In this mod file: https://github.com/BlueBrain/nmodldb/blob/downloader/models/db/modeldb/239421/mod/cdp5.mod

We try to solve with the sparse method with the sympy solver but:

[NMODL] [warning] :: SympySolverVisitor :: solve_lin_system python exception: nonlinear term encountered: -BTC*b1*ca*dsqvol*dt/(diam**2*vrat)

And so, next the NeuronSolver but, sparse method is not yet supported with it in nmodl.

[NMODL] [error] :: NeuronSolveVisitor :: solver method 'sparse' not supported

[cattabiani]

For sure that is not unintentional. I can see in the code:

    /// solve linear system (for "sparse" and "LINEAR")
    void solve_linear_system(const std::vector<std::string>& pre_solve_statements = {});

and

        if (solve_method == codegen::naming::SPARSE_METHOD) {
            solve_linear_system(pre_solve_statements);
        } else if (solve_method == codegen::naming::DERIVIMPLICIT_METHOD) {
            solve_non_linear_system(pre_solve_statements);
        } else {
            logger->error("SympySolverVisitor :: Solve method {} not supported", solve_method);
        }

Nonlinear systems are meant to be solved with the derivimplicit method. There is a check that is actually the case and we throw an error because we selected the sparse method. Now, I do not know how it works in mod2c. I suspect there is a check under-the-hood and it does actually use a derivimplicit method without telling you (or maybe gives just a warning). Otherwise it could just compute some stuff as if it were linear with a 0 approximation.

Quick stuff that you can check:

• do you get the same results with mod2c with derivimplicit?
• and nmodl with derivimplicit?

If they are all the same (numerically equal) there is a good chance that mod2c is correcting the method under-the-hood.

I would like to find more documentation on the sparse method (because it is some time and I forgot what it does exactly). However, I cannot find it. Do you know where I can look? @alkino

[cattabiani]

As by @alkino comment:

when I replace sparse with derivimplicit with mod2c:

Method derivimplicit can't be used with Block state at line 296 in file /gpfs/bbp.cscs.ch/home/cornu/TEST/nmodldb/models/db/modeldb/239421/mod/cdp5.mod

With nmodl the generation is without error.

[cattabiani]

Line 296 is the ond of the file. It does not make much sense to me

[nrnhines]

do not know how it works in mod2c.

The sparse_thread callback to int state(so, rhs, _threadargs_) where #define nrn_state _nrn_state__cdp5is indirect because gpus (at least in the past) could not use pointers to functions.
The key comment is in coreneuron/mechanism/mech/dimplic.cpp

/** When we use solve methods like euler, newton or kinetic schemes,
 *  the state/current updates function need to call solver methods
 *  defined in coreneuron. This is typically done via function pointers.
 *  But for GPU implementation using OpenACC, we can't pass function
 *  pointers.The temporary "workaround" for this was to generate
 *  switch case to select the proper callback function. This is implemented
 *  using python script that look into translated file and generate
 *  _kinderiv.h which has cases for steer functions defined below.
 *  This allows OpenACC to select gpu implementations at compile
 *  time.
 *  \todo: eulerfun/difun are legacy macros and can be replaced with
 *         actual steer function for euler/derivimplicit methods.
 */

The nrn_state function calls

sparse_thread((SparseObj*)_thread[_spth1]._pvoid, 19, _slist1, _dlist1, &t, dt, _kinetic_state_cdp5, _linmat1, _threadargs_);

with #define _linmat1 0 and the corresponding arg in sim/scopmath/sparse_thread.cpp is linflag. The latter is used in a loop to setup the (linearized) matrix equation with spfun(fun, so, so->rhs) solved by matsol(so, _iml) until the solution is close to 0. Here, the fun arg is _kinetic_state_cdp5 which from _kinderiv.h steers the call to the above int state ( function in the translated mod file. _kinderiv.h was generated by coreneuron/coreneuron/kinderiv.py which was called by coreneuron/extra/nrnivmodl_core_makefile.in

#define _NRN_KINETIC_CASES \
  case _kinetic_state_cdp5: state_cdp5(so, rhs, _threadargs_); break; \

The function that does the steering is coreneuron/mechanism/mech/dimplic.cpp:int nrn_kinetic_steer(int fun, SparseObj* so, double* rhs, _threadargsproto_) {
and coreneuron/mechanism/mech/mod2c_core_thread.hpp:#define spfun(arg1, arg2, arg3) nrn_kinetic_steer(arg1, arg2, arg3, _threadargs_);

[cattabiani]

So, if I understand correctly, the sparse method uses a newton iteration linearizing the non-linear problem right? In that case the nonlinear solver (it uses the newton method) should be used in nmodl not the linear one. I still do not remember why we decided to have the sparse solver with the linear system solver...

[nrnhines]

@cattabiani Correct for the first. And yes, if you fallback to the classic mod2c solver then it should be with arg linmod=0.
@pramodk Is it still the case on gpu that one cannot pass function pointers with openacc? Our current work around is quite obscure.

[cattabiani]

I turn the second question to @iomaganaris since pramod is on vacation (and will remain on vacation for a while)

[iomaganaris]

Is it still the case on gpu that one cannot pass function pointers with openacc? Our current work around is quite obscure.

I created a small example with a function pointer being used in an OpenACC region and it seems like this is still the case