27
27
28
28
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
29
29
30
+ from __future__ import print_function
30
31
import math
31
32
import numpy
32
33
import pylab
50
51
INT = 0.101
51
52
52
53
# The number of molecules to be injected into the centre
53
- NINJECT = 10000
54
+ NINJECT = 10000
54
55
55
56
# The number of tetrahedral elements to sample data from.
56
- SAMPLE = 2000
57
+ SAMPLE = 2000
57
58
58
59
# The diffusion constant for our diffusing species (m^2/s)
59
60
DCST = 20.0e-12
70
71
71
72
def printtime (end_time ):
72
73
73
- totsecs = int (end_time - beg_time )
74
- sec = totsecs % 60
75
- totmin = totsecs / 60
76
- min = totmin % 60
77
- hours = totmin / 60
78
-
79
- print 'Simulation time: %d h, %d min, %d sec' % (hours , min , sec )
74
+ totsecs = int (end_time - beg_time )
75
+ sec = totsecs % 60
76
+ totmin = totsecs / 60
77
+ min = totmin % 60
78
+ hours = totmin / 60
79
+
80
+ print ( 'Simulation time: %d h, %d min, %d sec' % (hours , min , sec ) )
80
81
81
82
########################################################################
82
83
@@ -92,17 +93,17 @@ def gen_model():
92
93
93
94
def gen_geom ():
94
95
95
- print "Loading mesh..."
96
+ print ( "Loading mesh..." )
96
97
mesh = smeshio .loadMesh ('meshes/sphere_rad10_11Ktets' )[0 ]
97
- print "Mesh Loaded"
98
+ print ( "Mesh Loaded" )
98
99
99
- # Find the total number of tetrahedrons in the mesh
100
+ # Find the total number of tetrahedrons in the mesh
100
101
ntets = mesh .countTets ()
101
102
# Create a compartment containing all tetrahedron
102
103
comp = stetmesh .TmComp ('cyto' , mesh , range (ntets ))
103
104
comp .addVolsys ('cytosolv' )
104
105
105
- print "Finding tetrahedron samples..."
106
+ print ( "Finding tetrahedron samples..." )
106
107
# Fetch the central tetrahedron index and store:
107
108
ctetidx = mesh .findTetByPoint ([0.0 , 0.0 , 0.0 ])
108
109
tetidxs [0 ] = ctetidx
@@ -155,7 +156,7 @@ def gen_geom():
155
156
# Store the radial distance (in microns):
156
157
tetrads [i ] = r * 1.0e6
157
158
158
- print "Tetrahedron samples found"
159
+ print ( "Tetrahedron samples found" )
159
160
160
161
return mesh
161
162
@@ -181,57 +182,57 @@ def gen_geom():
181
182
182
183
# Run NITER number of iterations:
183
184
for i in range (NITER ):
184
- sim .reset ()
185
- print "Running iteration" , i
185
+ sim .reset ()
186
+ print ( "Running iteration" , i )
186
187
# Inject all molecules into the central tet:
187
- sim .setTetCount (ctetidx , 'A' , NINJECT )
188
- for j in range (ntpnts ):
189
- sim .run (tpnts [j ])
188
+ sim .setTetCount (ctetidx , 'A' , NINJECT )
189
+ for j in range (ntpnts ):
190
+ sim .run (tpnts [j ])
190
191
# Loop over the tetrahedrons we are saving data for
191
- for k in range (SAMPLE ):
192
+ for k in range (SAMPLE ):
192
193
# Save the concentration in the tetrahedron, in uM
193
- res [i , j , k ] = sim .getTetConc (int (tetidxs [k ]), 'A' )* 1.0e6
194
- printtime (time .time ())
194
+ res [i , j , k ] = sim .getTetConc (int (tetidxs [k ]), 'A' )* 1.0e6
195
+ printtime (time .time ())
195
196
196
197
res_mean = numpy .mean (res , axis = 0 )
197
198
198
199
########################################################################
199
200
200
201
def plotres (tidx ):
201
- if (tidx >= INT / DT ):
202
- print "Time index is out of range."
203
- return
204
-
205
- pylab .scatter (tetrads , res_mean [tidx ], s = 2 )
206
- pylab .xlabel ('Radial distance of tetrahedron ($\mu$m)' )
207
- pylab .ylabel ('Concentration in tetrahedron ($\mu$M)' )
208
- t = tpnts [tidx ]
209
- pylab .title ('Unbounded diffusion. Time: ' + str (t ) + 's' )
210
- plotanlyt (t )
211
- pylab .xlim (0.0 , 10.0 )
212
- pylab .ylim (0.0 )
213
- pylab .show ()
202
+ if (tidx >= INT / DT ):
203
+ print ( "Time index is out of range." )
204
+ return
205
+
206
+ pylab .scatter (tetrads , res_mean [tidx ], s = 2 )
207
+ pylab .xlabel ('Radial distance of tetrahedron ($\mu$m)' )
208
+ pylab .ylabel ('Concentration in tetrahedron ($\mu$M)' )
209
+ t = tpnts [tidx ]
210
+ pylab .title ('Unbounded diffusion. Time: ' + str (t ) + 's' )
211
+ plotanlyt (t )
212
+ pylab .xlim (0.0 , 10.0 )
213
+ pylab .ylim (0.0 )
214
+ pylab .show ()
214
215
215
216
########################################################################
216
217
217
- def plotanlyt (t ):
218
- segs = 100
219
- anlytconc = numpy .zeros ((segs ))
220
- radialds = numpy .zeros ((segs ))
221
- maxrad = 0.0
222
- for i in tetrads :
223
- if (i > maxrad ): maxrad = i
224
- maxrad *= 1e-6
225
- intervals = maxrad / segs
226
- rad = 0.0
227
- for i in range ((segs )):
228
- # Find the conc from analytical solution, and convert to mol/L
229
- anlytconc [i ]= 1.0e3 * (1 / 6.022e23 )* \
230
- ((NINJECT / (math .pow ((4 * math .pi * DCST * t ),1.5 )))* \
231
- (math .exp ((- 1.0 * (rad * rad ))/ (4 * DCST * t ))))
232
- radialds [i ] = rad * 1e6
233
- rad += intervals
234
- pylab .plot (radialds , anlytconc , color = 'red' )
218
+ def plotanlyt (t ):
219
+ segs = 100
220
+ anlytconc = numpy .zeros ((segs ))
221
+ radialds = numpy .zeros ((segs ))
222
+ maxrad = 0.0
223
+ for i in tetrads :
224
+ if (i > maxrad ): maxrad = i
225
+ maxrad *= 1e-6
226
+ intervals = maxrad / segs
227
+ rad = 0.0
228
+ for i in range ((segs )):
229
+ # Find the conc from analytical solution, and convert to mol/L
230
+ anlytconc [i ]= 1.0e3 * (1 / 6.022e23 )* \
231
+ ((NINJECT / (math .pow ((4 * math .pi * DCST * t ),1.5 )))* \
232
+ (math .exp ((- 1.0 * (rad * rad ))/ (4 * DCST * t ))))
233
+ radialds [i ] = rad * 1e6
234
+ rad += intervals
235
+ pylab .plot (radialds , anlytconc , color = 'red' )
235
236
236
237
########################################################################
237
238
0 commit comments