Daniel,
Looks like I did what you described. I used _getFaceAreas() instead of
getFaceAreas() to calculate the face areas. I aslo used both the method
Jonathan described to calculate fluxes across Type 1 boundarues,
calculating face gradients using
'_FixedBCFaceGradVariable(head,boundaryConditions=GwBCs)'
And also extimating this flux by calulating the divergence in cells
associated with these boundary faces and subtracting the storage change.
FluxThroughBndryFace=Div(Discharge)-ChngInStorage
This latter method results in a better balance than the method Jonathan
described. I intially thought this was a problem with when I calculate
fluxes, caluclating them at the end of the time step instead of averaging
the fluxes at the start and of the time step, but this doesn'tseem to be
the case.
I'm likely just doing something stupid...
Andrew Reeve
Associate Professor
Dept. of Earth Sciences
University of Maine
207-581-2353
On Tue, 19 May 2009, Daniel Wheeler wrote:
Andrew,
I don't think I have time right now to dig through your code. I'll
keep it in my inbox and hopefully I'll get to it at some point. One
thing I would say is that I think the mass balance should be checked
via
var * mesh.getCellVolumes() - var.getOld() * mesh.getCellVolumes()
the value returned by that calculation should be equal to
allSourceTerms * mesh.getCellVolumes() * dt + allFluxesIn *
mesh.getFaceAreas() * dt
where you need to construct allSourceTerms and allFluxesIn as a
CellVariable and FaceVariable, respectively. You may have in fact done
this in another way. Is this basically what you did?
On Mon, May 18, 2009 at 11:56 AM, A.S.Reeve <[email protected]> wrote:
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.
--
Daniel Wheeler
