I have been trying to get a "minimal" working example together for an implementation of gap junctions over MPI before moving on to a more complex system. Well actually, I got the more complex system working using MPI-aware syntax (i.e. pc.source_var, pc.target_var) on a single process, but when I increase the number of processes it throws a segfault, so I have tried to get a simpler model working to help understand where the problem is. However, I am still running into trouble and I am not even sure whether it is the same trouble I was having before:)
Here is the MWE code I have written (NB: at this stage it is a "rectified gap junction" where the current only flows one way for testing purposes, whereas usually you will want to have a bidirectional gap junction. In that case you will need another gap junction connection going from the "post-synaptic" cell to the "pre-synaptic" cell):
Code: Select all
"""
A minimum working example of a NEURON gap junction over MPI
Author: Tom Close
Date: 8/1/2013
Email: tclose@oist.jp
"""
import os
import argparse
import numpy as np
# This is a hack I use on our cluster, to get MPI initialised=True. There is probably something
# wrong with our setup but I can't be bothered trying to work out what it is at this point. All
# suggestions welcome :)
try:
from mpi4py import MPI #@UnresolvedImport @UnusedImport
except:
print "mpi4py was not found, MPI will remain disabled if MPI initialized==False on startup"
from neuron import h, load_mechanisms
# Not sure this is necessary, or whether I can just use h.finitialize instead of h.stdinit
h.load_file('stdrun.hoc')
# The GID used for the gap junction connection. NB: this number is completely independent from the
# GID's used for NEURON sections.
GID_FOR_VAR = 0
# Arguments to the script
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--plot', action='store_true', help="Plot the data instead of saving it")
parser.add_argument('--output_dir', type=str, default=os.getcwd(),
help="The directory to save the output files into")
parser.add_argument('--gap_mechanism_dir', type=str, default=os.getcwd(),
help="The directory to load the gap mechanism from")
args = parser.parse_args()
# Load gap mechanism from another directory if required
if args.gap_mechanism_dir is not os.getcwd():
load_mechanisms(args.gap_mechanism_dir)
# Get the parallel context and related parameters
pc = h.ParallelContext()
num_processes = int(pc.nhost())
mpi_rank = int(pc.id())
print "On process {} of {}".format(mpi_rank+1, num_processes)
print "Creating test network..."
# The pre-synaptic cell is created on the first node and the post-synaptic cell on the last node
# (NB: which will obviously be the same if there is only one node)
if mpi_rank == 0:
print "Creating pre-synaptic cell on process {}".format(mpi_rank)
# Create the pre-synaptic cell
pre_cell = h.Section()
pre_cell.insert('pas')
# Connect the voltage of the pre-synaptic cell to the gap junction on the post-synaptic cell
pc.source_var(pre_cell(0.5)._ref_v, GID_FOR_VAR)
# Stimulate the first cell to make it obvious whether gap junction is working
stim = h.IClamp(pre_cell(0.5))
stim.delay = 50
stim.amp = 10
stim.dur = 100
# Record Voltage of pre-synaptic cell
pre_v = h.Vector()
pre_v.record(pre_cell(0.5)._ref_v)
if mpi_rank == (num_processes - 1):
print "Creating post-synaptic cell on process {}".format(mpi_rank)
# Create the post-synaptic cell
post_cell = h.Section()
post_cell.insert('pas')
# Insert gap junction
gap_junction = h.gap(0.5, sec=post_cell)
gap_junction.g = 1.0
# Connect gap junction to pre-synaptic cell
pc.target_var(gap_junction._ref_vgap, GID_FOR_VAR)
# Record Voltage of post-synaptic cell
post_v = h.Vector()
post_v.record(post_cell(0.5)._ref_v)
# Finalise construction of parallel context
pc.setup_transfer()
# Record time
rec_t = h.Vector()
rec_t.record(h._ref_t)
print "Finished network construction on process {}".format(mpi_rank)
# Run simulation
print "Setting maxstep on process {}".format(mpi_rank)
pc.set_maxstep(10)
print "Finitialise on process {}".format(mpi_rank)
#h.finitialize(-60)
h.stdinit()
print "Solving on process {}".format(mpi_rank)
pc.psolve(100)
print "Running worker on process {}".format(mpi_rank)
pc.runworker()
print "Completing parallel context on process {}".format(mpi_rank)
pc.done()
print "Finished run on process {}".format(mpi_rank)
# Convert recorded data into Numpy arrays
t_array = np.array(rec_t)
if mpi_rank == 0:
pre_v_array = np.array(pre_v)
if mpi_rank == (num_processes - 1):
post_v_array = np.array(post_v)
# Either plot the recorded values
if args.plot and num_processes == 1:
print "Plotting..."
import matplotlib.pyplot as plt
if mpi_rank == 0:
pre_fig = plt.figure()
plt.plot(t_array, pre_v_array)
plt.title("Pre-synaptic cell voltage")
plt.xlabel("Time (ms)")
plt.ylabel("Voltage (mV)")
if mpi_rank == (num_processes - 1):
pre_fig = plt.figure()
plt.plot(t_array, post_v_array)
plt.title("Post-synaptic cell voltage")
plt.xlabel("Time (ms)")
plt.ylabel("Voltage (mV)")
plt.show()
else:
# Save data
print "Saving data..."
if mpi_rank == 0:
np.savetxt(os.path.join(args.output_dir, "pre_v.dat"),
np.transpose(np.vstack((t_array, pre_v_array))))
if mpi_rank == (num_processes - 1):
np.savetxt(os.path.join(args.output_dir, "post_v.dat"),
np.transpose(np.vstack((t_array, post_v_array))))
print "Done."
Code: Select all
: gap.mod
: This is a conductance based gap junction model rather
: than resistance because Traub occasionally likes to
: set g=0 which of course is infinite resistance.
NEURON {
SUFFIX gap
POINT_PROCESS gap
RANGE g, i
POINTER vgap
ELECTRODE_CURRENT i
}
PARAMETER { g = 1e-10 (1/megohm) }
ASSIGNED {
v (millivolt)
vgap (millivolt)
i (nanoamp)
}
BREAKPOINT { i = (vgap - v)*g }
Code: Select all
MPI_Initialized==true, enabling MPI functionality.
numprocs=1
NEURON -- Release 7.2 (562:42a47463b504) 2011-12-21
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2008
See http://www.neuron.yale.edu/credits.html
--------------------------------------------------------------------------
An MPI process has executed an operation involving a call to the
"fork()" system call to create a child process. Open MPI is currently
operating in a condition that could result in memory corruption or
other system errors; your MPI job may hang, crash, or produce silent
data corruption. The use of fork() (or system() or other calls that
create child processes) is strongly discouraged.
The process that invoked fork was:
Local host: tombo20501 (PID 23294)
MPI_COMM_WORLD rank: 0
If you are *absolutely sure* that your application will successfully
and correctly survive a call to fork(), you may disable this warning
by setting the mpi_warn_on_fork MCA parameter to 0.
--------------------------------------------------------------------------
Additional mechanisms from files
adexp.mod alphaisyn.mod alphasyn.mod expisyn.mod gap.mod gsfa_grr.mod hh_traub.mod netstim2.mod new_gap.mod refrac.mod reset.mod stdwa_guetig.mod stdwa_softlimits.mod stdwa_songabbott.mod stdwa_symm.mod tmgsyn.mod tmisyn.mod vecstim.mod
On process 1 of 1
Creating test network...
Creating pre cell on process 0
Creating post cell on process 0
Finished network construction on process 0
Setting maxstep on process 0
Finitialise on process 0
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 23294 on node tombo20501 exited on signal 11 (Segmentation fault).
--------------------------------------------------------------------------
I get another error, which is probably unrelated, when I run the script with -np 2,
Code: Select all
MPI_Initialized==true, enabling MPI functionality.
NEURON -- Release 7.2 (562:42a47463b504) 2011-12-21
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2008
See http://www.neuron.yale.edu/credits.html
MPI_Initialized==true, enabling MPI functionality.
numprocs=2
--------------------------------------------------------------------------
An MPI process has executed an operation involving a call to the
"fork()" system call to create a child process. Open MPI is currently
operating in a condition that could result in memory corruption or
other system errors; your MPI job may hang, crash, or produce silent
data corruption. The use of fork() (or system() or other calls that
create child processes) is strongly discouraged.
The process that invoked fork was:
Local host: tombo20502 (PID 8277)
MPI_COMM_WORLD rank: 0
If you are *absolutely sure* that your application will successfully
and correctly survive a call to fork(), you may disable this warning
by setting the mpi_warn_on_fork MCA parameter to 0.
--------------------------------------------------------------------------
Additional mechanisms from files
adexp.mod alphaisyn.mod alphasyn.mod expisyn.mod gap.mod gsfa_grr.mod hh_traub.mod netstim2.mod new_gap.mod refrac.mod reset.mod stdwa_guetig.mod stdwa_softlimits.mod stdwa_songabbott.mod stdwa_symm.mod tmgsyn.mod tmisyn.mod vecstim.mod
On process 2 of 2
Creating test network...
Creating post cell on process 1
On process 1 of 2
Creating test network...
Creating pre cell on process 0
Finished network construction on process 1
Setting maxstep on process 1
Finished network construction on process 0
Setting maxstep on process 0
Finitialise on process 0
Finitialise on process 1
python: netpar.cpp:469: void nrn_spike_exchange_init(): Assertion `usable_mindelay_ >= dt && (usable_mindelay_ * dt1_) < 255' failed.
python: netpar.cpp:469: void nrn_spike_exchange_init(): Assertion `usable_mindelay_ >= dt && (usable_mindelay_ * dt1_) < 255' failed.
--------------------------------------------------------------------------
mpirun noticed that process rank 1 with PID 8278 on node tombo20502 exited on signal 11 (Segmentation fault).
--------------------------------------------------------------------------
Thanks,
Tom