<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body>
<div class="moz-cite-prefix">Hi Mark,</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">thanks for the input, however I didn't quite get what you mean.<br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">Maybe I should be a bit more precisce what I want to achieve and why:<br>
<br>
So this specific form of block system arises in some biomedical application that colleagues and me published in
<a class="moz-txt-link-freetext" href="https://www.sciencedirect.com/science/article/pii/S0045782521004230">
https://www.sciencedirect.com/science/article/pii/S0045782521004230</a> (the intersting part is Appendix B.3)<br>
<br>
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<br>
<br>
|F_1(u,p)| | 0 |<br>
| | = | |<br>
|F_2(u,p)| | 0 |<br>
<br>
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). <br>
<br>
After linearization, we arrive exactly at the aforementioned blocksystem, that we solve at the moment by a Schurcomplement approach based on K.<br>
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).<br>
<br>
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.<br>
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).<br>
<br>
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.<br>
<br>
Thanks<br>
Elias<br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">Am 30.03.23 um 19:41 schrieb Mark Adams:<br>
</div>
<blockquote type="cite" cite="mid:CADOhEh60Cd-r7FZHGpLv2VuFxt7C8HMxiLKv-n38V8C5TUYfkQ@mail.gmail.com">
<div dir="ltr">You can lag the update of the Schur complement and use your solver as a preconditioner.
<div>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.<br>
<div><br>
</div>
<div>Mark</div>
</div>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Thu, Mar 30, 2023 at 11:56 AM Karabelas, Elias (<a href="mailto:elias.karabelas@uni-graz.at" target="_blank" moz-do-not-send="true" class="moz-txt-link-freetext">elias.karabelas@uni-graz.at</a>) <<a href="mailto:elias.karabelas@uni-graz.at" target="_blank" moz-do-not-send="true" class="moz-txt-link-freetext">elias.karabelas@uni-graz.at</a>>
wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Dear Community,<br>
<br>
I have a linear system of the form<br>
<br>
|K B| du f1<br>
<br>
=<br>
<br>
|C D| dp f2<br>
<br>
where K is a big m x m sparsematrix that comes from some FE <br>
discretization, B is a coupling matrix (of the form m x 4), C is of the <br>
form (4 x m) and D is 4x4.<br>
<br>
I save B and C as 4 Vecs and D as a 4x4 double array. D might be <br>
singular so at the moment I use the following schur-complement approach <br>
to solve this system<br>
<br>
1) Calculate the vecs v1 = KSP(K,PrecK) * f1, invB = [ KSP(K, PrecK) * <br>
B[0], KSP(K, PrecK) * B[1], KSP(K, PrecK) * B[2], KSP(K, PrecK) * B[3] ]<br>
<br>
2) Build the schurcomplement S=[C ^ K^-1 B - D] via VecDots (C ^ K^-1 B <br>
[i, j] = VecDot(C[i], invB[j])<br>
<br>
3) invert S (this seems to be mostly non-singular) to get dp<br>
<br>
4) calculate du with dp<br>
<br>
So counting this, I need 5x to call KSP which can be expensive and I <br>
thought of somehow doing the Schur-Complement the other way around, <br>
however due to the (possible) singularity of D this seems like a bad <br>
idea (even using a pseudoinverse)<br>
<br>
Two things puzzle me here still<br>
<br>
A) Can this be done more efficiently?<br>
<br>
B) In case my above matrix is the Jacobian in a newton method, how do I <br>
make sure with any form of Schur Complement approach that I hit the <br>
correct residual reduction?<br>
<br>
Thanks<br>
<br>
Elias<br>
<br>
-- <br>
Dr. Elias Karabelas<br>
Universitätsassistent | PostDoc<br>
<br>
Institut für Mathematik & Wissenschaftliches Rechnen | Institute of Mathematics & Scientific Computing<br>
Universität Graz | University of Graz<br>
<br>
Heinrichstraße 36, 8010 Graz<br>
Tel.: +43/(0)316/380-8546<br>
E-Mail: <a href="mailto:elias.karabelas@uni-graz.at" target="_blank" moz-do-not-send="true" class="moz-txt-link-freetext">
elias.karabelas@uni-graz.at</a><br>
Web: <a href="https://ccl.medunigraz.at" rel="noreferrer" target="_blank" moz-do-not-send="true" class="moz-txt-link-freetext">
https://ccl.medunigraz.at</a><br>
<br>
Bitte denken Sie an die Umwelt, bevor Sie dieses E-Mail drucken. Danke!<br>
Please consider the environment before printing this e-mail. Thank you!<br>
<br>
</blockquote>
</div>
</blockquote>
<p><br>
</p>
<pre class="moz-signature" cols="72">--
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: <a class="moz-txt-link-abbreviated" href="mailto:elias.karabelas@uni-graz.at">elias.karabelas@uni-graz.at</a>
Web: <a class="moz-txt-link-freetext" href="https://ccl.medunigraz.at">https://ccl.medunigraz.at</a>
Bitte denken Sie an die Umwelt, bevor Sie dieses E-Mail drucken. Danke!
Please consider the environment before printing this e-mail. Thank you!</pre>
</body>
</html>