Setting and Confirming NEURON's Solver

[ajpc6d]

I'm trying to compare the results of different solvers offered in NEURON. My strategy goes something like this:

Code: Select all

moi_switch = 1
match moi_switch:
    case 0:
        h.CVode().use_daspk()
        self.cell.integration_method = 'SUNDIALS - DASPK'
    case 1:
        h.secondorder = 2
        self.cell.integration_method = 'Crank-Nicholson'
    case 2:
        h.CVode().use_mxb(0)
        self.cell.integration_method = 'SUNDIALS - MXB (structured)'
    case 3:
        h.CVode().use_mxb(1)
        self.cell.integration_method = 'SUNDIALS - MXB (sparse)'
    case 4:
        h.CVode().stiff(2)
        self.cell.integration_method = 'SUNDIALS - backwards Euler'
    case 5:
        h.CVode().stiff(0)
        self.cell.integration_method = 'Adams-Bashforth'
    case _:
    	self.cell.integration_method = 'Backwards Euler'
        pass

But I don't know how I would confirm that what I'm setting is the exclusive solver at work. Is there a query I can send to NEURON to confirm this? I think setting h.secondorder = 2 is correct, because the results are noticeably different. All the others give identical output, as far as I can tell - but is this due to convergence, or because my code isn't doing what I think it's doing?

[ajpc6d]

Update: I ran a percent error calculation comparing the results of all the solver settings in the post above. For all conditions except Crank-Nicholson, the error between any two solver settings is 0%. I interpret this to mean I have failed to actually change the solver, since I would expect at least negligible differences in results between different methods.

Update 2: I also ran several comparison tests between the backward Euler method and the Crank-Nicholson method. Using constant time steps down to 10 ns duration, the two methods never converged toward a common answer. To my understanding, both become more stable with decreasing step size. So I'm at a loss for what could be happening here?

[ajpc6d]

Update 3: Based on documentation here, I think (part of) the correct command order is:

Code: Select all

h.CVode().active(True)
h.CVode().use_daspk(True) # or whatever solver you want..

This appears to produce (negligible) differences in simulation output with constant time step methods, suggesting the solver has now been changed.
However, the output remains exactly identical for variable time step methods, suggesting... je ne sais quoi. It seems unlikely that time step changes issued through Cvode would be incompatible with solvers set through Cvode.

[ajpc6d]

Update 4:
I believe I've solved this. For anyone else facing the same problem, I recommend the following procedure:

Code: Select all

from neuron import h
cvode = h.CVode()

# let's say we want to activate the DASPK solver
cvode.active(True)
cvode.use_daspk(True)
cvode.use_mxb(1) # I think this is optional, in that NEURON [i]should[/i] adjust the setting automatically
covde.stiff(0) # activating Adams-Bashforth integration -- also optional, to be determined by your integration needs

# This is what I was looking for originally: some way to confirm the solver has actually been tuned to the desired settings
# this function prints out a numerical code corresponding to the solver settings. See [url=https://nrn.readthedocs.io/en/8.2.4/python/simctrl/cvode.html#CVode.current_method]this page[/url] for the key
cvode.current_method()

# do stuff in the code
# run a simulation in NEURON

cvode.statistics() # a nifty feature that prints out a slew of information about what CVode just did under the hood

# now, if we want to run a second simulation in the same code, everything must be deactivated manually
cvode.active(False)
cvode.stiff(2) # back to the default value
cvode.use_mxb(0) # back to the default value
cvode.use_daspk(False)
cvode.re_init()

# now we're ready to choose another solver option
h.secondorder = 2

If any of this is found to be incorrect, please do let me know!