Re: Integral term
On Thu, Jul 12, 2018 at 9:59 AM, Pavel Aleynikov wrote: > How should such equation be arranged? fp.numerix.cumsum seems to return a numpy array so won't automatically update when f1 changes. Probably best to update in the loop like this, import fipy as fp mesh = fp.Grid1D(dx=0.002, nx=100) mesh = mesh + [0.002] f1 = fp.CellVariable(name = "solution",mesh = mesh,value = 1.,hasOld=True) x = mesh.x M1 = lambda f1: fp.numerix.cumsum(f1*mesh.cellVolumes) M2 = lambda f1: fp.numerix.cumsum(x**2*f1*mesh.cellVolumes) Conv_pf = fp.CellVariable(mesh=mesh, rank=1) Diff_pf = fp.CellVariable(mesh=mesh) eq_steady1 = fp.DiffusionTerm(coeff=Diff_pf) + fp.PowerLawConvectionTerm(coeff=Conv_pf) f1.constrain(1,mesh.facesLeft) for sweep in range(10): Conv_pf[:] = M1(f1) / x**2 Diff_pf[:] = M2(f1) / x**3 eq_steady1.sweep(var=f1) print f1 f1.updateOld() -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Integral term
Hello, I would like to find a steady-state solution of an integro-differential equation with a spacial variable in the integration limit. So far it does not work. The minimal working example below gives wrong result. I understand that M1 and M2 are not updated in sweeps at all. If I set up M1 and M2 initially with the correct solution then the f1 solution is also correct. if I plug the integration directly to the term "fp.DiffusionTerm(coeff=fp.numerix.cumsum(f1*mesh.cellVolumes))” the error is: "IndexError: diffusion coefficent tensor is not an appropriate shape for this mesh" How should such equation be arranged? --- import fipy as fp mesh = fp.Grid1D(dx=0.002, nx=100) mesh = mesh + [0.002] f1 = fp.CellVariable(name = "solution",mesh = mesh,value = 1.,hasOld=True) x = mesh.x M1 = fp.numerix.cumsum(f1*mesh.cellVolumes) M2 = fp.numerix.cumsum(x**2*f1*mesh.cellVolumes) Conv_pf = [[1.0]] * (M1/x**2) Diff_pf = (M2/x**3) eq_steady1 = fp.DiffusionTerm(coeff=Diff_pf) + fp.PowerLawConvectionTerm(coeff=Conv_pf) f1.constrain(1,mesh.facesLeft) for sweep in range(10): eq_steady1.eq.sweep(var=f1) f1.updateOld() Thank you, Pavel ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
