On Oct 12, 2012, at 1:33 AM, Kris Kuhlman wrote:
> I am trying to simulate a non-linear problem with piecewise-defined
> parameters, depending on the range of the dependent variable.
>
> What is the best way to specify this? In my example below, I used a
> function and an intermediate CellVariable, using the "where" argument.
> I have also tried specifying the functions as one-liners using the
> numerix.where(h>hc,<true result>, <false result>) function, but this
> produces an error because the function does not return the right type
> of result.
You're working too hard. As you have found, setValue() assigns the value *right
then* to the CellVariable, with no further changes ever seen. Instead, just
write the expression you want explicitly:
param = <true result> * (h > hc) + <false result> * (h <= hc)
> I am also confused because putting a debugging print statement into
> the functions that define the parameters (e.g., C(h)) only prints out
> the message once when the problem is first set up. Are these functions
> not getting called each time to re-compute the parameters as the main
> variables are changing? Is FiPy stripping out or caching the print
> statements?
The functions are not called each time. When you write
DiffusionTerm(coeff=K(h)/C(h))
then K() and C() are called right then, each returning a constant CellVariable,
and the constant quotient is assigned to the coefficient of the DiffusionTerm.
FiPy has know way of knowing you want a dependency on h.
The script below illustrates the way I would do this.
Note that it crashes once h == hs because C is then zero and the DiffusionTerm
gets an infinite coefficient. Multiplying the equation through by C does not
work for some reason; Wheeler might know why.
--
from fipy import *
L = 10.0
nx = 75
dx = L/nx
timeStep = dx/10.0
alpha = 0.5
hs = -1.0/alpha
n = 0.2
expt = -(2 + 3*n)
Ks = 10.0
qs = 0.33
qr = 0.13
steps = 150
mesh = Grid1D(dx=dx, nx=nx)
h = CellVariable(name="$\psi$", mesh=mesh, value=-100.0, hasOld=True)
h.constrain(hs, where=mesh.facesLeft)
h.faceGrad.constrain((5*Ks,), where=mesh.facesRight)
h_lo = (h < hs)
h_hi = (h >= hs)
C = (-(qs-qr)*n*(-alpha*h)**(-n)/h) * h_lo
C.name = "C"
K = (Ks*(-alpha*h)**expt) * h_lo + Ks * h_hi
K.name = "K"
Se = (qr + (qs-qr)*(-alpha*h)**(-n)) * h_lo + qs * h_hi
Se.name = "$S_e$"
# Richards' equation (positive up/right)
Req = (TransientTerm(coeff=1.0) ==
DiffusionTerm(coeff=K/C) +
VanLeerConvectionTerm(coeff=[[1.0]]))
viewer = Viewer(vars=Se, datamin=qr, datamax=qs+0.02)
viewer.plot()
for step in range(steps):
h.updateOld()
Req.solve(dt=timeStep,var=h)
res = 100.0
while res > 1.0E-5:
res = Req.sweep(dt=timeStep,var=h)
viewer.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
