On Tue, Aug 19, 2008 at 07:44:00PM +0200, Jed Brown wrote: > On Tue 2008-08-19 18:32, Anders Logg wrote: > > On Tue, Aug 19, 2008 at 03:40:13PM +0200, Jed Brown wrote: > > > On Tue 2008-08-19 14:06, Anders Logg wrote: > > > > On Tue, Aug 19, 2008 at 01:49:22PM +0200, Jed Brown wrote: > > > > > On Tue 2008-08-19 13:40, Anders Logg wrote: > > > > > > On Tue, Aug 19, 2008 at 12:12:50PM +0200, Jed Brown wrote: > > > > > > > On Tue 2008-08-19 11:59, Anders Logg wrote: > > > > > > > > On Thu, Aug 14, 2008 at 10:10:03PM +0000, Jed Brown wrote: > > > > > > > > > One way to implement this is to allocate a vector for > > > > > > > > > Dirichlet values, > > > > > > > > > a vector for Homogeneous values, and a Combined vector. The > > > > > > > > > Homogeneous > > > > > > > > > vector is the only one that is externally visible. > > > > > > > > > > > > > > > > Isn't this problematic? I want the entire vector visible > > > > > > > > externally > > > > > > > > (and not the homogeneous part). It would make it difficult to > > > > > > > > plot > > > > > > > > solutions, saving to file etc. > > > > > > > > > > > > > > > > Maybe the Function class could handle the wrapping but it would > > > > > > > > involve a > > > > > > > > complication. > > > > > > > > > > > > > > Right, by `externally visible' I mean to the solution process, > > > > > > > that is > > > > > > > time-stepping, nonlinear solver, linear solvers, preconditioners. > > > > > > > The > > > > > > > vector you are concerned about is the post-processed state which > > > > > > > you can > > > > > > > get with zero communication. It is inherently tied to the mesh > > > > > > > and > > > > > > > anything you do with it likely needs to know mesh connectivity. > > > > > > > I don't > > > > > > > think it is advantageous to lump this in with the global state > > > > > > > vector. > > > > > > > > > > > > > > Jed > > > > > > > > > > > > I don't understand. What is the global state vector? > > > > > > > > > > The global state vector is the vector that the solution process sees. > > > > > Every entry in this vector is a real degree of freedom (Dirichlet > > > > > conditions have been removed). This is the vector used for computing > > > > > norms, applying matrices, etc. When writing a state to a file, this > > > > > global vector is scattered to a local vector and boundary conditions > > > > > are > > > > > also scattered into the local vector. The local vector is serialized > > > > > according to ownership of the mesh (you have to do this anyway). > > > > > > > > > > Jed > > > > > > > > I'm only worried about how to create a simple interface. Now, one may > > > > do > > > > > > > > u = Function(...); > > > > A = assemble(a, mesh) > > > > b = assemble(L, mesh) > > > > bc.apply(A, b) > > > > solve(A, u.x(), b) > > > > plot(u) > > > > > > > > How would this look if we were to separate out Dirichlet dofs? > > > > > > How about a FunctionSpace object which manages this distinction. > > > Something like the following should work. > > > > > > V = FunctionSpace(mesh, bcs); // Is this name clearer? > > > > We've discussed introducing a FunctionSpace concept earlier (on ufl-dev) > > to handle boundary conditions, and to enable sharing of function space > > data like meshes and dof maps. This might be a good idea, but it has > > to be something like > > > > V = FunctionSpace(element, mesh, bcs) > > Yes, of course. > > > > u = V.function(); > > > > I think this should be > > > > u = Function(V) > > Sure, but the representation is fairly closely tied to information in V > which might not otherwise need to be public. > > > > A = V.matrix(a); > > > P = V.matrix(p); // preconditioning matrix, optional [1] > > > > What do these accomplish? Return a matrix of appropriate size? I don't > > think that's necessary since the assembler can set the size. > > It's no problem to have the assembler check if the matrix is allocated > and create it if necessary. My example was just being explicit that > such allocation only needs to happen once. > > > > solver = LinearSolver(A, P); > > > u0 = u.copy(); // if nonlinear, set initial guess > > > > > > V.assemble(A, a, u0); // builds Jacobian matrix [2] > > > V.assemble(P, p, u0); // preconditioning form, optional > > > V.assemble(b, L, u0); // the ``linear form'' (aka nonlinear residual) > > > solver.solve(u, b); > > > > The problem here is that the solver needs to be aware of Functions. > > Now, a solver only knows about linear algebra, which is nice. > > Whatever syntax you like. The solver is just going to grab the global > representation and work with that, the line could just as well read > > solver.solve(u.x(), b.x()); > > My suggestion (relating to your other message about regarding > simplicity) is that solver components always call (in the language > presented below) .dofs() on their arguments while output components > always call .vector(). Then the user only sees the high lever Function > object. > > > > Note that u0 disappears if the problem is linear and P disappears if you > > > use the real Jacobian as the preconditioning matrix. The key point is > > > that boundary conditions are built into the function space. It is no > > > more work, it just happens in a different place. > > > > > > [1] It is essential that the test/trial spaces for the preconditioning > > > matrix are the same as for the Jacobian. > > > > > > [2] I don't like the implicit wiring of the state vector into > > > assemble(). It is really hard to track dependencies and the current > > > state is not a property of the Bilinear/Linear forms. > > > > Maybe one could do > > > > (A, b) = assemble(a, L, mesh, bc) > > solve(A, u.dofs(), b) > > > > This would be in addition to the current > > > > (A, b) = assemble(a, L, mesh) > > bc.apply(A, b) > > solve(A. u.vector(), b) > > > > In the first version, the assembler knows about the boundary > > conditions and in the second it doesn't. This would require having two > > members in the Function class that return a Vector: > > > > u.vector() // Returns entire vector > > u.dofs() // Returns non-Dirichlet dofs > > The Function object can keep track of which version is current and do > the necessary scatter to give the version you ask for. > > We seem to have a slightly different philosophy of `simple'. I would > rather write a couple extra lines which reveal the control flow and > enable doing more complicated things. For instance, I'd rather avoid > magic such as updating a special vector, put into the Bilinear form at > some earlier point, to base the new Jacobian. Similarly, if you solve > more than one problem then the (A,b) = assemble(a,L,mesh) version cannot > reuse the matrix. You can always support both versions but I think this > level of overloading makes it harder, rather than simpler, for the user > to keep track of what is happening. When someone says `how do I do ...' > and the answer is to basically refactor their code rather than just > insert the required line, it looks less like simplicity and more like > magic. > > Jed
Yes, I have a tendency towards the "magic" end of the spectrum. I admit this may cause trouble sometimes. On the other hand, I don't want things to look like PETSc where one needs to write large chunks of code for simple things. -- Anders
signature.asc
Description: Digital signature
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
