John Peterson wrote:
> On Sun, May 31, 2009 at 12:52 PM, Jed Brown <jed@...> wrote:
>> I'm using libmesh to provide a reference implementation of Q_2 elements
>> for Poisson, and Q_2Q_1 for Stokes. If possible, I would also like to
>> compare with thirdorder elements, but HIERARCHIC appears to produce
>> insanely illconditioned matrices. For example,
>>
>> ./ex4opt d 3 n 4 f HIERARCHIC o THIRD ksp_monitor ksp_converged_reason pc_type lu pc_factor_mat_solver_package superlu
>>
>> does not converge. All of petsc, umfpack, mumps, and spooles are
>> unsuccessful as well. Apparently there is no thirdorder Lagrange
>> implementation so perhaps I'm out of luck here.
>>
>>
>> For the Stokes solver, I modified ex11 for 3D, but no Schur complement
>> preconditioners are available (I could plug it into the new
>> factorization version of PCFieldSplit, but don't have time at the
>> moment). Since BoomerAMG [1] and ML fail when naively given indefinite
>> matrices, the libmesh Stokes implementation isn't really getting a fair
>> shake (I probably just won't include it). I know some people on this
>> list solve (Navier)Stokes, so I'm curious if anyone has an
>> algorithmically scalable implementation they'd like to compare.
>
> Did you modify the boundary conditions from ex11 to account for a
> noninterpolary basis? Looks like we aren't doing the full
> L2projection in ex11... the BC's in ex4 are about what you'd want.
The issue with HIERARCHIC elements occurs for ex4, I never bothered
trying for Stokes. Skipping the L^2 projection just means solving a
slightly different problem (and it won't impact existance of a solution
for the liddriven cavity). In other words, I didn't change this, but I
would be surprised if it made a difference.
> Also, did you pin a value of the pressure (see e.g. ex13)? This
> should not be strictly necessary, but sometimes it helps in practice
> to do it.
It shouldn't make any difference as long as you remove the null space
and you're not trying to use a direct solver. The issue with applying
normal preconditioners to indefinite problems is pretty wellknown.
Interestingly, libmesh's ordering is remarkably good for ILU. For
example, run
./ex11opt ksp_monitor_true_residual ksp_rtol 1e12 pc_type ilu pc_factor_mat_ordering_type rcm
This performs ILU(0) with an ordering computed by reverse CuthillMcKee.
It takes care of the penalty and then stagnates, but with the natural
ordering it converges completely. If you play with various orderings,
you'll find that it's really easy to make ILU fail.
General support for "block factorization" preconditioners is coming soon
in PCFieldSplit. For incompressible flow, you would provide an IS for
all the velocity degrees of freedom and we'll be able to offer a full
suite of preconditioners based on block factorization with various
choices for the Schur complement (e.g. variants of BFBt).
>> Heads up on a change in petscdev: The arguments to MatGetSubMatrix used
>> to reflect implementation details for MPIAIJ/MPIBAIJ and was
>> fundamentally unscalable (though fine for the sizes we've seen so far).
>> We decided to drop the 'csize' argument and take a parallel index set
>> for 'iscol' (instead of the serial gathered one used through 3.0.0).
>> The current implementations still do the gather so there is still a
>> scalability issue when taking submatrices with number of columns similar
>> to singlenode memory, but at least the interface is now scalable (I
>> take submatrices of MatShell so this is good). Speak up if you run into
>> this issue, it can be made scalable but it's not a high priority at the
>> moment.
>
> Thanks for the info! If something is in petscdev, does that mean it
> will go into a patchlevel change or will this take effect only when
> 3.0.1 is released?
Just 3.0.1 (I don't know when that will be). Since I build against
petscdev, I made this oneline change locally.
Jed
