Daniel, Jonathan and other fipy users:
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.
However, when I calculate the cellSat term as:
kond*cellSatH+cellSatV*kond*[0,0,1]
my computer model zips along quite nicely and seems to work correctly,
for the few situations I've tried (cellSatV is either close to
one or much larger than cellSatH, at least for most situations). This
seems backwards to me, and I'm wondering why the first approach
does not work?
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.
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've attached the code I'm using. The print statements in the sweeps
prints out values for one vertical column of cells.
Andy
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
from fipy import *
from time import time
startTime=time()
nx=41
ny=3
nz=14
dx=[0.0675*1.3**abs(i-17) for i in range(nx)]
dy=[0.0675*1.3**abs(i-(ny-1)/2.) for i in range(ny)]
dz=[0.5 for i in range(nz)]
mesh=Grid3D(dx=dx,dy=dy,dz=dz,nx=nx,ny=ny,nz=nz)
X,Y,Z=mesh.getFaceCenters()
BCs = (FixedValue(mesh.getExteriorFaces()&(X==0)&(Z<5.6),5.6),
FixedValue(mesh.getExteriorFaces()&(X==max(X))&(Z<2.9),2.9))
top=CellVariable(name='top of cell',mesh=mesh,
value=Z[mesh._getCellFaceIDs().data[5]])
bot=CellVariable(name='bottom of cell',mesh=mesh,
value=Z[mesh._getCellFaceIDs().data[4]])
X,Y,Z=mesh.getCellCenters()
head=CellVariable(name='head',mesh=mesh,value=8.,hasOld=1)
'''
These terms control teh desaturation of a cell, slowly decreasing
the permeability in the x and y directions as the cell dewaters.
permeability in the z direction is increased as the cell dewaters
to allow sorce/sink terms in the top layer to continue to influence
the saturated zone.
'''
###need to add this to satState variable
#minDesat=-3. ## need to specify minimum satState or will get overflow
residThick=1.e-2
f=6.0
g=1./residThick-f
satState=(head<top)*(head-bot)/(top-bot)+(head>top)
cellSatH=((satState>=0)*(satState+residThick*exp(-g*satState))+
(satState<0)* residThick*exp(f*satState))
cellSatV=(satState>=.01)/(satState)+ (satState<.01)*100
spStore=1.e-3
spYield=.2
desat=1.e-8
storage=desat+(head>=top)*spStore+((head>bot)&(head<top))*spYield
kond=CellVariable(name='hyd.cond.',mesh=mesh,rank=1,value=1.e-5)
print 'a',time()-startTime
###second method
#cellSat=kond*[cellSatH,cellSatH,cellSatV]
print 'b',time()-startTime
###recharge across top of model
rech=CellVariable(name='recharge',mesh=mesh,
value=where(Z==Z.max(),1.e-8,0))
print 'c',time()-startTime
###second method
#GwEqInj=TransientTerm(coeff=storage)==DiffusionTerm(coeff=cellSat)+rech
GwEqInj=TransientTerm(coeff=storage)==DiffusionTerm(coeff=kond*cellSatH+cellSatV*kond*[0,0,1])+rech
print 'd',time()-startTime
delT=1.e7
for i in range(2):
resid=1.
print '****************'
while resid>1.e-3:
head.updateOld()
print satState()[where((X<25.1)&(X>25)&(Y>.12)&(Y<.13))]
print cellSatH()[where((X<25.1)&(X>25)&(Y>.12)&(Y<.13))]
resid=GwEqInj.sweep(var=head,
solver= LinearPCGSolver(tolerance=1.0e-10,
iterations=1000),
dt=delT,
underRelaxation=None,
boundaryConditions=BCs
)
print resid
