# [petsc-users] Apply operator to linearised system before solving?

Asbjørn Nilsen Riseth riseth at maths.ox.ac.uk
Fri Jun 12 06:22:34 CDT 2015

Thank you Matt,

I'll read through the paper and see what I can map onto my problem.
The CPR method is also about isolating the pressure equation, so SIMPLE
looks relevant.

Ozzy

On Fri, 12 Jun 2015 at 04:54 Matthew Knepley <knepley at gmail.com> wrote:

> On Wed, Jun 10, 2015 at 7:23 PM, Asbjørn Nilsen Riseth <
> riseth at maths.ox.ac.uk> wrote:
>
>> Hi Matt,
>>
>> 1) K =
>> [               I                   0 ;
>> -A10*inv(diag(A00))     I  ]
>>
>
> Okay, it looks like this is some sort of weird approximate Schur
> complement. The right thing to
> do I think is to get everything written out in linear algebra language,
> since I think you can just
> use FieldSplit+PCComposite to do this from the command line. This one
> looks similar to SIMPLE.
> The paper by Shuttleworth, Elman, Shadid, etc. does a nice job of this
> kind of classification. It
> might help in the process.
>
>
>> There is a typo in my first email: S = Atilde11 = A11 -
>> A10*inv(diag(A00))*A01
>>
>> 2) I'd like to have the linear solver close to what these people
>> currently use. My ultimate goal is to look at nonlinear solver
>> strategies, not improve their linear solver. I already have a PC that seems
>> to work fine for my purposes, but was hoping I could implement this version
>> of CPR to get even closer to their current systems.
>>
>
> Okay.
>
>    Matt
>
>
>> Ozzy
>>
>> On Thu, 11 Jun 2015 at 00:54 Matthew Knepley <knepley at gmail.com> wrote:
>>
>>> On Wed, Jun 10, 2015 at 6:30 PM, Asbjørn Nilsen Riseth <
>>> riseth at maths.ox.ac.uk> wrote:
>>>
>>>> Dear PETSc community,
>>>>
>>>> I'm trying to implement a preconditioner used in reservoir modelling,
>>>> called Constrained Pressure Residual.
>>>> They apply a transformation to the linearised system we get from each
>>>> Newton step, before solving it. Then a 2-stage multiplicative
>>>> preconditioner on the transformed system.
>>>>
>>>
>>> 1) What is K?
>>>
>>> 2) Do you care about the 2-stage thing, or would any Stokes solver do?
>>>
>>>   Thanks,
>>>
>>>     Matt
>>>
>>>
>>>> The main problem is implementing step 1 below.
>>>>
>>>> Let A be the Jacobian, b the residual and K my transformation.
>>>> A = [A00 A01; A10 A11].
>>>> The process is roughly like this:
>>>> 1) Set Atilde = KA,  btilde = Kb.
>>>>  - We want to solve Atilde x =  btilde
>>>>
>>>> 2) Create a 2-stage multiplicative preconditioner using Atilde, btilde
>>>> pc0: This is only applied to the "fieldsplit 1" block of my system.
>>>> B_1 = [0 0; 0 S^-1]
>>>> Where S is a selfp Schur approximation from Atilde
>>>> S =  Atilde11 - Atilde10 * inv(diag(Atilde00))* Atilde01
>>>> pc1: This is a standard ILU on the whole system Atilde
>>>>
>>>> Currently I'm doing something like this
>>>> Step 1:
>>>> -ksp_type richardson -ksp_max_it 1
>>>> -pc_type python
>>>> Then I create Atilde from KA in PCSetup, and a FGMRES ksp to take care
>>>> of step 2 with PCApply
>>>>
>>>> Step 2:
>>>> pc_type composite
>>>> pc_composite_type multiplicative
>>>> pc_composite_pcs python,ilu,
>>>>
>>>> Are there better ways of dealing with this transformation?
>>>> To me it looks similar to a 2-step right preconditioner on top of a
>>>> left preconditioner.
>>>>
>>>> Regards,
>>>> Ozzy
>>>>
>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> -- Norbert Wiener
>>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their