Andrew,
Rather than trying to debug what's happening in your script line-by-
line, I thought it would be more productive for me to show how I would
extract the quantities you want from your equations. Unfortunately, I
got stumped pretty quickly even doing that.
My first step was to figure out exactly what `rech` represents, where
I found
topFaceAreas=mesh._getFaceAreas()[mesh._getCellFaceIDs().data[5,:]]
This is not at all how we intended FiPy to work. I won't say it never
happens, but the user should never have to call any routine that
starts with `_` [*] and certainly shouldn't try to slice particular
indices out of any FiPy field.
If rech is doing what I think it's doing, I think it should be written
as the divergence of a flux, where the flux is only non-zero at the
inlet surfaces. Basically:
X, Y, Z = mesh.getFaceCenters()
rech = ((Z == Z.max()) * [[0.], [someFlux]]).getDivergence()
where `someFlux` might be your 1.5e-8, but I don't understand what
"need to change units from L/T to (Vol H2O)/(Vol Cell)/time" means.
You should not generally need to account for face areas and cell
volumes; FiPy's operators usually do that for you.
If this isn't right (or even if it is), what's the mathematical
expression for this recharge flux?
As an aside, I've been trying to establish a convention in my own
codes and in the examples of using
x, y, z = mesh.getCellCenters()
X, Y, Z = mesh.getFaceCenters()
Both sets of coordinates are then available and it's immediately clear
whether you're talking about cell quantities or face quantities.
[*] if you find you need access to an internal `_` method, please file
a trouble ticket or ask on the list so that we can either direct you
to an alternative or we can shift the call into the public API.
