usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
504 stars 148 forks source link

term multiplication changes result #195

Closed guyer closed 9 years ago

guyer commented 9 years ago

As highlighted in http://thread.gmane.org/gmane.comp.python.fipy/1711, multiplying terms by 1 seems to change their behavior. E.g.,

from fipy import *

mesh = Grid2D(nx=2, ny=2)

x, y = mesh.getCellCenters()

phi1 = CellVariable(mesh=mesh)
phi2 = CellVariable(mesh=mesh)
phi3 = CellVariable(mesh=mesh)

BCs = (FixedValue(faces=mesh.getFacesLeft(), value=0),
FixedValue(faces=mesh.getFacesRight(), value=0))

eq1 = TransientTerm() h2. DiffusionTerm(coeff=1.) + PowerLawConvectionTerm(coeff=[[1], [0]]) + 1.0
eq2 = TransientTerm() DiffusionTerm(coeff=1.) + PowerLawConvectionTerm(coeff=[[1], [0]]) * 1.0 + 1.0
eq3 = TransientTerm() == DiffusionTerm(coeff=1.) * 1.0 + PowerLawConvectionTerm(coeff=[[1], [0]]) + 1.0

eq1.sweep(var=phi1, boundaryConditions=BCs)
eq2.sweep(var=phi2, boundaryConditions=BCs)
eq3.sweep(var=phi3, boundaryConditions=BCs)

print "phi1:", phi1
print "phi2:", phi2
print "phi3:", phi3

returns

phi1: [ 0.35395662  0.29216836  0.35395662  0.29216836]
phi2: [ 0.33333333  0.33333333  0.33333333  0.33333333]
phi3: [ 0.75  0.5   0.75  0.5 ]

Imported from trac ticket #291, created by guyer on 04-14-2010 at 14:52, last modified: 11-17-2011 at 14:34

fipymigrate commented 9 years ago

Also encountered this issue. Solving: A.equation = TransientTerm() h2. DiffusionTerm() - A.getFaceGrad().getDivergence() returned the expected result of no change (transient = 0).

Solving what should be an identical equation: A.equation = TransientTerm() 1*DiffusionTerm() - A.getFaceGrad().getDivergence() and did not get an identical result.

I'm happy to provide the full script and results if they will be helpful. In summary, I started with a step function (1's, then 0's). At the step the second equation produces these results (compared with the expected result of no change):

Expected Observed 1 1 1 0.999999989 1 1.000001074 1 0.999894781 1 1.010310363 0 -0.010310363

Note also that similar behavior occurred with multiples other than 1* (2, B, etc).

Trac comment by gathrw@rpi.edu on 05-09-2010 at 12:42

wd15 commented 9 years ago

Check that this is fixed???

Trac comment by wd15 on 04-27-2011 at 16:13

wd15 commented 9 years ago

This seems to have been fixed! No idea how or when this occurred. Made some minor modifications as the script above didn't run against trunk. See r4959.

Trac comment by wd15 on 11-17-2011 at 14:34