On Nov 17, 2008, at 10:18 AM, asreeve wrote:
I've been struggling for quite some time trying to get FiPy to update
a rank 1 CellVariable term as it sweeps to a solution. I think I've
finally got it working in a reasonable way, but don't understand why
the original way I was doing this was so slow and failed to update.
The original, apparently incorrect way, involved creating a rank 1
cell variable of permeabilities:
kond=CellVariable(name='hyd.cond.',mesh=mesh,rank=1,value=1.e-5)
and multiplying this by a 'saturation phase' variables, with the
impact from the saturation states being different in the horizontal
(cellSatH) and vertical (cellSatV) directions.
cellSat=kond*[cellSatH,cellSatH,cellSatV]
This single calculation takes up most of the time in the simulation
and this variable never changes as cellSatH and cellSatV change.
The reason it doesn't change as the problem evolves is because FiPy
doesn't see cellSatH and cellSatV changing. From this statement, FiPy
composes a binOp of the CellVariable kond and a constant list
containing the values of cellSatH and cellSatV at the time of
definition. The list is not a FiPy Variable and so it never sees the
contents change.
I suspect that the reason it's slow is because it's using Python list
semantics to iterate through the values of
[cellSatH,cellSatH,cellSatV], instead fast NumPy array semantics.
You could define a Variable:
cellSatCoeff = Variable(value=[cellSatH0,cellSatH0,cellSatV0])
and then change the coefficient values by updating appropriate values
in the [0], [1], and [2] indices of cellSatCoeff, but the easiest
thing is probably more or less what you've already done:
kond * (cellSatH * [1, 1, 0] + cellSatV * [0, 0, 1])
The second approach isn't quite right, as the vertical kond term is
being doubled in fully saturated cases, but is easy to fix if I just
set the vertical kond term to half its real value.
I believe what I've done above should avoid this problem
In my many attempts to sort this out, I varied the solver I was using
and noted this simulation solves when I select the 'linearPCG' solver,
but never converges when I use 'linearGMRES', and returns NaNs when I
use the 'linearCGS' solver. I was surprised that the solver choice was
so important for such a simple problem. Would you have expected this
behavior? Does this indicate somethings wrong with pysparse, at least
on my computer? What do I need to do to try the scipy sparse solvers
to see how they perform?
I'll let Daniel field your solver questions. I see the same results
you do, but note that PCG is the default.
I've attached the code I'm using. The print statements in the sweeps
prints out values for one vertical column of cells.
Just so you know, this:
print satState()[where((X<25.1)&(X>25)&(Y>.12)&(Y<.13))]
print cellSatH()[where((X<25.1)&(X>25)&(Y>.12)&(Y<.13))]
can be simplified to this:
print satState[(X<25.1)&(X>25)&(Y>.12)&(Y<.13)]
print cellSatH[(X<25.1)&(X>25)&(Y>.12)&(Y<.13)]
