FiPy list,
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.
I have not seen any examples with this type of behavior. Is there a
best way to do this? The way I am doing it seems sub-optimal.
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?
I appreciate any help or pointers to applicable documentation and examples.
Kris
--
from fipy import *
import numpy as np
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)
Se = CellVariable(name="$S_e$", mesh=mesh, value=0.0)
h.constrain(hs, where=mesh.facesLeft)
h.faceGrad.constrain((5*Ks,), where=mesh.facesRight)
def C(h):
r = CellVariable(mesh=mesh)
result = -(qs-qr)*n*(-alpha*h)**(-n)/h
r.setValue(result)
r.setValue(0.0, where=h>=hs)
return r
def K(h):
r = CellVariable(mesh=mesh)
result = Ks*(-alpha*h)**expt
r.setValue(result)
r.setValue(Ks, where=h>=hs)
return r
def fSe(h):
r = CellVariable(mesh=mesh)
result = qr + (qs-qr)*(-alpha*h)**(-n)
r.setValue(result)
r.setValue(qs, where=h>=hs)
return r
# Richards' equation (positive up/right)
Req = (TransientTerm(coeff=1.0) ==
DiffusionTerm(coeff=K(h)/C(h)) +
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)
Se.setValue(fSe(h))
viewer.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
