I'll look into this in detail when I return next week. Cheers. On Wed, Apr 14, 2010 at 2:55 PM, Jonathan Guyer <[email protected]> wrote: > > > 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. > >
-- Daniel Wheeler
