[petsc-users] Block preconditioning for 3d problem

[Jed Brown]

Dave Lee <davelee2804 at gmail.com> writes:

> Hey Jed,
>
> Right now they are uncoupled purely for testing purposes. I just want to
> get the uncoupled problem right before moving on to the coupled problem.

Okay, well you may not want to change your dof ordering to put slabs
together since you'll likely have some parts of your code that do
vertical processing (column integrals/physics or PDE stiffness in the
vertical).  It's common for climate/weather codes to use some horizontal
index (within a blocking scheme or an entire subdomain array index) as
the innermost to facilitate vectorization of the column processing
(which has high enough arithmetic intensity that vectorization matters).

> When I add the vertical dynamics back in both the outer nonlinear
> convergence and the inner linear convergence are roughly the same as for
> when I solve for the 2D slabs only as a 3D system.

Which you've said is terrible compared to separate 2D slabs, so we have
to get to the bottom of that.  We outlined a number of tests you can do
to understand why it isn't converging as expected, but haven't heard
back on the results so it'll be hard to advise further.

> I like your scaling idea. I tried using the pressure multiplied by the mass
> matrix as a preconditioner for the inner linear problem, but this didn't
> help. Perhaps some sort of scaling is the way to go though.
>
> Cheers, Dave.
>
> On Fri, Oct 11, 2019 at 2:45 PM Jed Brown <jed at jedbrown.org> wrote:
>
>> Why are your slabs decoupled at present?  (Have you done a transform in
>> the vertical?)  Is the linear convergence significantly different when
>> you include the multiple layers?
>>
>> Dave Lee <davelee2804 at gmail.com> writes:
>>
>> > Hi Jed and Mark,
>> >
>> > thanks for your helpful comments. Yes the nonlinear outer problem is
>> > uncoupled between the slabs, it is only the linear inner problem where
>> they
>> > are coupled.
>> >
>> > I've tried to make the slab DOFs close in memory, and also tried using a
>> > tight tolerance on the outer KSP (1.0e-20), but without success.
>> >
>> > You make some really good points about scaling issues Jed and Mark. This
>> is
>> > a solve for a global atmosphere, and each 2D slab is a horizontal layer
>> of
>> > the atmosphere, so the pressure (which the linear solve is for) will vary
>> > dramatically between slabs. Perhaps I can additionally precondition the
>> > linear problem to normalise the pressure in each slab so that they stay
>> > close to unity.
>> >
>> > Cheers, Dave.
>> >
>> > On Fri, Oct 11, 2019 at 2:52 AM Mark Adams <mfadams at lbl.gov> wrote:
>> >
>> >> I can think of a few sources of coupling in the solver: 1) line search
>> and
>> >> 2) Krylov, and 3) the residual test (scaling issues). You could turn
>> >> linesearch off and use Richardson (with a fixed number of iterations) or
>> >> exact solves as Jed suggested. As far as scaling can you use the same NL
>> >> problem on each slab? This should fix all the problems anyway. Or, on
>> the
>> >> good decouple solves, if the true residual is of the same scale and
>> *all*
>> >> of the slabs converge well then you should be OK on scaling. If this
>> works
>> >> then start adding stuff back in and see what breaks it.
>> >>
>> >> On Thu, Oct 10, 2019 at 11:01 AM Jed Brown via petsc-users <
>> >> petsc-users at mcs.anl.gov> wrote:
>> >>
>> >>> Dave Lee via petsc-users <petsc-users at mcs.anl.gov> writes:
>> >>>
>> >>> > Hi PETSc,
>> >>> >
>> >>> > I have a nonlinear 3D problem for a set of uncoupled 2D slabs.
>> (Which I
>> >>> > ultimately want to couple once this problem is solved).
>> >>> >
>> >>> > When I solve the inner linear problem for each of these 2D slabs
>> >>> > individually (using KSPGMRES) the convergence of the outer nonlinear
>> >>> > problem is good. However when I solve the inner linear problem as a
>> >>> single
>> >>> > 3D problem (with no coupling between the 2D slabs, so the matrix is
>> >>> > effectively a set of uncoupled blocks, one for each 2D slab) the
>> outer
>> >>> > nonlinear convergence degrades dramatically.
>> >>>
>> >>> Is the nonlinear problem also decoupled between slabs?
>> >>>
>> >>> If you solve the linear problem accurately (tight tolerances on the
>> >>> outer KSP, or global direct solve), is the outer nonlinear convergence
>> >>> good again?  If not, test that your Jacobian is correct (it likely
>> isn't
>> >>> or you have inconsistent scaling leading to ill conditioning).  SNES
>> has
>> >>> automatic tests for that, but you aren't using SNES so you'd have to
>> >>> write some code.
>> >>>
>> >>> What happens if you run the 2D problem (where convergence is currently
>> >>> good) with much smaller subdomains (or -pc_type pbjacobi)?
>> >>>
>> >>> > Note that I am not using SNES, just my own quasi-Newton approach for
>> the
>> >>> > outer nonlinear problem.
>> >>> >
>> >>> > I suspect that the way to recover the convergence for the 3D coupled
>> >>> > problem is to use some sort of PCBJACOBI or PCFIELDSPLIT
>> preconditioner,
>> >>> > but I'm not entirely sure. I've tried following this example:
>> >>> >
>> >>>
>> https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex7.c.html
>> >>> > but without any improvement in the convergence.
>> >>> >
>> >>> > On each processor the slab index is the major index and the internal
>> >>> slab
>> >>> > DOF index is the minor index. The parallel decomposition is within
>> the
>> >>> 2D
>> >>> > slab dimensions only, not between slabs.
>> >>>
>> >>> For convergence, you usually want the direction of tight coupling
>> >>> (sounds like that is within slabs) to be close in memory.
>> >>>
>> >>> In general, use -ksp_monitor_true_residual -ksp_converged_reason.
>> >>>
>> >>> > I'm configuring the 3D inner linear solver as
>> >>> >
>> >>> > KSPGetPC(ksp, &pc);
>> >>> > PCSetType(pc, PCBJACOBI);
>> >>> > PCBJacobiSetLocalBlocks(pc, num_slabs, NULL);
>> >>> > KSPSetUp(ksp);
>> >>> > PCBJacobiGetSubKSP(pc, &nlocal, &first_local, &subksp);
>> >>> > for(int ii = 0; ii < nlocal; ii++) {
>> >>> >     KSPGetPC(subksp[ii], &subpc);
>> >>> >     PCSetType(subpc, PCJACOBI);
>> >>> >     KSPSetType(subksp[ii], KSPGMRES);
>> >>> >     KSPSetTolerances(subksp[ii], 1.e-12, PETSC_DEFAULT,
>> PETSC_DEFAULT,
>> >>> > PETSC_DEFAULT);
>> >>> > }
>> >>> >
>> >>> > However this does not seem to change the outer nonlinear convergence.
>> >>> When
>> >>> > I run with ksp_view the local block sizes look correct - number of
>> >>> local 2D
>> >>> > slab DOFs for each block on each processor.
>> >>> >
>> >>> > Any ideas on what sort of KSP/PC configuration I should be using to
>> >>> recover
>> >>> > the convergence of the uncoupled 2D slab linear solve for this 3D
>> linear
>> >>> > solve? Should I be rethinking my node indexing to help with this?
>> >>> >
>> >>> > Cheers, Dave.
>> >>>
>> >>
>>