Hi Erik, I did a few more efficient improvements to FiPy on a branch.
<http://matforge.org/fipy/browser/branches/efficiency>
I have also attached and updated script. The new script (500x500) and
gives the following timings with OLD being old script + trunk and NEW
being new script + branches/efficiency.
OLD (s)
NEW (s)
--pysparse 175
66
--trilinos 691
119
mpirun -np 4 254
48
mpirun -np 12 106
20
(on slow network)
This is a lot faster now. We haven't done a lot of our own research
with version 2.1 of FiPy yet and hadn't been through any serious
optimization cycles. I think there are still some major efficiencies
particularly in populating the trilinos matrix, which can be improved
on, but that will take more time. I will try and merge these changes
back to version-2_1 ASAP and maybe do another release.
Cheers
On Thu, Oct 7, 2010 at 12:55 PM, Daniel Wheeler
<[email protected]> wrote:
>
> Erik, FYI <http://matforge.org/fipy/browser/branches/erik> has moved
> to <http://matforge.org/fipy/browser/branches/efficiency> and is an
> offshoot of version-2_1 rather than trunk. Cheers
>
> On Wed, Oct 6, 2010 at 5:41 PM, Erik Sherwood <[email protected]> wrote:
>>
>> Thanks very much for looking into this, and for taking the time to edit
>> FiPy! I'm heading to a conference and won't be able to check out the new
>> FiPy for several days, but I will try it out once I'm back and let you know
>> what difference it makes on my system, and how it performs in the "adaptive"
>> time-stepping loop that I implemented.
>>
>> Best,
>>
>> Erik
>>
>> On Oct 5, 2010, at 3:09 PM, Daniel Wheeler wrote:
>>
>>> 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 <[email protected]> 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 <[email protected]> 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
>>> <erik.py>
>>
>>
>>
>
>
>
> --
> Daniel Wheeler
>
>
>
--
Daniel Wheeler
import os
os.system('echo ${PYTHONPATH}')
from fipy import *
from time import clock
# Setting up the mesh
nx = 500
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)
solverU = LinearPCGSolver(tolerance=1e-8, precon=None)
solverV = LinearPCGSolver(tolerance=1e-8, precon=None)
solverW = LinearPCGSolver(tolerance=1e-8, precon=None)
eqns = ((uu, eqXuu, solverU), (vv, eqXvv, solverV), (ww, eqXww, solverW))
## 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 = 10
desiredResidual = 1e-6
# Loop through timesteps
tic = clock()
if '--profile' in sys.argv:
profileON = True
print 'profiling...'
else:
profileON = False
for step in range(steps):
if profileON:
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, solver 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
print 'residual',residual
# 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()
if profileON:
profile.stop()
print "Finished in %.3f secs."%(toc-tic)
