Hi Erik, I've attached a script that seems to run faster. It took
57.60s to do 200 steps as opposed to your original script, which took
183.50. In order to get some of this speed up, I had to hack the FiPy
code a little. The hack is on a branch called erik
<http://matforge.org/fipy/browser/branches/erik>
<http://matforge.org/fipy/changeset?new=3893%40branches%2Ferik&old=3892%40branches%2Ferik>
I'll review this and get it back into trunk ASAP.
A few other things:
* turned off the viewers
* reduced the solver tolerance (1e-7)
* using implicit source terms that don't seem to cause any issues at
least up to 200 time steps.
I can run the script out to 300 steps without any evidence of stiffness.
Cheers
On Sat, Oct 2, 2010 at 2:44 PM, Erik Sherwood <wes...@bu.edu> wrote:
>
> Thanks for taking a look at my code.
>
> For setting the desiredResidual, I just mimicked what I saw in
> examples.diffusion.mesh1D (p.81 of the manual), which has a similar form for
> the sweeps loop.
>
> Using ImplicitSourceTerms for Suu, Svv, Sww, the time to take 150 steps
> (dt=2.0) is cut from 141 secs to 99 secs, which is great. But if I go
> further, right around step 165 I get the warning
>
> mimura2.py:96: StagnatedSolverWarning: The solver stagnated. Iterations: 1.
> Relative error: nan
> residual = eqn.sweep(var=var, dt=dt)
>
> And from then on the simulation produces nothing of value. In this scenario,
> the FiPy makes 4-5 sweeps per step on Suu and Svv, and 1 sweep for Sww, up
> until the 130th step. From there the number of sweeps for Suu climbs, until
> it is about 70 or so per step, until the solver stagnates.
>
> If I try the same run parameters with Suu, Svv, Sww not as
> ImplicitSourceTerms, then I can go for up to about 220 steps before the
> simulation gets stuck (that is, it takes more than 300 sweeps without
> achieving the desiredResidual). In this case, FiPy makes 5-6 sweeps per step
> on Suu and Svv, and 2 sweeps for Sww, up to the 216th step. From there, the
> number of sweeps for Suu climbs again.
>
> I suppose one thing to do would be to keep track of how many sweeps are
> needed for each step, and once that gets to be too high, reduce dt. That
> should work whether I'm using ImplicitSourceTerms or not.
>
> Are there other things to try?
>
> Thanks,
>
> Erik
>
> On Sep 30, 2010, at 2:00 PM, Jonathan Guyer wrote:
>
>>
>> I'm running a workshop right now, so won't have time to look at this for a
>> day or two and it sounds like Daniel is on it.
>>
>> One thing that occurs to me is that you should be able to use an
>> ImplicitSourceTerm for Suu and Svv, which might reduce the number of sweeps.
>>
>> How many sweeps are you doing?
>>
>> How did you pick desiredResidual? Generally, you want to compare the
>> residual to the initial residual, rather than some arbitrary number.
>>
>> On Sep 30, 2010, at 12:55 PM, Erik Sherwood wrote:
>>
>>>
>>> Thanks for taking a look at the code. I'm using FiPy 2.1.
>>>
>>> Erik
>>>
>>> On Sep 30, 2010, at 11:21 AM, Daniel Wheeler wrote:
>>>
>>>>
>>>> Eric, Thanks for your interest in FiPy. I'll benchmark your code, but
>>>> before I do that, can you tell me which version of FiPy you are using
>>>> (trunk or version-2_1?) so we are comparing apples with apples. Thanks
>>>>
>>>> On Wed, Sep 29, 2010 at 5:03 PM, Erik Sherwood <wes...@bu.edu> wrote:
>>>>>
>>>>> Hi all,
>>>>>
>>>>> I am quite new to FiPy, which I have been using to simulate a
>>>>> reaction-diffusion model of bacterial colony growth:
>>>>>
>>>>> # Simulation parameters
>>>>> d = 0.12
>>>>> a0 = 1.; a1 = 1./2400; a2 = 1./120;
>>>>> v0 = 0.071
>>>>>
>>>>> # Source terms
>>>>> Suu = uu*vv - uu/((1+uu/a1)*(1+vv/a2))
>>>>> Svv = -uu*vv
>>>>> Sww = uu/((1+uu/a1)*(1+vv/a2))
>>>>>
>>>>> # Equations
>>>>> eqXuu = TransientTerm() - ImplicitDiffusionTerm(coeff=d) - Suu
>>>>> eqXvv = TransientTerm() - ImplicitDiffusionTerm(coeff=1.) - Svv
>>>>> eqXww = TransientTerm() - Sww
>>>>>
>>>>> There are three PDEs, two of which have linear diffusion terms, and
>>>>> the
>>>>> equations are coupled via nonlinear reaction terms. I hope to be
>>>>> able to use
>>>>> Fipy on a fairly large mesh to explore the pattern formation
>>>>> possibilities
>>>>> of the system as certain parameters are varied, but the simulations
>>>>> run
>>>>> significantly slower than I'd expected, about 4-5 minutes per run
>>>>> (with
>>>>> --inline and the pysparse solver on a 2008 Mac Pro) on a 150x150
>>>>> square
>>>>> mesh; 600x600 is more like what I'd like to be able to use.
>>>>>
>>>>> Before using FiPy, I'd implemented the system in Matlab, just using
>>>>> finite
>>>>> differences with a 9 point stencil. I haven't gotten much, if any,
>>>>> speed up
>>>>> using FiPy, but I think I'm not using FiPy's full capabilities, or
>>>>> I've not
>>>>> implemented the model in FiPy as efficiently as I could have. The
>>>>> script I'm
>>>>> using is below. I suspect that the looping I'm doing for the sweeps
>>>>> of each
>>>>> variable in succession is not the right way to go about things.
>>>>> I've played
>>>>> around with sweeps and time step sizes, but haven't gotten much
>>>>> improvement.
>>>>>
>>>>> If the looping isn't the problem, should I try moving to trilinos?
>>>>> Or should
>>>>> I not expect particularly fast solution of this kind of system?
>>>>>
>>>>> Any suggestions for better performance would be welcome.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Erik
>>>>>
>>>>>
>>>>> from fipy import *
>>>>> from time import clock
>>>>>
>>>>> # Setting up the mesh
>>>>> nx = 150
>>>>> ny = nx
>>>>>
>>>>> dx = 1.0
>>>>> dy = dx
>>>>> L = dx*nx
>>>>>
>>>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
>>>>>
>>>>>
>>>>> # Simulation parameters
>>>>> d = 0.12
>>>>> a0 = 1.; a1 = 1./2400; a2 = 1./120;
>>>>> v0 = 0.071
>>>>>
>>>>> # Simulation variables
>>>>> uu = CellVariable(name="active bacteria",
>>>>> mesh=mesh,
>>>>> value=0.0, hasOld=1)
>>>>> ww = CellVariable(name="inactive bacteria",
>>>>> mesh=mesh,
>>>>> value=0.0, hasOld=1)
>>>>> vv = CellVariable(name="nutrient concentration",
>>>>> mesh=mesh,
>>>>> value=v0, hasOld=1)
>>>>>
>>>>> # Source terms
>>>>> Suu = uu*vv - uu/((1+uu/a1)*(1+vv/a2))
>>>>> Svv = -uu*vv
>>>>> Sww = uu/((1+uu/a1)*(1+vv/a2))
>>>>>
>>>>> # Equations
>>>>> eqXuu = TransientTerm() - ImplicitDiffusionTerm(coeff=d) - Suu
>>>>> #ImplicitSourceTerm(Suu)
>>>>> eqXvv = TransientTerm() - ImplicitDiffusionTerm(coeff=1.) - Svv
>>>>> #ImplicitSourceTerm(Svv)
>>>>> eqXww = TransientTerm() - Sww #ImplicitSourceTerm(Sww)
>>>>>
>>>>>
>>>>> # Initial conditions
>>>>> x,y = mesh.getCellCenters()
>>>>> radius = 10*dx
>>>>> C = (nx*dx/2, ny*dy/2)
>>>>> C =(0,0)
>>>>> uu.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.09,
>>>>> maximum=0.115))
>>>>> uu.setValue(0.0,where=(0.05*(x-C[0])**2 + (y-C[1])**2) > radius**2)
>>>>>
>>>>> eqns = ((uu, eqXuu), (vv, eqXvv), (ww, eqXww))
>>>>>
>>>>> if __name__ == '__main__':
>>>>> # Do plotting if executing script from commandline
>>>>> Uview = uu
>>>>> Uview.setName('Active')
>>>>> Uviewer = Viewer(vars=uu, datamin=0., datamax=0.1)
>>>>>
>>>>> Vview = vv
>>>>> Vview.setName('Nutrient')
>>>>> Vviewer = Viewer(vars=Vview, datamin=0., datamax=0.1)
>>>>>
>>>>> Wview = ww
>>>>> Wview.setName('Inactive')
>>>>> Wviewer = Viewer(vars=Wview, datamin=0., datamax=0.1)
>>>>>
>>>>> dt = 2.0
>>>>> steps = 1000
>>>>> desiredResidual = 1e-6
>>>>>
>>>>> # Loop through timesteps
>>>>> tic = clock()
>>>>> for step in range(steps):
>>>>> print "Step %d"%(step)
>>>>>
>>>>> # For each variable, equation pair, do sufficient sweeps to
>>>>> # get residual to acceptable level. When all three variables
>>>>> # have been taken care of, then we take the next time step
>>>>> for var, eqn in eqns:
>>>>> residual = 10
>>>>> var.updateOld()
>>>>> sweepcount = 0
>>>>> while residual > desiredResidual and sweepcount < 300:
>>>>> #print "Sweep %d"%(sweepcount+1)
>>>>> residual = eqn.sweep(var=var, dt=dt)
>>>>> sweepcount += 1
>>>>> # Bail if residual remains too large
>>>>> if residual > desiredResidual:
>>>>> print "Unable to reduce residual!"
>>>>> break
>>>>> if residual > desiredResidual:
>>>>> break
>>>>>
>>>>> # Update the plots
>>>>> if __name__ == '__main__':
>>>>> Uviewer.plot()
>>>>> Vviewer.plot()
>>>>> Wviewer.plot()
>>>>> #
>>>>> toc = clock()
>>>>> print "Finished in %.3f secs."%(toc-tic)
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Daniel Wheeler
>>>>
>>>>
>>>>
>>>
>>>
>>
>>
>>
>>
>
>
>
--
Daniel Wheeler
from fipy import *
from time import clock
# Setting up the mesh
nx = 150
ny = nx
dx = 1.0
dy = dx
L = dx*nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
# Simulation parameters
d = 0.12
a0 = 1.; a1 = 1./2400; a2 = 1./120;
v0 = 0.071
# Simulation variables
uu = CellVariable(name="active bacteria",
mesh=mesh,
value=0.0, hasOld=1)
ww = CellVariable(name="inactive bacteria",
mesh=mesh,
value=0.0, hasOld=1)
vv = CellVariable(name="nutrient concentration",
mesh=mesh,
value=v0, hasOld=1)
# Source terms
## Suu = uu*vv - uu/((1+uu/a1)*(1+vv/a2))
## Svv = -uu*vv
## Sww = uu/((1+uu/a1)*(1+vv/a2))
Suu = vv - 1/((1+uu/a1)*(1+vv/a2))
Svv = -uu
Sww = uu/((1+uu/a1)*(1+vv/a2))
# Equations
eqXuu = TransientTerm() - ImplicitDiffusionTerm(coeff=d) - ImplicitSourceTerm(Suu)
eqXvv = TransientTerm() - ImplicitDiffusionTerm(coeff=1.) - ImplicitSourceTerm(Svv)
eqXww = TransientTerm() - Sww #ImplicitSourceTerm(Sww)
# Initial conditions
x,y = mesh.getCellCenters()
radius = 10*dx
C = (nx*dx/2, ny*dy/2)
C =(0,0)
uu.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.09, maximum=0.115))
uu.setValue(0.0,where=(0.05*(x-C[0])**2 + (y-C[1])**2) > radius**2)
eqns = ((uu, eqXuu), (vv, eqXvv), (ww, eqXww))
solver = LinearPCGSolver(tolerance=1e-7)
if __name__ == '__main__':
# Do plotting if executing script from commandline
Uview = uu
Uview.setName('Active')
Uviewer = Viewer(vars=uu, datamin=0., datamax=0.1)
Vview = vv
Vview.setName('Nutrient')
Vviewer = Viewer(vars=Vview, datamin=0., datamax=0.1)
Wview = ww
Wview.setName('Inactive')
Wviewer = Viewer(vars=Wview, datamin=0., datamax=0.1)
dt = 2.0
steps = 300
desiredResidual = 1e-6
# Loop through timesteps
tic = clock()
for step in range(steps):
## if step == 1:
## from profiler import Profiler
## from profiler import calibrate_profiler
## fudge = calibrate_profiler(10000)
## profile = Profiler('profile', fudge=fudge)
print "Step %d"%(step)
# For each variable, equation pair, do sufficient sweeps to
# get residual to acceptable level. When all three variables
# have been taken care of, then we take the next time step
for var, eqn in eqns:
residual = 10
var.updateOld()
sweepcount = 0
while residual > desiredResidual and sweepcount < 300:
residual = eqn.sweep(var=var, dt=dt, solver=solver)
sweepcount += 1
# Bail if residual remains too large
print "Sweep %d"%(sweepcount+1)
if residual > desiredResidual:
print "Unable to reduce residual!"
break
if residual > desiredResidual:
break
# Update the plots
## if __name__ == '__main__':
## Uviewer.plot()
## Vviewer.plot()
## Wviewer.plot()
#
toc = clock()
##profile.stop()
print "Finished in %.3f secs."%(toc-tic)
