I'm attempting to solve a nonlinear problem using sweeps. I'm solving for
hydraulic head, and the flux is influenced by the head gradient and the
cross sectional area of a cell, which is influenced by hydraulic head:
S dh/dt = d/dx ( K*A(h)*dh/dx) +Rech
'Rech' is a constant influx across the top of the model domain (assigned to the
top row of cells).
As I can't change the area of the cell face directly, I'm making a
cell variable ('SatX') that represents the proportion cell area needs
to be decreased by:
S dh/dt = d/dx( K*SatX(h)*A dh/dx) + Rech
I can get the my code to run but I'm getting the wrong answer. I
import the appropriate modules from fipy and setting up a 2D uniform
mesh. I assign no flow boundaries to the top and bottom and type 1
boundaries to the left and right edges of the model domain. A gigantic
time step is assigned to simulate steady-state conditions and I take
one time step. I create several variables:
head=CellVariable(name='hydraulicHead',mesh=mesh,value=10.,hasOld=1)
spStore=CellVariable(name='specific storage',mesh=mesh,value=1.e-4)
rech=CellVariable(name='horizontal saturation',mesh=mesh,value=0.)
rech[-nx:,...]=1.e-6/delZ ##set rech for top row
kond=VectorCellVariable(name='hydraulicConductivity',mesh=mesh)
kondX=1e-3 ##base cond values for model
kondY=1e+1
and define a function to calculate 'SatX':
minThick=1e-10 ## based on Doherty, 2001
residThick=1.e-2
f=6.0
g=(1./residThick)-f
topBot=mesh.getFaceCenters()[...,1]
def SatX(head,topBot,residThick,minThick,f,g):##get top and bottom elevation of
cells
satX=CellVariable(name='horizontal saturation',mesh=mesh,value=1.,hasOld=1)
satThick=CellVariable(name='saturated thickness',mesh=mesh,hasOld=1)
satThick=N.minimum((head[:]-topBot[:nx*nz])/(topBot[nx:nx*(nz+1)]-topBot[:nx*nz]),1.)
satX.setValue(residThick*N.exp(-g*satThick[:])+satThick[:],
where=((head[:]<topBot[nx:nx*(nz+1)])
&(head[:]>topBot[:nx*nz])))
satX.setValue(N.maximum(residThick*N.exp(f*satThick[:]),minThick),
where=((head[:]<=topBot[:nx*nz])))
return satX
Solving my equation:
eqn=TransientTerm(coeff=spStore)==ImplicitDiffusionTerm(coeff=kond)+rech
for tstep in xrange(nSteps):
res=1e+5
head.updateOld()
while res>1.e-10:
kond[...,0]=kondX*SatX(head,topBot,residThick,minThick,f,g)
kond[...,1]=kondY
res=eqn.sweep(var=head,boundaryConditions=boundCond,dt=delT)
The 'head' at the center of the model should be 27.8, and FiPy yields 35.4.
Does 'SatX' need to be included in the FiPy equation? What am I doing
wrong?
Andy
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
(207)581-2353
[EMAIL PROTECTED]
