[petsc-users] Augmented Linear System
Elias Karabelas
elias.karabelas at uni-graz.at
Sat Apr 1 08:26:31 CDT 2023
Dear Mark,
I thought so that all bets would be off with my attempt. So just to
illustrate what I meant by my second approach.
I can start from the second line of the blocksystem and rewrite it to
D dp = f2 - C du
then I tried a pesudoinverse (for lack of alternatives I did this with
an SVD)
and get (denoting the pseudoinverse with D^+)
dp = D^+(f2 - C du)
plugging that into the first line I get
(K - B D^+ C) du = f1 - B D^+ f2
The left hand side I put into a MATSHELL were I store D^+ and the 8
vectors for B and C as a user ctx and update them on the fly during newton.
However as you said using this matrix a Amat and pure K as a Pmat in a
KSP fails more often then it suceeds so I guess this version of the
Schur-Complement is not really well-behaved.
And thanks for the advice with Walker et al, what I meant here was that
in the outer newton I wanted to do an inexact Newton satisfying
|| F(x) + F'(x)dx || <= eta_k || F(x) ||
with F(x) being my residual comprised of the two vectors and F'(x) my
block-system. Now I can select eta_k according to SOME rules, but I
haven't found anything on howto choose the tolerances in my KSP to
finally arrive at this inequality. Reason for doing that in the first
place was, that I wanted to safe iteration counts on my total solves, as
K is non-symmetric and I have to use a GMRES (preferred) or BiCGSTAB
(not so preferred).
Cheers
Elias
Am 31.03.2023 um 15:25 schrieb Mark Adams:
>
>
> On Fri, Mar 31, 2023 at 4:58 AM Karabelas, Elias
> (elias.karabelas at uni-graz.at) <elias.karabelas at uni-graz.at> wrote:
>
> Hi Mark,
>
> thanks for the input, however I didn't quite get what you mean.
>
> Maybe I should be a bit more precisce what I want to achieve and why:
>
> So this specific form of block system arises in some biomedical
> application that colleagues and me published in
> https://www.sciencedirect.com/science/article/pii/S0045782521004230
> (the intersting part is Appendix B.3)
>
> It boils down to a Newton method for solving nolinear mechanics
> describing the motion of the human hear, that is coupled on some
> Neumann surfaces (the surfaces of the inner cavities of each
> bloodpool in the heart) with a pressure that comes from a
> complicated 0D ODE model that describes cardiovascular physiology.
> This comes to look like
>
> |F_1(u,p)| | 0 |
> | | = | |
> |F_2(u,p)| | 0 |
>
> with F_1 is the residual of nonlinear mechanics plus a nonlinear
> boundary coupling term and F_2 is a coupling term to the ODE
> system. In this case u is displacement and p is the pressure
> calculated from the ODE model (one for each cavity in the heart,
> which gives four).
>
> After linearization, we arrive exactly at the aforementioned
> blocksystem, that we solve at the moment by a Schurcomplement
> approach based on K.
> So using this we can get around without doing a MATSHELL for the
> schurcomplement as we can just apply the KSP for K five times in
> order to approximate the solution of the schur-complement system.
>
>
> So you compute an explicit Schur complement (4 solves) and then the
> real solve use 1 more K solve.
> I think this is pretty good as is. You are lucky with only 4 of these
> pressure equations.
> I've actually do this on a problem with 100s of extra equations
> (surface averaging equations) but the problem was linear and would
> 1000s of times steps, or more, and this huge set up cost was amortized.
>
> However here it gets tricky: The outer newton loop is working
> based on an inexact newton method with a forcing term from Walker
> et al. So basically the atol and rtol of the KSP are not constant
> but vary, so I guess this will influence how well we actually
> resolve the solution to the schur-complement system (I tried to
> find some works that can explain how to choose forcing terms in
> this case but found none).
>
>
> Honestly, Walker is a great guy, but I would not get too hung up on this.
> I've done a lot of plasticity work long ago and gave up on Walker et
> al. Others have had the same experience.
> What is new with your problem is how accurately do you want the Schur
> complement (4) solves.
>
> This brought me to think if we can do this the other way around
> and do a pseudo-inverse of D because it's 4x4 and there is no need
> for a KSP there.
> I did a test implementation with a MATSHELL that realizes (K - B
> D^+ C) and use just K for building an GAMG prec however this fails
> spectaculary, because D^+ can behave very badly and the other way
> around I have (C K^-1 B - D) and this behaves more like a singular
> pertubation of the Matrix C K^-1 B which behaves nicer. So here I
> stopped investigating because my PETSc expertise is not bad but
> certainly not good enough to judge which approach would pay off
> more in terms of runtime (by gut feeling was that building a
> MATSHELL requires then only one KSP solve vs the other 5).
>
> However I'm happy to hear some alternatives that I could persue in
> order to speed up our current solver strategy or even be able to
> build a nice MATSHELL.
>
>
> OK, so you have tried what I was alluding to.
> I don't follow what you did exactly and have not worked it out, but
> there should be an iteration on the pressure equation with a (lagged)
> Schur solve as a preconditioner.
> But with only 4 extra solves in your case, I don't think it is worth
> it unless you want to write solver papers.
> And AMG in general really has to have a normal PDE, e.g. the K^-1
> solve, and if K is too far away from the Laplacian (or elasticity)
> then all bets are off.
> Good luck,
> Mark
>
>
> Thanks
> Elias
>
>
>
>
> Am 30.03.23 um 19:41 schrieb Mark Adams:
>> You can lag the update of the Schur complement and use your
>> solver as a preconditioner.
>> If your problems don't change much you might converge fast enough
>> (ie, < 4 iterations with one solve per iteration), but what you
>> have is not bad if the size of your auxiliary, p, space does not
>> grow.
>>
>> Mark
>>
>> On Thu, Mar 30, 2023 at 11:56 AM Karabelas, Elias
>> (elias.karabelas at uni-graz.at) <elias.karabelas at uni-graz.at> wrote:
>>
>> Dear Community,
>>
>> I have a linear system of the form
>>
>> |K B| du f1
>>
>> =
>>
>> |C D| dp f2
>>
>> where K is a big m x m sparsematrix that comes from some FE
>> discretization, B is a coupling matrix (of the form m x 4), C
>> is of the
>> form (4 x m) and D is 4x4.
>>
>> I save B and C as 4 Vecs and D as a 4x4 double array. D might be
>> singular so at the moment I use the following
>> schur-complement approach
>> to solve this system
>>
>> 1) Calculate the vecs v1 = KSP(K,PrecK) * f1, invB = [ KSP(K,
>> PrecK) *
>> B[0], KSP(K, PrecK) * B[1], KSP(K, PrecK) * B[2], KSP(K,
>> PrecK) * B[3] ]
>>
>> 2) Build the schurcomplement S=[C ^ K^-1 B - D] via VecDots
>> (C ^ K^-1 B
>> [i, j] = VecDot(C[i], invB[j])
>>
>> 3) invert S (this seems to be mostly non-singular) to get dp
>>
>> 4) calculate du with dp
>>
>> So counting this, I need 5x to call KSP which can be
>> expensive and I
>> thought of somehow doing the Schur-Complement the other way
>> around,
>> however due to the (possible) singularity of D this seems
>> like a bad
>> idea (even using a pseudoinverse)
>>
>> Two things puzzle me here still
>>
>> A) Can this be done more efficiently?
>>
>> B) In case my above matrix is the Jacobian in a newton
>> method, how do I
>> make sure with any form of Schur Complement approach that I
>> hit the
>> correct residual reduction?
>>
>> Thanks
>>
>> Elias
>>
>> --
>> Dr. Elias Karabelas
>> Universitätsassistent | PostDoc
>>
>> Institut für Mathematik & Wissenschaftliches Rechnen |
>> Institute of Mathematics & Scientific Computing
>> Universität Graz | University of Graz
>>
>> Heinrichstraße 36, 8010 Graz
>> Tel.: +43/(0)316/380-8546
>> E-Mail: elias.karabelas at uni-graz.at
>> Web: https://ccl.medunigraz.at
>>
>> Bitte denken Sie an die Umwelt, bevor Sie dieses E-Mail
>> drucken. Danke!
>> Please consider the environment before printing this e-mail.
>> Thank you!
>>
>
> --
> Dr. Elias Karabelas
> Universitätsassistent | PostDoc
>
> Institut für Mathematik & Wissenschaftliches Rechnen | Institute of Mathematics & Scientific Computing
> Universität Graz | University of Graz
>
> Heinrichstraße 36, 8010 Graz
> Tel.: +43/(0)316/380-8546
> E-Mail:elias.karabelas at uni-graz.at
> Web: https://ccl.medunigraz.at
>
> Bitte denken Sie an die Umwelt, bevor Sie dieses E-Mail drucken. Danke!
> Please consider the environment before printing this e-mail. Thank you!
>
--
Dr. Elias Karabelas
Universitätsassistent | PostDoc
Institut für Mathematik & Wissenschaftliches Rechnen | Institute of Mathematics & Scientific Computing
Universität Graz | University of Graz
Heinrichstraße 36, 8010 Graz
Tel.: +43/(0)316/380-8546
E-Mail:elias.karabelas at uni-graz.at
Web:https://ccl.medunigraz.at
Bitte denken Sie an die Umwelt, bevor Sie dieses E-Mail drucken. Danke!
Please consider the environment before printing this e-mail. Thank you!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230401/8616b9a8/attachment-0001.html>
More information about the petsc-users
mailing list