On Mon, Jun 27, 2011 at 5:11 AM, Michael McNeil Forbes <[email protected]> wrote: > > Dear FiPyers, > > A few problems/confusions I am having: > > 1) I am a bit confused about using ImplicitSourceTerm. I thought the idea > was that you could linearize a non-linear term and have control over what > parts were placed in the matrix.
Essentially, this is true, but it automates the reversion to explicit in the event that the source term is negative. We've debated whether this is the best approach and settled on what we have. In our work we found that we were doing this calculation all the time so we made it part of the source. Basically if the sign is negative on the LHS (with a positive equation) it gets added explicitly and if the sign is positive on the LHS (with a positive equation) it gets added implicitly. See <http://matforge.org/fipy/browser/trunk/fipy/terms/implicitSourceTerm.py#L71> Currently, there is no way to switch this behavior off short of editing implicitSourceTerm.py. Maybe we should add that as an option or have an UnadulteratedImplicitSourceTerm of some kind. To turn off the automation, replace (*) line 85 to 88 in <http://matforge.org/fipy/browser/trunk/fipy/terms/implicitSourceTerm.py#L85> with return {'diagonal' : numerix.ones(var.shape, 'd'), 'old value' : numerix.zeros(var.shape, 'd'), 'b vector' : numerix.zeros(var.shape, 'd'), 'new value' : numerix.zeros(var.shape, 'd')} The above is for trunk, which you are probably not using, but the idea will be the same in other revisions. (*) Actually, you can just subclass from ImplicitSourceTerm and overwrite _getWeight. > Thus, I would have thought that if I set: > > K = DiffusionTerm() > A = CellVariable(...) > B = CellVariable(...) > .. > X = CellVariable(..., hasOld=1) > > eq = K + A*X == B > > then eq.solve(X) solve etc. would do something like > > X -> K^{-1}(B - A*X) Seems correct. > but > > eq = K + ImplicitSourceTerm(A) == B > > would do > > X -> (K + A)^{-1}B Depends on the sign of A. Try making A all positive and all negative just to test. Make the modification to implicitSourceTermp.py and it should work as you want. > 2) The ability to couple equations with & (or possibly +), and by specifying > variables looks really nice, but does not seem to appear anywhere in the > documentation. Is there anything I should be aware about trying to use > this? We haven't released this yet and we still have to add the documentation. I've also just added the ability to do vector calculations. See <http://matforge.org/fipy/browser/trunk/examples/riemann/acoustics.py> as a prelude to adding Riemann flux updates in FiPy. > I.e. > > eq_x = DiffusionTerm(var=X) + ImplicitSourceTerm(B,var=Y) == D > eq_y = DiffusionTerm(var=Y) + ImplicitSourceTerm(C,var=Y) == E > eqs = eq_x & eq_y > eqs.solve() > > Again, I was hoping to have enough control to ensure that this would try to > solve > > [X] = [K, B]^{-1} [D] > [Y] [C, K] [E] This should work fine for you with the modifications if you subclass ImplicitSourceTerm. > 3) If FiPy itself cannot do this, is there an easy way of extracting the > matrices so I can solve them myself? In particular, I don't really want to > have to run through a whole solver iteration to generate the matrix in the > usual way. Is there some way of calling something like eqs.buildMatrix() to > just get the matrix so I can solve it outside of FiPy? Related question: > once I have the matrix, is there a simpler way of getting the full > non-sparse array? See if the above works for you and then we can go from there. > Basically, what is the most efficient way of doing something like: > > eqs.cacheMatrix() > eqs.solve() # but don't actually do any solving! > M = pysparse.sparse.PysparseMatrix(matrix=eqs.matrix.matrix).getNumpyArray() > > where > > M = [K, B] > [C, K] > > (As it stands, I think M = diag(K, K)*dx with B and C being combined with D > and E...) Using justResidualVector() does the same as solve, but doesn't solve the matrix. from fipy import * m = Grid1D() v = CellVariable(mesh=m) e = DiffusionTerm() e.cacheMatrix() e.justResidualVector(v) print e.matrix > I would like to formulate the problem in fipy (because later I will use it > for time evolution) but need to solve for the static configurations first. > (Also, is there a way of not using sparse solvers and just using numpy > arrays? For debugging small systems this would be very useful...) You can ask for "e.matrix.numpyArray" to get the matrix as a numpy array if you prefer that format. Hope this helps. Cheers. -- Daniel Wheeler
