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 ]
