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 <wes...@bu.edu> 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 <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 >> <erik.py> > > > -- Daniel Wheeler
