diff --git a/scripts/PassiveCable.jl b/scripts/PassiveCable.jl new file mode 100644 index 0000000000000000000000000000000000000000..7c50a6457eb0c28b34c2fb9f407b3c38442ab74d --- /dev/null +++ b/scripts/PassiveCable.jl @@ -0,0 +1,104 @@ +module PassiveCable + +export cable_normalize, cable, rallpack1 + +# Compute solution g(x, t) to +# +# ∂²g/∂x² - g - ∂g/∂t = 0 +# +# on [0, L] × [0,∞), subject to: +# +# g(x, 0) = 0 +# ∂g/∂x (0, t) = 1 +# ∂g/∂x (L, t) = 0 +# +# (This implementation converges slowly for small t) + +function cable_normalized(x, t, L; tol=1e-8) + if t<=0 + return 0.0 + else + ginf = -cosh(L-x)/sinh(L) + sum = exp(-t/L) + + for k = countfrom(1) + a = k*pi/L + e = exp(-t*(1+a^2)) + + sum += 2/L*e*cos(a*x)/(1+a^2) + resid_ub = e/(L*a^3*t) + + if resid_ub<tol + break + end + end + return ginf+sum + end +end + + +# Compute solution f(x, t) to +# +# λ²∂²f/∂x² - f - τ∂f/∂t = 0 +# +# on [0, L] x [0, ∞), subject to: +# +# f(x, 0) = V +# ∂f/∂x (0, t) = I·r +# ∂f/∂x (L, t) = 0 +# +# where: +# +# λ² = 1/(r·g) length constant +# τ = r·c time constant +# +# In the physical model, the parameters correspond to the following: +# +# L: length of cable +# r: linear axial resistivity +# g: linear membrane conductance +# c: linear membrane capacitance. +# V: membrane reversal potential +# I: injected current on the left end (x = 0) of the cable. + +function cable(x, t, L, lambda, tau, r, V, I; tol=1e-8) + scale = I*r*lambda; + if scale == 0 + return V + else + tol_n = abs(tol/scale) + return scale*cable_normalized(x/lambda, t/tau, L/lambda, tol=tol_n) + V + end +end + + +# Rallpack 1 test +# +# One sided cable model with the following parameters: +# +# RA = 1 Ω·m bulk axial resistivity +# RM = 4 Ω·m² areal membrane resistivity +# CM = 0.01 F/m² areal membrane capacitance +# d = 1 µm cable diameter +# EM = -65 mV reversal potential +# I = 0.1 nA injected current +# L = 1 mm cable length +# + +function rallpack1(x, t; tol=1e-8) + RA = 1 + RM = 4 + CM = 1e-2 + d = 1e-6 + EM = -65e-3 + I = 0.1e-9 + L = 1e-3 + + r = 4*RA/(pi*d*d) + lambda = sqrt(d/4 * RM/RA) + tau = CM*RM + + return cable(x, t, L, lambda, tau, r, EM, -I, tol=tol) +end + +end #module diff --git a/scripts/README.md b/scripts/README.md index e645a762eaf88c5f7196a995502589da44c2f89c..230d5be45661f32f8fc7dc213f9b90f71c8585f1 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -78,4 +78,60 @@ of the timeseries data. Use the `-o` or `--output` option to save the plot as an image, instead of displaying it interactively. +profstats +--------- +`profstats` collects the profiling data output from multiple MPI ranks and performs +a simple statistical summary. + +Input files are in the JSON format emitted by the profiling code. + +By default, `profstats` reports the quartiles of the times reported for each +profiling region and subregion. With the `-r` option, the collated raw times +are reported instead. + +Output is in CSV format. + +PassiveCable.jl +--------------- + +Compute analytic solutions to the simple passive cylindrical dendrite cable +model with step current injection at one end from t=0. + +This is used to generate validation data for the first Rallpack test. + +Module exports the following functions: + + * `cable_normalized(x, t, L; tol)` + + Compute potential _V_ at position _x_ in [0, _L_] at time _t_ ≥ 0 according + to the normalized cable equation with unit length constant and time constant. + + Neumann boundary conditions: _V'_(0) = 1; _V'_(L) = 0. + Initial conditions: _V_( _x_, 0) = 0. + + Absolute tolerance `tol` defaults to 1e-8. + + * `cable(x, t, L, lambda, tau, r, V, I; tol)` + + Compute the potential given: + * length constant `lambda` + * time constant `tau`, + * axial linear resistivity `r` + * injected current of `I` at the origin + * reversal potential `V` + + Implied units must be compatible, e.g. SI units. + + Absolute tolerance `tol` defaults to 1e-8. + + * `rallpack1(x, t; tol)` + + Compute the value of the potential in the Rallpack 1 test model at position + `x` and time `t`. + + Parameters for the underlying cable equation calculation are taken from the + Rallpack model description in SI units; as the cable length is 1 mm in this + model, `x` can take values in [0, 0.001]. + + Absolute tolerance `tol` defaults to 1e-8. diff --git a/scripts/profstats b/scripts/profstats index 33f37356836104e57758f0407e1da5fd0625a111..88f68c72e6253fa4c0d240d051d7ce1be1960604 100755 --- a/scripts/profstats +++ b/scripts/profstats @@ -1,4 +1,4 @@ -#!env python +#!/usr/bin/env python2 #coding: utf-8 import json diff --git a/scripts/tsplot b/scripts/tsplot index 9ad6e2a5c714a5f634e37c31c88411f63923f130..8e89cbda96451680bfc9de8087ba558c0282defc 100755 --- a/scripts/tsplot +++ b/scripts/tsplot @@ -1,4 +1,4 @@ -#!env python +#!/usr/bin/env python2 #coding: utf-8 import argparse