Diffusion: tests, fixed coefficients in the solver, & tutorial

[thorstenhater]

By careful re-analysis of the diffusion solver we found that the coeffiecients of the LHS need an additional factor of the CV Volume.

Scopes be creeping. During the analysis it became apparent that the during coupling to the cable equation proper raises at lot of potential problems.

Conversion of ionic current densities $i_X$ to mass transfer is $\dot c_X = \frac{i_XA}{q F V}$ with volume, area, Faraday's constant, and charge. Now, people are interested in using neutral species $n$, ie $q=0$ for which we also should expect $i_n=0$. Yet we need to special case on this (and assert zero current) to avoid ill-defined division. That's the smaller problem.

The more worrying issue is this:
We construct a coupling term from the cable equation to the ionic diffusion via $\dot c_X = \frac{i_XA}{q F V}$. Yet, there's no backreaction to the cable equation, unless the user explicitly constructs it. How? Well, the way ions should feed back to the voltage is this:
iX = g(U-eX)
this is the conductance model and should be found in an NMODL file and
eX $\sim \ln\frac{X_i}{X_o}$
which is the Nernst equation found, again, in a builtin NMODL file. But, for technical reasons we had to invent a special diffusive concentration Xd which is not identical to Xi. For proper behaviour, we should have used Xd instead of Xi in the Nernst term above. So, a non-standard Nernst module needs to be used.

This brings me to the great change in this context: remove the offending term and carefully document how to retrofit it in NMODL and add the proper Nernst model, too.

fixes #2145
requires #2209

[jlubo]

The functionality seems fine, at least on CPU single-core systems! I've checked this by running several tests. First of all, and most importantly, all tests from #2209 pass, including the previously commented ones (I guess @thorstenhater you checked this as well). Also the tests of my single- and multi-compartment network simulations (which are for a more particular model, but relatively extensive) all pass. The code too looks fine to me, too, as far as I can judge. Regarding GPU systems, I guess #2244 will be needed as well!? What else will be needed here?

It's also very nice to have the tutorial! However, I would suggest to split this into a documentation and a tutorial page. I have a draft for a documentation page on the diffusion mechanism that I started to write a while ago - maybe we could create another PR for that and put parts of the tutorial here into that one.

[jlubo]

This should not be 0 in the first place... Would be good to ensure this upstream and otherwise throw a warning.

[thorstenhater]

This can happen at certain points and would lead to div by zero.

[jlubo]

I just tested the latest version of this with the tests from python/test/unit/test_diffusion.py, and on CPU, the tests still pass.

However, as #2244 is still pending, I guess this is not supposed to run on GPU yet, right? I find the tests to fail on GPU with the following output:

======================================================================
FAIL: test_diffusion_different_length (test.unit.test_diffusion.TestDiffusion)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/jannik/arbor_repo/python/test/fixtures.py", line 40, in inject_fixture
    return func(*bound.args, **bound.kwargs)
  File "/home/jannik/arbor_repo/python/test/unit/test_diffusion.py", line 440, in test_diffusion_different_length
    self.simulate_and_test_diffusion(
  File "/home/jannik/arbor_repo/python/test/unit/test_diffusion.py", line 393, in simulate_and_test_diffusion
    self.assertAlmostEqual(
AssertionError: 0.0013333333333262613 != 1.5915494309189535 within 0.015915494309189534 delta (1.5902160975856272 difference)

======================================================================
FAIL: test_diffusion_different_radii (test.unit.test_diffusion.TestDiffusion)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/jannik/arbor_repo/python/test/fixtures.py", line 40, in inject_fixture
    return func(*bound.args, **bound.kwargs)
  File "/home/jannik/arbor_repo/python/test/unit/test_diffusion.py", line 455, in test_diffusion_different_radii
    self.simulate_and_test_diffusion(
  File "/home/jannik/arbor_repo/python/test/unit/test_diffusion.py", line 393, in simulate_and_test_diffusion
    self.assertAlmostEqual(
AssertionError: 0.0012307692307876126 != 0.48970751720583183 within 0.004897075172058319 delta (0.4884767479750442 difference)

======================================================================
FAIL: test_diffusion_equal_radii (test.unit.test_diffusion.TestDiffusion)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/jannik/arbor_repo/python/test/fixtures.py", line 40, in inject_fixture
    return func(*bound.args, **bound.kwargs)
  File "/home/jannik/arbor_repo/python/test/unit/test_diffusion.py", line 425, in test_diffusion_equal_radii
    self.simulate_and_test_diffusion(
  File "/home/jannik/arbor_repo/python/test/unit/test_diffusion.py", line 393, in simulate_and_test_diffusion
    self.assertAlmostEqual(
AssertionError: 0.0013333333333262613 != 1.5915494309189535 within 0.015915494309189534 delta (1.5902160975856272 difference)

----------------------------------------------------------------------
Ran 3 tests in 0.695s

FAILED (failures=3)

So if we want to merge this now for CPU only, I'd suggest to throw a warning when someone tries to run this on GPU.

Apart from that, considering the code and the tutorial, I don't see any obvious issues (just suggested some small adjustments).

[jlubo]

Suggested change

	into the cable model. However, for various technical reasons, this has proven
	into the cable model (linked to the internal concentration, see :ref:`here <cablecell-decoration>`). However, for various technical reasons, this has proven

[jlubo]

Suggested change

	diffusive ion dynamics with a cable cell neuron. While we focus on a single
	diffusive particle dynamics with a cable cell neuron. While we focus on a single

[jlubo]

Suggested change

	neuron, this should simply work as is in a network of cells.
	neuron, this will naturally work as well in a network of cells.

[thorstenhater]

Regarding GPU, I did run it on one of our machines and the tests seemed to pass. However,
I cannot claim surprise, as stochastic failures are precisely what is expected.

As for the fix or warning: I'd rather have a temporary and slow fix via atomics. I'll think about it.
Let's assume for now that we either get a simple fix or the issue will stand as 'Known Issue'
until the next release.

[thorstenhater]

That said, I wanted to take a long, hard look at reduce-by-key for a while anyhow.

[jlubo]

As a note, we should include fixes for the documentation to state the correct units for ion/particle concentrations and diffusivity constants (see issue #2084).