Closed alkino closed 2 years ago
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:
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
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.
Line 296 is the ond of the file. It does not make much sense to me
do not know how it works in mod2c.
The sparse_thread
callback to int state(so, rhs, _threadargs_)
where #define nrn_state _nrn_state__cdp5
is 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_);
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...
@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.
I turn the second question to @iomaganaris since pramod is on vacation (and will remain on vacation for a while)
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
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:And so, next the NeuronSolver but,
sparse
method is not yet supported with it in nmodl.