Test usage of the experimental PDE solver

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)


Total Subdomain Volume: Right(A)=0.338902497332 Left(B)=0.243929901352
Create params: A=138.934380514 B=100
Net create rate: A=47.0852085215 B=24.3929901352
Attempting to compile
Compilation complete in 259.804553986s
Beginning simulation
Simulation complete in 270.856409073s

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'])


Plotting solution
Out[7]:
<matplotlib.legend.Legend at 0x1134626d0>

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]:
<matplotlib.text.Text at 0x1141dfc50>

In [ ]:


In [ ]:


In [11]:
sd = model.get_solver_datastructure()
print sd.keys()


['D', 'G', 'vol', 'tspan', 'K', 'u0', 'dofvolumes', 'N', 'p', 'report', 'data', 'sd']

In [13]:
sd['data'].shape


Out[13]:
(1, 1063)

In [17]:
Data = numpy.zeros((2,1063))

In [18]:
Data[:,155]


Out[18]:
array([ 0.,  0.])

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'])


Plotting solution
Out[4]:
<matplotlib.legend.Legend at 0x110936650>

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))


Attempting to compile
Beginning simulation:
0123456789Simulation of 10 trajectories: avg=3.17652721405 std=0.457474561977

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))


NSM: Attempting to compile
NSM: Beginning simulation:0123456789 
NSM: Simulation of 10 trajectories: avg=3.61504948139 std=0.214400820178

In [ ]:


In [ ]: