On Apr 12, 2010, at 9:45 AM, Igor wrote:
> Could you tell me if it is correct to break the problem into
> subproblems in the following way:
>
> "basic" problem:
> eq = TransTerm == DiffTerm_x - ConvTerm_x + ConvTermA_y - ConvTermB_y + Source
> ..
> while time <= timeStop:
> eq.solve(var=phi,...)
>
> "broken" problem:
> eqX = TransTerm == DiffTerm_x - ConvTerm_x + Source
> eqAY = TransTerm == ConvTermA_y
> eqBY = TransTerm == -ConvTermB_y
> ..
> while time <= timeStop:
> eqX.solve(var=phi,...)
> eqAY.solve(var=phi,...)
> eqBY.solve(var=phi,...)
>
> I tested it and the results agreed up to 4th digit after coma.
My gut reaction is that this is not valid (and that it certainly should be
"swept" to convergence so that all the terms see each other), but Wheeler is
probably better able to determine that. Also, even if this works, operator
splitting usually requires that you take fractional time steps, which you don't
seem to be doing.
> New fipy lets me do multiplication directly in "basic" problem:
> eq = TransTerm == DiffTerm_x - ConvTerm_x + ConvTermA_y*f(param) -
> ConvTermB_y + Source
>
> so I tried it and checked the results with "broken" problem results.
> Surprisingly, they didn't agree! So, I questioned the way the problem
> is broken!
I don't believe this has anything to do with your "broken" problem, as this
more minimal script demonstrates:
{{{
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() == 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 ]
}}}
Something certainly appears to be wrong with "term multiplication". Dan Wheeler
wrote that and he's out until next week. When he gets back, we will try to sort
out what is going on here.
I've filed http://matforge.org/fipy/ticket/291 to track this.
