Thanks a lot for this example ! Your example is a very good representation of what I want to do if Ny > 1(where Ny is a number of species, we obviously have Ncell>1, see comments below). Vectorizing your example is not straightforward to me; I mean the case where Ny > 1, since the case where N\phi >1 (number of PDE) is well documented in your examples. Would you have a simple example to get me started? Bruno
I add some comments to help describing my problem further: In fact, I just didn't setup F to work witth the spatial dimension, i.e. F is for me a chemistry solver and does not know about space. \phi and y are like components and species (N\phi is the number of components and Ny is the number of species, and we have in general N\phi > 1 and always Ny > N\phi because of chemical reactions) Therefore the "local " F is in my case a vector of length Ny whereas the "globa" F is Ny*Ncell. Then , the "local " jacobian of F is full square matrix Ny*Ny whereas the "global" jacobian of F is a block diagonal matrix (Ny*Ncell)*(Ny*Ncell), whose blocks are indeed the Ny*Ny local Jacobian. Jonathan Guyer <[email protected]> Envoyé par : <[email protected]> 05/22/2012 07:50 PM Veuillez répondre à <[email protected]> A FIPY <[email protected]> cc Objet Re: Coupled PDE and Non-Liner Algebraic Equations On May 22, 2012, at 4:46 AM, <[email protected]> <[email protected]> wrote: > I'm trying to solve the following system of PDE and Non-Liner Algebraic Equations (NLAE): > > the PDEs read: > \frac{\partial{\phi}}{\partial t} = \nabla ( D(y) \nabla y) > where both phi and y are arrays of cellVariable and the number of PDE equals the length of the \phi array > > and the set of NLAE > > F(\phi, y)= 0 at every node in the domain > where y is the local root of F, i.e. \frac{\partial F}{\partial y} is a square full matrix. > This local problem is non-linear and I use SciPy.optimize.fsolve as solver. > > > My problem: > > I don't know an explict way of writing y = G{\phi} (and its derivatives) as suggested for example here: http://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.simple.html?highlight=expansion%20source The equations you have posed are fully explicit, so I'm not sure any of that is germane. > Do you have any suggestion for improving the performance of this resolution? and avoiding the use of setValue getValue for this kind of coupled problem? It's hard to say, as I'm not sure I accurately understand what you are trying to do. If I assume that D(y) = y^2 F(\phi, y) = \phi + y is this a fair representation of your code? import scipy.optimize from fipy import * mesh = Grid1D(nx=10, dx=1.) x, = mesh.getCellCenters() phi = CellVariable(mesh=mesh, name=r"$\phi$", value=x, hasOld=True) y = CellVariable(mesh=mesh, name="y", value=1., hasOld=True) D = y.getFaceValue()**2 eq = TransientTerm() == (D * y.getFaceGrad()).getDivergence() viewer = Viewer(vars=(phi, y)) def F(x, *args): return phi.getValue() + x for step in range(100): phi.updateOld() y.updateOld() res = 1e100 while res > 1e-3: res = eq.sweep(var=phi, dt=1e-4) y.setValue(scipy.optimize.fsolve(func=F, x0=y.getValue())) viewer.plot() _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] "This e-mail is confidential and may contain legally privileged information. If you are not the intended recipient, you should not copy, distribute, disclose or use the information it contains. Please e-mail the sender immediately and delete this message from your system. E-mails are susceptible to corruption, interception and unauthorised amendment; we do not accept liability for any such changes, or for their consequences. You should be aware, that the company may monitor your emails and their content."
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
