Diffusion: tests, fixed coefficients in the solver, & tutorial #2226
Conversation
…pointing to the fact that the diffusion across segments of different radius is still erroneous; see arbor-sim#2145)
…ternal-catalogues
…pointing to the fact that the diffusion across segments of different radius is still erroneous; see arbor-sim#2145)
…de readability improved
…ent_radii' again (for as long as the related fix is not merged)
|
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. |
arbor/fvm_layout.cpp
Outdated
| @@ -436,14 +429,15 @@ fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global | |||
| cv_length += embedding.integrate_length(cable); | |||
| } | |||
|
|
|||
| if (cv_length>0) { | |||
There was a problem hiding this comment.
This should not be 0 in the first place... Would be good to ensure this upstream and otherwise throw a warning.
There was a problem hiding this comment.
This can happen at certain points and would lead to div by zero.
There was a problem hiding this comment.
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).
| ====================================== | ||
|
|
||
| Originally we intended for the axial diffusion process to be fully integrated | ||
| into the cable model. However, for various technical reasons, this has proven |
There was a problem hiding this comment.
| 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 |
| is subject to diffusion as opposed to the internal concentration. | ||
|
|
||
| This tutorial will show you how to build a simulation that fully integrates | ||
| diffusive ion dynamics with a cable cell neuron. While we focus on a single |
There was a problem hiding this comment.
| 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 |
|
|
||
| This tutorial will show you how to build a simulation that fully integrates | ||
| diffusive ion dynamics with a cable cell neuron. While we focus on a single | ||
| neuron, this should simply work as is in a network of cells. |
There was a problem hiding this comment.
| neuron, this should simply work as is in a network of cells. | |
| neuron, this will naturally work as well in a network of cells. |
jlubo
changed the title
|
Regarding GPU, I did run it on one of our machines and the tests seemed to pass. However, As for the fix or warning: I'd rather have a temporary and slow fix via atomics. I'll think about it. |
|
That said, I wanted to take a long, hard look at reduce-by-key for a while anyhow. |
Review suggestions Co-authored-by: Jannik Luboeinski <[email protected]>
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:$\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:
$\sim \ln\frac{X_i}{X_o}$
We construct a coupling term from the cable equation to the ionic diffusion via
iX = g(U-eX)this is the conductance model and should be found in an NMODL file and
eXwhich is the Nernst equation found, again, in a builtin NMODL file. But, for technical reasons we had to invent a special diffusive concentration
Xdwhich is not identical toXi. For proper behaviour, we should have usedXdinstead ofXiin 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