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 fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]