[petsc-users] Block preconditioning for 3d problem

[Jed Brown]

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.
>>>
>>