Thanks Jon! So I guess as long as I stay under 4th order I'll be fine then.
On Wed, Jan 18, 2012 at 7:45 AM, Jonathan Guyer <[email protected]> wrote: > Thanks for your report. This example of solving the biharmonic equation > definitely has some issues with stability. I've filed a ticket about it ( > http://matforge.org/fipy/ticket/423), but I'd say we're unlikely to make > fixing it a priority unless you need it for your research. If you're just > exploring the different FiPy examples, I'd recommend you just move on to > the next one. This equation is a bit strange and it doesn't come up in the > sorts of problems that we typically solve. > > On Jan 16, 2012, at 7:58 PM, Yun Tao wrote: > > > Hi all, > > > > I've been working through the examples and, since my research focus is > on the dynamics before the solution reaches steady-state (if at all), I've > been re-coding everything with the transient term not set to zero. However, > on the 4th order diffusion example ( > http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.nthOrder.input4thOrder1D.html), > the program kept stalling/crashing everytime it's run (see attachment for > screenshots under different target residual values). This is after I > implemented continuous-sweeping-until-target-residual-is-met method that > worked elsewhere. What I did is reproduced below: > > > > ---------------------------------------------------------------- > > # mesh config > > L = 100. > > nx = 500. > > dx = L / nx > > mesh = Grid1D(nx=nx, dx=dx) > > > > > > # temporal config > > timeStepDuration = 0.5 > > steps = 100 > > > > > > # boundary conditions > > alpha1 = 2. > > alpha2 = 1. > > alpha3 = 4. > > alpha4 = -3. > > > > BCs = (FixedValue(faces=mesh.getFacesLeft(), value = alpha1), > > FixedFlux(faces=mesh.getFacesRight(), value = alpha2), > > NthOrderBoundaryCondition(faces=mesh.getFacesLeft(), value = > alpha3, order = 2), > > NthOrderBoundaryCondition(faces=mesh.getFacesRight(), value > = alpha4, order = 3)) > > > > # Variable > > phiS = CellVariable(name="transient", > > mesh=mesh, > > value=0., > > hasOld=1) > > > > # equation > > eqS = TransientTerm() == DiffusionTerm(coeff = (1,1)) > > > > > > # analytical solution > > analytical = CellVariable(mesh=mesh, name='analytical value') > > > > x = mesh.getCellCenters()[0] > > analytical.setValue(alpha4 / 6. * x**3 + alpha3 / 2. * x**2 + \ > > (alpha2 - alpha4/2. * L**2 - alpha3*L) * x + alpha1) > > > > # viewer config > > if __name__ == '__main__': > > viewer = Viewer(vars=(phiS,analytical)) > > > > # sweep until accurate > > for step in range(steps): > > > > phiS.updateOld() > > print 'step: ', step, '\n', 'phiS: ', phiS, '\n', 'phiS.old: ', > phiS.old() > > > > res = 1e+10 > > while res > 1e-6: > > res = eqS.sweep(var=phiS, > > boundaryConditions=BCs, > > dt = timeStepDuration) > > > > if __name__ == '__main__': > > viewer.plot() > > > > > > > ----------------------------------------------------------------------------- > > > > > > Is there's some intrinsic limitation to fipy's ability to handle these > systems of pdes or is my code simply flawed? Thank you for your generous > help. > > > > > > - Yun > > > > > > > > > > > > > > > > > > -- > > Yun Tao > > > > Graduate Group of Ecology Doctoral Candidate > > Department of Environmental Science and Policy > > Center for Population Biology > > University of California, Davis > > > > > > <nth_diff_res1e-6.jpg><nth_diff_res1e-2.jpg><ATT00001..c> > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Yun Tao Graduate Group of Ecology Doctoral Candidate Department of Environmental Science and Policy Center for Population Biology University of California, Davis
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
