<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<div class="moz-cite-prefix">Dear Mark,</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">I thought so that all bets would be off
with my attempt. So just to illustrate what I meant by my second
approach.</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">I can start from the second line of the
blocksystem and rewrite it to</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">D dp = f2 - C du</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">then I tried a pesudoinverse (for lack
of alternatives I did this with an SVD)</div>
<div class="moz-cite-prefix">and get (denoting the pseudoinverse
with D^+)<br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">dp = D^+(f2 - C du)</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">plugging that into the first line I get</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">(K - B D^+ C) du = f1 - B D^+ f2</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">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.</div>
<div class="moz-cite-prefix">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.</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">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</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">|| F(x) + F'(x)dx || <= eta_k ||
F(x) ||</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">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).</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">Cheers</div>
<div class="moz-cite-prefix">Elias<br>
</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">Am 31.03.2023 um 15:25 schrieb Mark
Adams:<br>
</div>
<blockquote type="cite"
cite="mid:CADOhEh5eUMJkej6-yCCnun-mMt4basK85uiUvx2FT_65YMu74Q@mail.gmail.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<div dir="ltr">
<div dir="ltr"><br>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Fri, Mar 31, 2023 at
4:58 AM Karabelas, Elias (<a
href="mailto:elias.karabelas@uni-graz.at"
moz-do-not-send="true" class="moz-txt-link-freetext">elias.karabelas@uni-graz.at</a>)
<<a href="mailto:elias.karabelas@uni-graz.at"
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">
<div>
<div>Hi Mark,</div>
<div><br>
</div>
<div>thanks for the input, however I didn't quite get what
you mean.<br>
</div>
<div><br>
</div>
<div>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
href="https://www.sciencedirect.com/science/article/pii/S0045782521004230"
target="_blank" moz-do-not-send="true"
class="moz-txt-link-freetext">
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. </div>
</div>
</blockquote>
<div><br>
</div>
<div>So you compute an explicit Schur complement (4 solves)
and then the real solve use 1 more K solve.</div>
<div>I think this is pretty good as is. You are lucky with
only 4 of these pressure equations.</div>
<div>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.</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div>
<div>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>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>Honestly, Walker is a great guy, but I would not get too
hung up on this.</div>
<div>I've done a lot of plasticity work long ago and gave up
on Walker et al. Others have had the same experience.</div>
<div>What is new with your problem is how accurately do you
want the Schur complement (4) solves.</div>
<div> <br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div>
<div>
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>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>OK, so you have tried what I was alluding to.</div>
<div>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.</div>
<div>But with only 4 extra solves in your case, I don't think
it is worth it unless you want to write solver papers.</div>
<div>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.</div>
<div> </div>
<div>Good luck,</div>
<div>Mark</div>
<div><br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div>
<div>
<br>
Thanks<br>
Elias<br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div>Am 30.03.23 um 19:41 schrieb Mark Adams:<br>
</div>
<blockquote type="cite">
<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 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 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>
Web: <a href="https://ccl.medunigraz.at" target="_blank" moz-do-not-send="true" class="moz-txt-link-freetext">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>
</div>
</blockquote>
</div>
</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>