Hi Andrew,
Apologies for the delay in responding.
On Jun 26, 2007, at 7:10 PM, ASReeve wrote:
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
Generally, it is unusual for geometric quantities that are dependent
on the mesh to appear in the
undiscretized form of the equations. It makes me wonder if there is a
better way to write the equation.
How can the undiscretized form of the equations be dependent on
something associated with the mesh?
Certainly, when the divergence theorem is applied then that will
include the area of the face. You also
have the adjusted area of the face within the diffusion term.
'Rech' is a constant influx across the top of the model domain
(assigned to the top row of cells).
Again, that is a confusing way to think about the problem. Anything
being added to that equation should
be independent of the mesh. If a flux is required change Rech to be,
flux.getDivergence()
where "flux" flux is zero everywhere apart from the top surface.
'flux' should be initialized
as a VectorCellVariable in fipy. You can also use a FixedFlux
boundary condition.
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
Don't think that is a good way to set up the problem. I think we need
to define the equation as if
we were looking for an analytical solution. Also, why are we
including area in what looks like a 1D
solution?
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.
The equation above is 1D. Why do you need a 2D 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
Don't index directly into variables if at all possible. It shouldn't
be needed.
Besides, the ID ordering can change depending on the mesh. Rech
should be based on physical
position rather than grid cells.
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?
Let's sort out the problems I've highlighted above and hopefully it
will help
yield the correct answer.
Cheers
--
Daniel Wheeler
