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 ]
CylindricalGrid1D mesh volumes
Hi, How are "mesh.cellVolumes" defined in fp.CylindricalGrid1D case? The surface of a circle grid is not equal pi*r^2. import fipy as fp L= 1.; nx = 1000 mesh = fp.CylindricalGrid1D(dx=L/nx, nx=nx) print(mesh.cellVolumes.sum()) >> 0.5 Why not pi? How should I integrate Variables on such a grid? (2*pi*Var*mesh.cellVolumes).sum()? Thank you, Pavel PS. The solution of a diffusion equation on the CylindricalGrid1D seems fine, i.e. it converges to Bessel J_0 eigenvalue fairly well. PPS. I've accidentally sent this message from an unregistered email first and the cancelation link in the autoreply does not work. So, sorry if this is the duplicate. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: CylindricalGrid1D mesh volumes
Ok. I understand. I see the logic of saving on some operations. > Argulably, these quantities should be calculated in some > sort of ratio so the absolute value of the volume won't matter. Unless you have a source with a defined integral absolute value. But I agree that in general "theta = 1" is not a problem, as long as it's clear from the documentation. -- Pavel > On 15 Jan 2020, at 17:34, Daniel Wheeler wrote: > > On Wed, Jan 15, 2020 at 4:21 AM Pavel Aleynikov > wrote: >> >> Hi, >> >> How are "mesh.cellVolumes" defined in fp.CylindricalGrid1D case? >> The surface of a circle grid is not equal pi*r^2. >> >> import fipy as fp >> L= 1.; nx = 1000 >> mesh = fp.CylindricalGrid1D(dx=L/nx, nx=nx) >> print(mesh.cellVolumes.sum()) >>>> 0.5 >> >> Why not pi? > > I'm not sure. It's an arbitrary choice though. The angle was chosen as > "1" rather than "2 * pi". The volume of an element is "theta * r * dr" > where "theta" is the > angle, "r" are the cell centers and "dr" is the cell spacing. It's > possible that by choosing "theta=1", then the "theta" can be omitted > saving an extra operation. > >> How should I integrate Variables on such a grid? >> (2*pi*Var*mesh.cellVolumes).sum()? > > Makes sense. Argulably, these quantities should be calculated in some > sort of ratio so the absolute value of the volume won't matter. > > -- > Daniel Wheeler > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
