PyURDME includes a deterministic solver (PDE).
In [1]:
%matplotlib inline
from pyurdme.pde_solver import PDESolver
import pyurdme
import time
import scipy.fftpack
import numpy
import dolfin
import mshr
In [2]:
class RightEdge(dolfin.SubDomain): # Define the right subdomain
def inside(self, x, on_boundary):
return on_boundary and x[0] > 4.9
class LeftEdge(dolfin.SubDomain): # Define the left subdomain
def inside(self, x, on_boundary):
return on_boundary and x[0] < -4.9
In [3]:
class CylinderDemo(pyurdme.URDMEModel):
def __init__(self, model_name="CylinderDemo"):
pyurdme.URDMEModel.__init__(self, model_name)
# Define Species
A = pyurdme.Species(name="A", diffusion_constant=1.0)
B = pyurdme.Species(name="B", diffusion_constant=1.0)
self.add_species([A, B])
# Define Geometry
cylinder = mshr.Cylinder(dolfin.Point(5, 0, 0), dolfin.Point(-5, 0, 0), 1.0, 1.0)
self.mesh = pyurdme.URDMEMesh(mesh=mshr.generate_mesh(cylinder, 32))
# Define Subdomains
self.add_subdomain(RightEdge(), 2)
self.add_subdomain(LeftEdge(), 3)
# This section is a fix to ensure the creation rates are equal. Since
# The creation reactions are zero-th order, They are proportional to volume.
# Due to the discritization the sum volume of the voxels in subdomain 1 may
# not be the same as the sum volume of the voxels in subdomain 2. We adjust
# the parameters so that the net creation rate the same for both.
data = self.get_solver_datastructure()
vol = data['vol']
sd = data['sd']
right = numpy.sum(vol[sd == 2])
left = numpy.sum(vol[sd == 3])
print "Total Subdomain Volume: Right(A)={0} Left(B)={1}".format(left, right)
k_creat = 100
k_creatA_val = k_creat*left/right
k_creatB_val = k_creat
print "Create params: A={0} B={1}".format(k_creatA_val, k_creatB_val)
print "Net create rate: A={0} B={1}".format(k_creatA_val*left, k_creatB_val*right)
# Define Parameters
k_react = pyurdme.Parameter(name="k_react", expression=1.0)
k_create1 = pyurdme.Parameter(name="k_create1", expression=k_creatA_val)
k_create2 = pyurdme.Parameter(name="k_create2", expression=k_creatB_val)
self.add_parameter([k_react, k_create1, k_create2])
# Define Reactions
R1 = pyurdme.Reaction(name="R1",reactants=None, products={A:1}, rate=k_create1, restrict_to=2)
R2 = pyurdme.Reaction(name="R2",reactants=None, products={B:1}, rate=k_create2, restrict_to=3)
R3 = pyurdme.Reaction(name="R3",reactants={A:1, B:1}, products=None, rate=k_react)
self.add_reaction([R1, R2, R3])
# Define simulation timespan
self.timespan(range(200))
In [4]:
model = CylinderDemo()
sol = PDESolver(model)
print "Attempting to compile"
t1 = time.time()
sol.compile()
print "Compilation complete in {0}s".format(time.time()-t1)
print "Beginning simulation"
t1 = time.time()
result = sol.run()
print "Simulation complete in {0}s".format(time.time()-t1)
In [7]:
print "Plotting solution"
# Plot of the time-average spatial concentration.
x_vals = model.mesh.coordinates()[:, 0]
A_vals = result.get_species("A",-1, concentration=True)
B_vals = result.get_species("B",-1, concentration=True)
plt.plot(x_vals,A_vals,'.r',x_vals,B_vals,'.b')
plt.legend(['A', 'B'])
Out[7]:
In [8]:
x_vals = model.mesh.coordinates()[:, 0]
#A_vals = numpy.sum(result.get_species("A", concentration=True), axis=0)
#B_vals = numpy.sum(result.get_species("B", concentration=True), axis=0)
A_vals = result.get_species("A",-1, concentration=True)
B_vals = result.get_species("B",-1, concentration=True)
plt.figure(figsize=(12,6), dpi=100)
num_bins = 25
plt.hist(x_vals, bins=num_bins, weights=B_vals, facecolor='red', alpha=0.5)
plt.hist(x_vals, bins=num_bins, weights=A_vals, facecolor='blue', alpha=0.5)
plt.legend(['B','A'])
plt.xlabel('x-coordinate')
plt.ylabel('concentration')
plt.title('Concentration binned by x-coordinate')
Out[8]:
In [ ]:
In [ ]:
In [11]:
sd = model.get_solver_datastructure()
print sd.keys()
In [13]:
sd['data'].shape
Out[13]:
In [17]:
Data = numpy.zeros((2,1063))
In [18]:
Data[:,155]
Out[18]:
In [ ]:
In [4]:
print "Plotting solution"
# Plot of the time-average spatial concentration.
x_vals = model.mesh.coordinates()[:, 0]
A_vals = numpy.mean(result.get_species("A", concentration=True), axis=0)
B_vals = numpy.mean(result.get_species("B", concentration=True), axis=0)
plt.plot(x_vals,A_vals,'.r',x_vals,B_vals,'.b')
plt.legend(['A', 'B'])
Out[4]:
In [5]:
sol = PDESolver(model, report_level=0)
print "Attempting to compile"
sol.compile()
N = 10
run_times = []
print "Beginning simulation:"
for n in range(N):
sys.stdout.write(str(n))
sys.stdout.flush()
t1 = time.time()
result = sol.run()
elapsed_time = time.time()-t1
run_times.append(elapsed_time)
print "Simulation of {0} trajectories: avg={1} std={2}".format(N, numpy.mean(run_times), numpy.std(run_times))
In [6]:
from pyurdme.nsmsolver import NSMSolver
nsm = NSMSolver(model, report_level=0)
print "NSM: Attempting to compile"
nsm.compile()
N = 10
nsm_times = []
print "NSM: Beginning simulation:",
for n in range(N):
sys.stdout.write(str(n))
sys.stdout.flush()
t1 = time.time()
result = nsm.run()
elapsed_time = time.time()-t1
nsm_times.append(elapsed_time)
print ""
print "NSM: Simulation of {0} trajectories: avg={1} std={2}".format(N, numpy.mean(nsm_times), numpy.std(nsm_times))
In [ ]:
In [ ]: