[petsc-users] Augmented Linear System

Fri Mar 31 03:58:31 CDT 2023

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

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.

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<mailto:elias.karabelas at uni-graz.at>) <elias.karabelas at uni-graz.at<mailto: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<mailto: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<mailto: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/20230331/844e8d0e/attachment-0001.html>