2008/5/19 Jed Brown <[EMAIL PROTECTED]>: > On Mon 2008-05-19 08:41, Johan Hoffman wrote: >> > 2008/5/19 Johan Hoffman <[EMAIL PROTECTED]>: >> >>> On Sun 2008-05-18 22:55, Johan Hoffman wrote: >> >>>> > On Sun 2008-05-18 21:54, Johan Hoffman wrote: >> >>>> >> > On Sat, May 17, 2008 at 04:40:48PM +0200, Johan Hoffman wrote: >> >>>> >> > >> >>>> >> > 1. Solve time may dominate assemble anyway so that's where we >> >>>> should >> >>>> >> > optimize. >> >>>> >> >> >>>> >> Yes, there may be such cases, in particular for simple forms >> >>>> (Laplace >> >>>> >> equation etc.). For more complex forms with more terms and >> >>>> coefficients, >> >>>> >> assembly typically dominates, from what I have seen. This is the >> >>>> case >> >>>> >> for >> >>>> >> the flow problems of Murtazo for example. >> >>>> > >> >>>> > This probably depends if you use are using a projection method. If >> >>>> you >> >>>> > are >> >>>> > solving the saddle point problem, you can forget about assembly >> >>>> time. >> >>>> >> >>>> Well, this is not what we see. I agree that this is what you would >> >>>> like, >> >>>> but this is not the case now. That is why we are now focusing on the >> >>>> assembly bottleneck. >> >>> >> >>> This just occurred to me. If you have a Newtonian fluid, then the >> >>> momentum >> >>> equations are block diagonal, but this is not reflected in the matrix >> >>> structure. >> >>> Sure enough, run the stokes demo with -mat_view_draw -draw_pause -1 and >> >>> note >> >>> that the off-diagonal blocks of the momentum equations are cyan which >> >>> means they >> >>> are set, but have value zero. This almost doubles the number of >> >>> insertions into >> >>> the global matrix. >> >> >> >> Good. You are right, this piece of information is not used. >> >> >> >> I guess the most general thing is to have ffc delete zero matrix entries >> >> in computing the sparsity pattern. I do not think this is done today? >> > >> > We could add it with an appropriate >> > ufc::dof_map::tabulate_sparsity_pattern(...), >> > since the form compiler can figure out which entries are always zero. >> > Currently, this information is hidden from dolfin, and therefore it must >> > simply use the local blocks. >> >> Yes, this would be very useful to add. >> >> > But A.add(...) calls a single block-addition function in PETSc >> > or Trilinos, does anyone know how these will perform if the values >> > contain zeros that are outside the initialized sparsity pattern? >> >> Worst case scenario is that a reallocation is triggered, but maybe it is >> dealt with in a less drastic way? > > It is important to preallocate for all the entries you will be inserting. If > the option MAT_IGNORE_ZERO_ENTRIES is set then zero entries will not be > inserted. I don't think you want to rely on this as a general mechanism for > eliminating non-existant coupling since it would also eliminate entries that > just happened to be missing for the current Jacobian. > > Note that MatSetValues() needs to be called with a ``full'' block so if we > eliminate the coupling zeros, we must make separate calls for the rows of each > vector component. However, this should not make much performance difference
Ok, this was what I wanted to know. And if the ongoing work with assembling semi-local matrices over patches of the mesh and inserting one row at a time goes well, then we'll use single row insertion anyway. Epetra::FECrsMatrix (Trilinos) has a function to add a single row which is documented to ignore values that aren't part of the sparsity pattern (Epetra::FECrsGraph). Regarding adding sparsity pattern information to an ufc::form, it seems to me it should be part of the *integral classes, one tabulate_sparsity_pattern to match each tabulate_tensor, plus a num_nonzeros() function. > since the only overhead from splitting the rows is the function call. Of > course, if we use a hashed cache then this is definitely a non-issue. > >> >> Of course, if you really care about speed, you could >> >>> implement the momentum equations matrix-free with a custom >> >>> preconditioner >> >>> that >> >>> only uses one block (since they are the same). >> >>> >> >>> Also, it seems like a strange design decision for general coupled >> >>> systems >> >>> for >> >>> the blocks to be separate like this. (Actually, ordering the unknowns >> >>> this way >> >>> makes some factorizations ``work'' even when the matrix is indefinite.) >> >> >> >> Ok. We should look in to this. >> > >> > The reasoning behind this decision was that it makes access to subvectors >> > for components of general hierarchial mixed elements easy. Vector elements >> > and mixed elements are treated as basically the same thing. >> > >> > Thus, renumbering should preferably happen in dolfin, which will have to >> > be >> > done for many occasions anyway. >> >> Yes, it seems natural to have the option to renumber freely in dolfin, in >> particular to switch between "vertex-oriented" and "vector-component >> oriented" numbering. > > I think this would be good. Note that the non-square coupling matrices > between > mixed elements never need to be formed (for most saddle point preconditioners) > so it would be natural to write them in separate form files. > > >> Iterating once is similar to a projection method, but we iterate until >> convergence in the outer non linear residual. Of course, after some time >> the number of iterations is typically 1, so the method collapses into a >> projection type method. > > Okay, cool. Thanks for clarifying. How fast is the convergence the first > time > for the nonlinear problem? My (limited) experience and impression from the > literature is that you only see quadratic convergence from the Newton method > if > you solve the linear saddle point problem instead of replacing it with the > projection preconditioner. Of course, if it only takes one iteration after > some > time, then $\Delta t \eta$ must be small in which case we expect the > projection > method to be good (this is not the case for the slow non-Newtonian fluids I > work > with). > > Jed > > _______________________________________________ > DOLFIN-dev mailing list > [email protected] > http://www.fenics.org/mailman/listinfo/dolfin-dev > > -- Martin _______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
