{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Parallel simulations\n", "\n", "
\n", "**Topics**: Parallel solvers, Mesh partitioning.\n", "
\n", "\n", "Parallel simulations are run with the `'TetOpSplit'` solver, which supports parallel stochastic spatial reaction-difusion-EField simulations. Note however that this solver does not support dynamic load balancing.\n", "\n", "For more details about the accuracy and performace of the parallel TetOpSplit solver, please check the following papers:\n", "\n", "Hepburn, I., Chen, W., and De Schutter, E. (2016). **Accurate reaction-diffusion operator splitting on tetrahedral meshes for parallel stochastic molecular simulations**. J. Chem. Phys. 145, 054118–22. doi:10.1063/1.4960034.\n", "\n", "Chen W and De Schutter E (2017). **Parallel STEPS: Large Scale Stochastic Spatial Reaction-Diffusion Simulation with High Performance Computers**. Front. Neuroinform. 11:13. doi: 10.3389/fninf.2017.00013.\n", "\n", "## Converting serial simulations to parallel\n", "\n", "We will focus here on models involving tetrahedral meshes. So far, in these cases, we created the simulation object with the `'Tetexact'` serial solver. In order to parallelize these simulations, we need to decide on a way to partition the tetrahedral mesh so that each part is run in a separate process. The first thing that we need to do is thus to declare this partition.\n", "\n", "### Partitioning the mesh\n", "\n", "#### Linear mesh partition\n", "\n", "Mesh partitions are usually created by using one of the available dedicated functions. The simplest one consists in partitioning the mesh along a 3D grid:\n", "\n", "```python\n", "import steps.interface\n", "\n", "from steps.geom import *\n", "\n", "# ...\n", "\n", "part = LinearMeshPartition(mesh, xbins, ybins, zbins)\n", "```\n", "\n", "In this example, we removed the mesh declaration code in order to focus on the partitioning. The `LinearMeshPartition` function takes the mesh as first parameter and then 3 integers that represent the number of bins along the x, y, and z axis respectively. Assuming the mesh is much longer on the x-axis than on the two others, we could parition the mesh in 10 bins along the x-axis only with `LinearMeshPartition(mesh, 10, 1, 1)`.\n", "\n", "This partition function first paritions tetrahedrons according to their position. If patches are present in the model, their triangles will automatically be partitioned according to the tetrahedron partitions.\n", "\n", "#### Complex mesh partition\n", "\n", "The grid-based partitioning approach may not be suitable for complex geometries, such as a dendritic tree, in this case, we suggest using third party partitioning tools such as [Metis](http://glaros.dtc.umn.edu/gkhome/metis/metis/overview).\n", "To partition a STEPS TetMesh using Metis, we first need to convert the tetrahedron connectivity information in the mesh to Metis format:\n", "\n", "```python\n", "mesh.ConvertToMetis('/path/to/file/mymesh.metis')\n", "```\n", "\n", "We then call Metis from a bash terminal:\n", "\n", "```bash\n", "mpmetis -ncommon=3 -minconn -niter=1000 mymesh.metis 10\n", "```\n", "\n", "Metis will partition the mesh into 10 segments and store the partition information to \"mymesh.metis.epart.10\". Details of Metis parameters can be found in the [Metis user manual](http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/manual.pdf). After the partitioning, you can read the partition information from the file with:\n", "\n", "```python\n", "part = MetisPartition('/path/to/file/mymesh.metis.epart.10')\n", "```\n", "\n", "#### Manual mesh partition\n", "\n", "The `LinearMeshPartition` and `MetisPartition` functions both return a `MeshPartition` object. It is possible to directly construct a `MeshPartition` by explicitely specifying which tetrahedrons and triangles should be attributed to which partition:\n", "\n", "```python\n", "\n", "def getTetPartition(tet):\n", " # ...\n", " \n", "def getTriPartition(tri):\n", " # ...\n", "\n", "tet_hosts = [getTetPartition(tet) for tet in mesh.tets]\n", "\n", "tri_hosts = {}\n", "for patch in mesh.ALL(Patch):\n", " for tri in patch.tris:\n", " tri_hosts[tri.idx] = getTriPartition(tri)\n", "\n", "part = MeshPartition(mesh, tet_hosts, tri_hosts)\n", "```\n", "\n", "In this example, we assume that we have two functions : `getTetPartition` that takes a tetrahedron as a parameter and returns an integer partition index ; and `getTriPartition` that takes a triangle as a parameter and returns an integer partition index. We then create the `tet_hosts` list that will associate each tetrahedron to an integer partition index. All tetrahedrons need to be partitioned so the length of the list is equal to the number of tetrahedrons in the mesh. Triangles only need to be partitioned if they are part of a `Patch`. We thus loop over all `Patch` objects in the mesh and build the `tri_hosts` dictionary that associates a triangle index to an integer partition index.\n", "\n", "Finally, we create the partition object by simply passing the mesh along with `tet_hosts` and `tri_hosts` to the `MeshPartition` constructor. More details are available in the [documentation](API_geom.rst#steps.API_2.geom.MeshPartition).\n", "\n", "#### Visualizing partitions\n", "\n", "If PyQtGraph and PyOpenGL are installed, partitions can be visualized using the visual module (see [corresponding guide](STEPS_Tutorial_Visual.ipynb)):\n", "\n", "```python\n", "from steps.visual import *\n", "\n", "# ...\n", "\n", "sc = SimControl()\n", "with sc:\n", " PartitionDisplay(part)\n", "sc.run()\n", "```\n", "\n", "This will display the tetrahedron partitions, `PartitionDisplay(part, 'tri')` will instead display the triangles partitions. The image below is a comparison between a `LinearMeshPartition` (left) and a `MetisPartition` (right), each color represents a partition to be assigned to an MPI process.\n", "\n", "\n", "\n", "### Using parallel solvers\n", "\n", "After having partitioned the mesh, we only need to change the solver from `'Tetexact'` to `'TetOpSplit'`:\n", "\n", "```python\n", "# ...\n", "sim = Simulation('TetOpSplit', model, mesh, rng, True, part)\n", "# ...\n", "```\n", "\n", "In this example, we took the `Simulation` construction from the previous chapter, changed `'Tetexact'` to `'TetOpSplit'` and added the mesh partition `part` as 6th parameter. These changes alone (create a mesh partition, use the `'TetOpSplit'` solver) are enough to convert a script from serial to parallel.\n", "\n", "If, during the execution, we need to manually save data or print messages to screen, this should only be done in one of the MPI processes. We can check in which MPI process we are by using the `MPI` class fomr the `sim` module:\n", "\n", "```python\n", "from steps.sim import *\n", "\n", "#...\n", "\n", "if MPI.rank == 0:\n", " print(...)\n", "```\n", "`MPI.rank` returns the rank of the current MPI process. We can get the total number of processes with `MPI.nhosts`.\n", "\n", "#### Parallel EField simulations\n", "\n", "The `'TetOpSplit'` solver can compute membrane potentials in different ways. As we saw in the previous chapter, membrane potential computation can be turned on by passing `True` as 5th parameter of the `Simulation` constructor. It is however also, when using `'TetOpSplit'`, to pass any of the following values:\n", "\n", "- `MPI.EF_NONE` (equivalent to `False`): No EField solver is needed;\n", "- `MPI.EF_DEFAULT` (equvalent to `True`): Run serial EField simulation (Tetexact version) on process 0;\n", "- `MPI.EF_DV_BDSYS`: Use parallel SuperLU EField solver;\n", "- `MPI.EF_DV_PETSC`: Use parallel PETSc EField solver.\n", "\n", "For small scale EFIeld simulation we recommand the `MPI.EF_DV_BDSYS` parameter, for large scale EField simulation we recommand the `MPI.EF_DV_PETSC` parameter.\n", "\n", "### Running parallel simulations\n", "\n", "You can run a parallel `'TetOpSplit'` simulation by calling 'mpirun' or similar executable provided by your MPI distribution.\n", "\n", "```bash\n", "mpirun -n NPROCS python3 STEPS_Tutorial_MPI.py\n", "```\n", "\n", "`NPROCS` is the number of MPI processes that should be used.\n", "\n", "### Example\n", "\n", "We now turn to converting the previous chapter simulation from serial to parallel. We use the `'TetOpSplit'` solver with the SuperLU Efield solver:\n", "\n", "The corresponding python scripts: [STEPS_Tutorial_MPI.py](https://github.com/CNS-OIST/STEPS_Example/tree/master/user_manual/source/API_2/scripts/STEPS_Tutorial_MPI.py) and [STEPS_Tutorial_MPI_plot.py](https://github.com/CNS-OIST/STEPS_Example/tree/master/user_manual/source/API_2/scripts/STEPS_Tutorial_MPI_plot.py)\n", "\n", "```python\n", "# Example: Parallel simulations\n", "# http://steps.sourceforge.net/manual/API_2/STEPS_Tutorial_MPI.html\n", "\n", "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #\n", "\n", "import steps.interface\n", "\n", "from steps.model import *\n", "from steps.geom import *\n", "from steps.sim import *\n", "from steps.saving import *\n", "from steps.rng import *\n", "\n", "import numpy as np\n", "import math\n", "\n", "# Potassium conductance = 0.036 S/cm2\n", "\n", "# Potassium single-channel conductance\n", "K_G = 20.0e-12 # Siemens\n", "\n", "# Potassium channel density\n", "K_ro = 18.0e12 # per square meter\n", "\n", "# Potassium reversal potential\n", "K_rev = -77e-3 # volts\n", "\n", "# Sodium conductance = 0.120 S/cm2\n", "\n", "# Sodium single-channel conductance\n", "Na_G = 20.0e-12 # Siemens\n", "\n", "# Sodium channel density\n", "Na_ro = 60.0e12 # per square meter\n", "\n", "# Sodium reversal potential\n", "Na_rev = 50e-3 # volts\n", "\n", "# Leak single-channel conductance\n", "L_G = 0.3e-12 # Siemens\n", "\n", "# Leak density\n", "L_ro = 10.0e12 # per square meter\n", "\n", "# Leak reveral potential\n", "leak_rev = -54.4e-3 # volts\n", "\n", "# A table of potassium channel population factors:\n", "# n0, n1, n2, n3, n4\n", "K_facs = [ 0.21768, 0.40513, 0.28093, 0.08647, 0.00979 ]\n", "\n", "# A table of sodium channel population factors\n", "# m0h0, m1h0, m2h0, m3h0, m0h1, m1h1, m2h1, m3h1:\n", "Na_facs = [[0.34412, 0.05733, 0.00327, 6.0e-05],\n", " [0.50558, 0.08504, 0.00449, 0.00010]]\n", "\n", "# Temperature for gating kinetics\n", "celsius = 20.0\n", "\n", "# Current injection\n", "Iclamp = 50.0e-12 # amps\n", "\n", "# Voltage range for gating kinetics in Volts\n", "Vrange = [-100.0e-3, 50e-3, 1e-4]\n", "\n", "def HHRateFunction(A, B, C, D, F, H, V):\n", " num = A + B * V\n", " denom = C + H * math.exp((V + D) / F)\n", " if num == denom == 0:\n", " return F * B / (H * math.exp((V + D) / F))\n", " else:\n", " return num / denom\n", "\n", "# The simulation dt\n", "DT_sim = 1.0e-4 # seconds\n", "\n", "# The time until which the simulation should be run\n", "ENDT = 4.0e-3\n", "\n", "#########################\n", "# Model setup\n", "#########################\n", "\n", "model = Model()\n", "\n", "r = ReactionManager()\n", "\n", "with model:\n", " ssys = SurfaceSystem.Create()\n", "\n", " # Potassium channel\n", " Ko, Kc = SubUnitState.Create()\n", " KSU = SubUnit.Create([Ko, Kc])\n", " VGKC = Channel.Create([KSU]*4)\n", "\n", " # Sodium channel\n", " Na_mo, Na_mc, Na_hi, Na_ha = SubUnitState.Create()\n", " NamSU, NahSU = SubUnit.Create(\n", " [Na_mo, Na_mc],\n", " [Na_hi, Na_ha]\n", " )\n", " VGNaC = Channel.Create([NamSU, NamSU, NamSU, NahSU])\n", "\n", " # Leak channel\n", " lsus = SubUnitState.Create()\n", " Leak = Channel.Create([lsus])\n", "\n", " thi = math.pow(3.0, ((celsius-6.3)/10.0))\n", "\n", " _a_n = VDepRate(\n", " lambda V: thi * 1e3 * HHRateFunction(-0.55, -0.01, -1, 55, -10, 1, V*1e3),\n", " vrange=Vrange\n", " )\n", " _b_n = VDepRate(\n", " lambda V: thi * 1e3 * HHRateFunction(1, 0, 0, 65, 80, 8, V*1e3),\n", " vrange=Vrange\n", " )\n", "\n", " _a_m = VDepRate(\n", " lambda V: thi * 1e3 * HHRateFunction(-4, -0.1, -1, 40, -10, 1, V*1e3),\n", " vrange=Vrange\n", " )\n", " _b_m = VDepRate(\n", " lambda V: thi * 1e3 * HHRateFunction(1, 0, 0, 65, 18, 0.25, V*1e3),\n", " vrange=Vrange\n", " )\n", "\n", " _a_h = VDepRate(\n", " lambda V: thi * 1e3 * HHRateFunction(1, 0, 0, 65, 20, 1 / 0.07, V*1e3),\n", " vrange=Vrange\n", " )\n", " _b_h = VDepRate(\n", " lambda V: thi * 1e3 * HHRateFunction(1, 0, 1, 35, -10, 1, V*1e3),\n", " vrange=Vrange\n", " )\n", "\n", " with ssys:\n", " with VGKC[...]:\n", " Kc.s Ko.s\n", " r[1].K = _a_n, _b_n\n", "\n", " with VGNaC[...]:\n", " Na_hi.s Na_ha.s\n", " r[1].K = _a_h, _b_h\n", "\n", " Na_mc.s Na_mo.s\n", " r[1].K = _a_m, _b_m\n", " \n", " VGKC_I = OhmicCurr.Create(VGKC[Ko, Ko, Ko, Ko], K_G, K_rev)\n", " VGNaC_I = OhmicCurr.Create(VGNaC[Na_mo, Na_mo, Na_mo, Na_ha], Na_G, Na_rev)\n", " Leak_I = OhmicCurr.Create(Leak[lsus], L_G, leak_rev)\n", "\n", "#########################\n", "# Geom setup\n", "#########################\n", "\n", "mesh = TetMesh.LoadAbaqus('../meshes/axon.inp', scale=1e-6)\n", "\n", "with mesh:\n", " facetris = TriList([tri for tri in mesh.tris if tri.center.z == mesh.bbox.min.z])\n", " injverts = facetris.verts\n", "\n", " memb_tris = mesh.surface - facetris\n", "\n", " # The points along (z) axis at which to record potential\n", " pot_pos = np.arange(mesh.bbox.min.z, mesh.bbox.max.z, 10e-6)\n", " pot_tet = TetList(mesh.tets[0, 0, z] for z in pot_pos)\n", "\n", " cyto = Compartment.Create(mesh.tets)\n", " patch = Patch.Create(memb_tris, cyto, None, ssys)\n", "\n", " # Create the membrane across which the potential will be solved\n", " membrane = Membrane.Create([patch])\n", "\n", "#########################\n", "# Simulation setup\n", "#########################\n", "\n", "rng = RNG('mt19937', 512, 1234)\n", "\n", "partition = LinearMeshPartition(mesh, 1, 1, MPI.nhosts)\n", "\n", "sim = Simulation('TetOpSplit', model, mesh, rng, True, partition)\n", "\n", "rs = ResultSelector(sim)\n", "\n", "NaCurrs = rs.TRIS(memb_tris).VGNaC_I.I\n", "KCurrs = rs.TRIS(memb_tris).VGKC_I.I\n", "CellPot = rs.TETS(pot_tet).V\n", "\n", "NaCurrs.metaData['trizpos'] = [tri.center.z for tri in memb_tris]\n", "KCurrs.metaData['trizpos'] = [tri.center.z for tri in memb_tris]\n", "\n", "NaCurrs.metaData['triarea'] = [tri.Area for tri in memb_tris]\n", "KCurrs.metaData['triarea'] = [tri.Area for tri in memb_tris]\n", "\n", "CellPot.metaData['tetzpos'] = pot_pos\n", "\n", "sim.toSave(NaCurrs, KCurrs, CellPot, dt=DT_sim)\n", "\n", "#########################\n", "# Run simulation\n", "#########################\n", "\n", "with HDF5Handler('Efield_MPI') as hdf:\n", " sim.toDB(hdf, 'TetOpSplitSim')\n", "\n", " sim.newRun()\n", " \n", " # Inject channels\n", " surfarea = sim.patch.Area\n", " \n", " for state in VGNaC:\n", " prop = Na_facs[state.Count(Na_ha)][state.Count(Na_mo)]\n", " sim.patch.VGNaC[state].Count = Na_ro * surfarea * prop\n", " \n", " for state in VGKC:\n", " prop = K_facs[state.Count(Ko)]\n", " sim.patch.VGKC[state].Count = K_ro * surfarea * prop\n", " \n", " sim.patch.Leak[lsus].Count = L_ro * surfarea\n", " \n", " # Set dt for membrane potential calculation to 0.01ms\n", " sim.EfieldDT = 1.0e-5\n", " \n", " # Initialize potential to -65mV\n", " sim.membrane.Potential = -65e-3\n", " \n", " # Set capacitance of the membrane to 1 uF/cm^2 = 0.01 F/m^2\n", " sim.membrane.Capac = 1.0e-2\n", " \n", " # Set resistivity of the conduction volume to 100 ohm.cm = 1 ohm.meter\n", " sim.membrane.VolRes = 1.0\n", " \n", " # Set the current clamp\n", " sim.VERTS(injverts).IClamp = Iclamp/len(injverts)\n", " \n", " # Run the simulation\n", " sim.run(ENDT)\n", "```\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }