Jonathan,
Thanks for the quick response. The '_FixedBCFaceGradVariable' method was
what I was looking for, but I still can't seem to get a reasonable mass
(water) balance for my simulation.
What I'd like to get is a list of all fluxes through sources and sinks to
the model, a list of all the changes in water stored within each model
cell, and then a comparision between these to verify that the amount of
water lost (or gained) through sources and sinks is equal to the change in
water stored within the model. I initially was doing this just to look at
where water was coming from within the model but noticed large mass
balance errors.
I've tried to do this several ways with no success in closing the
difference between the sum of sources/sinks and the sum of the changes in
stored water. I must be doing something wrong when I set up these
calculations.
I've attached my script for doing this. If you run it, you'll see the
following output:
1) change in water stored (volume of water released from storafe over
the time step)
2) rech (source term, volume of water entering the top row of cells
across the model.)
3) faceFluxes (fluxes across type 1 boundaries, calculated from
divergence of cells abutting fixedValue faces;
flux=divergence-storageChange)
4) sum of fluxes (sum of fluxes plus sum of storage changes, should
equal zero)
5) BC discharge (sum of fluxes across faces using
'_FixedBCFaceGradVariable' method.)
I consistently get large mass balance errors between the sum of sources
and change in storage.
Andy
Andrew Reeve
Associate Professor
Dept. of Earth Sciences
University of Maine
207-581-2353
On Tue, 12 May 2009, Jonathan Guyer wrote:
On May 12, 2009, at 11:00 AM, A.S.Reeve wrote:
Is there a simple way to calculate the change in mass stored in a cell
over a single time step?
Unless I'm missing something, I think this is just (mass - mass.getOld()),
isn't it?
I'm trying to do this using a flux assigned
to a FaceVariable with the operations (kond and head are CellVaribles):
flux=-head.getFaceGrad()*kond.getHarmonicFaceValue()
flux.getDivergence()
This seems to work at internal cells but I'm having trouble sorting this
out at cells with fixedValue boundary condtions assigned to one of their
faces. How does one seperate changes in storage in these cells from influx
across face with a fixedValue conditions? Does the divergence operator
in FiPy lump these two terms (storage and flux across face) together in
these cells?
Neither the gradient nor divergence operators see boundary conditions by
default. This is on our to-do list. There is a workaround for the gradient,
though, which we originally wrote up for Zhiwen Liang, and which has since
been included in FiPy 2.0. You would write:
> > > from fipy.variables.fixedBCFaceGradVariable import
> > > _FixedBCFaceGradVariable
> > > flux=-_FixedBCFaceGradVariable(head,
> > > boundaryConditions=???)*kond.getHarmonicFaceValue()
where the boundaryConditions argument takes the same tuple that you pass to
solve. The divergence should be fine.
That being said, I don't think you need any of this for this particular task.
from fipy import *
from fipy.variables.fixedBCFaceGradVariable import _FixedBCFaceGradVariable
import cPickle
def WaterBalance(head,delT,waterVol,fixedFlux,fixedValueCells,divFlux):
###calc volumes of water exchanged over a time step
storeChng=(waterVol-waterVol.getOld())*mesh.getCellVolumes()
sourceFlux=[delT*sum(param()*mesh.getCellVolumes()) for param in fixedFlux]
sourceName=[param.name for param in fixedFlux]
faceFlux= ([delT*divFlux()[fvc]-storeChng[fvc]
for fvc in fixedValueCells])
print '#####'
print 'Change in Water Stored=',sum(storeChng)
for i in range(len(sourceFlux)):
print sourceName[i],sourceFlux[i]
for i in range(len(faceFlux)):
print 'faceFlux ',i, sum(faceFlux[i])
print 'sum of fluxes', sum(sum(faceFlux))+sum(sourceFlux)-sum(storeChng)
print '#####'
return storeChng
###set paramters to change
kAcroTop=5.e-3; kAcroBot=2.e-5; kCato=1.e-6 #hydr. conductivities
kVkH=0.1 #anisotropy
###Desat Params
minDesat=-4.
residThick=1.e-2
f=6.0
g=1./residThick-f
maxSatV=100.
###grid params
nx=10;ny=1;nz=7
dx=[10. for i in range(nx)]
dy=[1. for i in range(ny)]
dz=[1. for i in range(nz)]
mesh=Grid3D(dx=dx,dy=dy,dz=dz,nx=nx,ny=ny,nz=nz)
X,Y,Z=mesh.getFaceCenters()
### ID type 1 boundaries
fixedValues=[5.1+sum(dx)*.001,5.1]
fixedValueFaces=[(mesh.getExteriorFaces()&(X==X.min())&(Z<fixedValues[0])),
(mesh.getExteriorFaces()&(X==X.max())&(Z<fixedValues[1]))]
### ID cells ajacent to fixed value faces
fixedValueCells=[]
for fvf in fixedValueFaces:
fixedValueCells.append(extract(fvf,mesh.getFaceCellIDs()[0]))
GwBCs = (FixedValue(fixedValueFaces[0],fixedValues[0]),
FixedValue(fixedValueFaces[1],fixedValues[1]))
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]])
thick=CellVariable(name='thick',mesh=mesh, value=top-bot)
topFaceAreas=mesh._getFaceAreas()[mesh._getCellFaceIDs().data[5,:]]
X,Y,Z=mesh.getCellCenters()
head=CellVariable(name='head',mesh=mesh,value=fixedValues[0]+.5,hasOld=1)
### Storage paramters
spStore=CellVariable(name='Spec.Store',mesh=mesh,value=5.e-3)
spYield=CellVariable(name='Spec.Yield',mesh=mesh,value=0.1)
desat=CellVariable(name='Spec.Store',mesh=mesh,value=1.e-20)
### Desaturation Variables
satState=(((head-bot)/(thick)<minDesat)*minDesat+
(((head-bot)/(thick)>=minDesat)&(head<=top))*(head-bot)/(thick)+
(head>top)*1.)
### storage due to confined and unconfined conditions, with added 'desaturated
storage'
Store=(desat+(satState>=0.)*(spStore+(satState<1.)*spYield))
cellSatH=CellVariable(name='Horiz. K. Mult',mesh=mesh,hasOld=0)
cellSatH.setValue((satState>=0)*(satState+residThick*exp(-g*satState))+
(satState<0)*residThick*exp(f*satState))
cellSatV=((satState<-1)*maxSatV+
((satState>-1)&(satState<0))*1/(satState+1.01)+
1.)
### Calc total amount of water to estimate change in stored water
waterVol=((head-bot)*(satState>0.)*(spStore+(satState<1.)*spYield))
### Enter hydraulic conductivity values for model
kond=CellVariable(name='hydr.cond.',mesh=mesh,rank=1,value=kAcroTop)
kond.setValue([where(Z<4.5,kAcroBot,kond[0]),
where(Z<4.5,kAcroBot,kond[1]),
where(Z<4.5,kVkH*kAcroBot,kond[2])])
kond.setValue([where(Z<4.,kCato,kond[0]),
where(Z<4.,kCato,kond[1]),
where(Z<4.,kVkH*kCato,kond[2])])
### recharge of water across top of model. need to change units from L/T to
(Vol H2O)/(Vol Cell)/time
rech=CellVariable(name='rech', mesh=mesh,value=
(Z==Z.max())*1.5e-8*topFaceAreas/mesh.getCellVolumes())
### calculate flux (specific discharge) and total discharge
flux=-head.getFaceGrad()*(kond*(cellSatH*[1,1,0]+cellSatV*[0,0,1])).getHarmonicFaceValue()
discharge=flux*mesh._getFaceAreas()###volume of H2O/time
### make list of all sources for water balance calculation
fixedFlux=[rech]
GwEqInj=(TransientTerm(coeff=Store)
==DiffusionTerm(coeff=(kond*(cellSatH*[1,1,0]+
cellSatV*[0,0,1])).getHarmonicFaceValue())
+rech)
delT=1.e2
for i in range(2):
head.updateOld()
resid=1.
while resid>1.e-14:
resid=GwEqInj.sweep(var=head,
solver= LinearGMRESSolver(tolerance=1.0e-16,
iterations=500),
dt=delT,
underRelaxation=None,
boundaryConditions=GwBCs
)
print resid
storeChng=WaterBalance(head,delT,waterVol,fixedFlux,fixedValueCells,discharge.getDivergence())
###calc boundary face fluxes using new method
dischargeBC=(mesh._getFaceAreas()*
_FixedBCFaceGradVariable(head,boundaryConditions=GwBCs)*
(kond*(cellSatH*[1,1,0]+cellSatV*[0,0,1])).getHarmonicFaceValue())
print 'BC discharge=',delT*sum(sum(dischargeBC()))
cPickle.dump(head,open('head.pkl','w'))
cPickle.dump(kond(),open('kond.pkl','w'))
