Daniel,
Thanks for the reply. After playing with my code for the last couple of
days, everything seems to be working fine. I'm not sure what I changed to
correct the problem. Despite this, you raised a few points I'd like to have
clarified.
I'm attempting to simulate flow in a vertical profile, with recharge across
the top of the model. I'm treating recharge as a simple source term. I was
under the impression that the 'FixedFlux' boundary condition only worked
for no-flow boundaries. Am I mistaken?
Writing the 2-D equation out:
S dh/dt = d/dx { K_x (Z_sat / Z_cell) dh/dx} + d/dz(K_z dh/dz) +- source
z is vertical, x is horizontal
S is a storage constant (units (1/length)
K is a conductivity (Length/time)
h is hydraulic head (Length)
Z_cell is the height of the cell
Z_sat is the water level relative to the cell (base of cell has 0 relative
height), and is related to hydraulic head. I don't think there is an easy
way to linearize the equation.
Z_sat/Z_cell is a proportion of saturation, and as the water level in the
cell drops, the area that water can flow through drops and the flux across
the cell face drops. The method I'm using never allows this to go to zero,
instead, this ratio is made very small when the cell dewaters.
In my situation, the source is recharge and only occurs at the top of
the model (where rain enters the system). Why did you suggest using a
'getDivergence' term in my model? Is this similar in some way to
assigning a source term, or was my previous e-mail so confusing that
you misunderstood what I was trying to do?
I'm unsure how to specify, in a generic way, where the top row of cells
are without directly indexing the cells. I understand I can use the
'where' functionality to define parameters in a cell variable, but don't
know how to specify the coordinates for the bottom and top of a particular
row of cells. I'm also using indexing to retrieve the cell
elvations using:
topBot=mesh.getFaceCenters()[...,1]
and then use topBot[:nx*nz] to define the 'bottom elevation of each cell
and topBot[nx:nx*(nz+1)] to define the top elevation of each cell. If this
indexing isn't the best way to do this, what should I be doing? For a
structured mesh, I could just assign recharge using something like
rech.setValue(1e-6/delZ, where=(mesh.getCellCenters[...,1]>=lenZ-delZ/2.))
but this won't work in unstructured meshes.
Thanks,
Andy
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
(207)581-2353
[EMAIL PROTECTED]
