On May 23, 2012, at 1:21 PM, <[email protected]>
<[email protected]> wrote:
> > Please give the simplest set of toy equations for your PDEs and for F that
> > represent what you want.
> Here's an example with N\phi = 2 and Ny = 3
> PDEs
> \frac{\partial{\phi_1}}{\partial t} = \nabla ( D(y_1) \nabla y_1) + \nabla (
> D(y_3) \nabla y_3)
> \frac{\partial{\phi_2}}{\partial t} = \nabla ( D(y_2) \nabla y_2) + \nabla (
> D(y_3) \nabla y_3)
> NLAEs (at every point in space)
> \phi_1 - exp(y1) - exp(y3) = 0
> \phi_2 - exp(y2) - exp(y3) = 0
> -a3(y1,y2,y3) - y3 + a2(y1,y2,y3) + y2 + a1(y1,y2,y3) + y1 + constant = 0
>
> Again thanks a lot !!!!!!!
OK, that's clearer to me, now, thanks.
For small values of N\phi and Ny, I think the following would do the trick. For
larger numbers of variables and equations, I have some ideas, but would need to
experiment a little. Either way, I recommend you read
http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#root-finding-large-problems
----
import scipy.optimize
from fipy import *
mesh = Grid1D(nx=10, dx=1.)
N = mesh.getNumberOfCells()
phi = [CellVariable(mesh=mesh, name=r"$\phi_1$", hasOld=True),
CellVariable(mesh=mesh, name=r"$\phi_2$", hasOld=True)]
y = [CellVariable(mesh=mesh, name="y_1", value=1., hasOld=True),
CellVariable(mesh=mesh, name="y_2", value=1., hasOld=True),
CellVariable(mesh=mesh, name="y_3", value=1., hasOld=True)]
eq = [TransientTerm() == (D(y[0]) * y[0].getFaceGrad()).getDivergence() +
(D(y[2]) * y[2].getFaceGrad()).getDivergence(),
TransientTerm() == (D(y[1]) * y[1].getFaceGrad()).getDivergence() +
(D(y[2]) * y[2].getFaceGrad()).getDivergence()]
viewer = Viewer(vars=phi + y)
def F(x, *args):
y = [x[0:N], x[N:2*N], x[2*N:3*N]]
phi, = args
return numerix.concatenate([phi[0].getValue() - numerix.exp(y[0]) -
numerix.exp(y[2]),
phi[1].getValue() - numerix.exp(y[1]) -
numerix.exp(y[2]),
-a3(y) - y[2] + a2(y) + y[1] + a1(y) + y[0] +
K])
for step in range(100):
for v in phi + y:
v.updateOld()
res = 1e100
while res > 1e-3:
res1 = eq[0].sweep(var=phi[0], dt=1e-3)
res2 = eq[1].sweep(var=phi[1], dt=1e-3)
res = max(res1, res2)
yval = scipy.optimize.fsolve(func=F,
x0=numerix.concatenate([yN.getValue() for
yN in y]),
args=(phi,))
y[0].setValue(yval[0:N])
y[1].setValue(yval[N:2*N])
y[2].setValue(yval[2*N:3*N])
viewer.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
