On Fri, Mar 31, 2023 at 4:58 AM Karabelas, Elias ( elias.karabe...@uni-graz.at) <elias.karabe...@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.karabe...@uni-graz.at) <elias.karabe...@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.karabe...@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.karabe...@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! > >
