<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body dir="ltr">
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Hi Matt, </div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
What you suggested in your last email was exactly what I did on my very first attempt at the problem, and while it "worked," convergence was not satisfactory due to the Newton step being fixed in step 6. This is the reason I would like to use the
<span style="font-family: Calibri, Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255); display: inline !important">
linesearch in </span>SNES instead. Indeed in your manual you "recommend most PETSc users work directly with SNES, rather than using PETSc for the linear problem within a nonlinear solver." Ideally I'd like to create a SNES solver, pass in the functions to
evaluate <b>K</b>, <b>F,</b> <b>a,</b> and k, and set up the underlying KSP object as in my first message. Is this possible?<br>
</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Thanks,</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
<div id="Signature">
<div style="font-family:Tahoma; font-size:13px">
<div>
<div><font style="font-size:13px"><font><b style="color:rgb(34,34,34); font-size:small; background-color:rgb(255,255,255)">Alexander (Olek) Niewiarowski</b>
</font></font>
<div style="color:rgb(34,34,34); background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">PhD Candidate, Civil & Environmental Engineering</font></div>
<div style="color:rgb(34,34,34); background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">Princeton University, 2020</font></div>
<div style="color:rgb(34,34,34); background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">Cell: +1 (610) 393-2978</font></div>
</div>
</div>
</div>
</div>
<div id="appendonsend"></div>
<hr style="display:inline-block;width:98%" tabindex="-1">
<div id="divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size:11pt" color="#000000"><b>From:</b> Matthew Knepley <knepley@gmail.com><br>
<b>Sent:</b> Thursday, February 6, 2020 5:33<br>
<b>To:</b> Olek Niewiarowski <aan2@princeton.edu><br>
<b>Cc:</b> Smith, Barry F. <bsmith@mcs.anl.gov>; petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?</font>
<div> </div>
</div>
<div>
<div dir="ltr">
<div dir="ltr">On Wed, Feb 5, 2020 at 8:53 PM Olek Niewiarowski <<a href="mailto:aan2@princeton.edu">aan2@princeton.edu</a>> wrote:<br>
</div>
<div class="x_gmail_quote">
<blockquote class="x_gmail_quote" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr">
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<span style="margin:0px; font-family:Arial,Helvetica,sans-serif">Hi Barry and Matt, </span></div>
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<br>
</div>
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<span style="margin:0px; font-family:Arial,Helvetica,sans-serif">Thank you for your input and for creating a new issue in the repo. </span></div>
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<span style="margin:0px; font-family:Arial,Helvetica,sans-serif">My initial question was more basic (how to configure the SNES's KSP solver as in my first message with
<b>a</b> and <b>k</b>), but now I see there's more to the implementation. To reiterate, f</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">or my problem's structure, a good solution algorithm (on the algebraic level) is the following "double
back-substitution": </span></div>
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<span style="margin:0px; font-family:Arial,Helvetica,sans-serif">For each nonlinear iteration:</span></div>
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<ol>
<li><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">define intermediate vectors u_1 and u_2 </span></li><li><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">solve Ku_1 = -F --> u_1 = -K^{-1}F (this solve is cheap, don't actually need K^-1)</span></li><li><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">solve Ku_2 = -a --> u_2 = </span><span style="margin:0px"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; background-color:rgb(255,255,255); display:inline">-K^{-1}a </span><span style="margin:0px; background-color:rgb(255,255,255); display:inline"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; background-color:rgb(255,255,255); display:inline">(ditto)</span></span></span></li><li><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">define \beta = 1/(1 + k a^Tu_2)</span></li><li><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">\Delta u = u_1 + \beta k u_2^T F u_2</span><br>
</li><li><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">u = u + \Delta u</span></li></ol>
</div>
</div>
</div>
</blockquote>
<div>This is very easy to setup: </div>
<div><br>
</div>
<div> 1) Create a KSP object KSPCreate(comm, &ksp)</div>
<div><br>
</div>
<div> 2) Call KSPSetOperators(ksp, K, K,)</div>
<div><br>
</div>
<div> 3) Solve the first system KSPSolve(ksp, -F, u_1)</div>
<div><br>
</div>
<div> 4) Solve the second system KSPSolve(ksp, a, u_2)</div>
<div><br>
</div>
<div> 5) Calculate beta VecDot(a, u_2, &gamma); beta = 1./(1. + k*gamma);</div>
<div><br>
</div>
<div> 6) Update the guess VecDot(u_2, F, &delta); VecAXPBYPCZ(u, 1.0, beta*delta, 1.0, u_1, u_2)</div>
<div><br>
</div>
<div>Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div><br>
</div>
<blockquote class="x_gmail_quote" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr">
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<div style="margin:0px; font-size:12pt; font-family:Calibri,Arial,Helvetica,sans-serif">
<div style="margin:0px"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">I don't need the Jacobian inverse, </span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)">[</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><b>K</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)">−</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)">k</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><b>aa</b></span><sup style="color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255); font-family:MJXc-TeX-main-B,MJXc-TeX-main-Bx,MJXc-TeX-main-Bw"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><i>T</i></span></sup><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)">]</span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">-1</span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"> </span></sup></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)">=</span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"> </span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><b>K</b></span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">-1 </span></sup></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"> -
(k</span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>K</b></span><sup style="font-family:Arial,Helvetica,sans-serif; background-color:white"><span style="margin:0px; font-size:12pt">-1</span></sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"> </span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>a</b></span><span style="margin:0px; background-color:white"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>a</b></span><sup style="font-family:MJXc-TeX-main-B,MJXc-TeX-main-Bx,MJXc-TeX-main-Bw; background-color:white"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><i>T</i></span></sup></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>K</b></span><sup style="font-family:Arial,Helvetica,sans-serif; background-color:white"><span style="margin:0px; font-size:12pt">-1</span></sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white">)</span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)">/(1+k</span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; color:rgb(34,34,34); word-spacing:normal; background-color:rgb(255,255,255)"><span style="margin:0px; background-color:white"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>a</b></span><sup style="font-family:MJXc-TeX-main-B,MJXc-TeX-main-Bx,MJXc-TeX-main-Bw; background-color:white"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><i>T</i></span></sup></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>K</b></span><sup style="font-family:Arial,Helvetica,sans-serif; background-color:white"><span style="margin:0px; font-size:12pt">-1</span></sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white"><b>a</b></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif; background-color:white">) </span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">just
the solution </span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">Δ</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"><b>u =</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"> </span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">[</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"><b>K</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">−</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">k</span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"><b>aa</b></span><sup style="word-spacing:normal; font-family:MJXc-TeX-main-B,MJXc-TeX-main-Bx,MJXc-TeX-main-Bw"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><i>T</i></span></sup><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">]</span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; word-spacing:normal"><sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">-1</span></sup></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"><b>F<span> </span></b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">= </span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"><b>K</b></span><span style="margin:0px; font-size:11pt; font-family:Arial,Helvetica,sans-serif,serif,EmojiFont; word-spacing:normal"><sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">-1</span></sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>F</b></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"> -
(k</span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>K</b></span><sup style="font-family:Arial,Helvetica,sans-serif"><span style="margin:0px; font-size:12pt">-1</span></sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"> </span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>a</b></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b><b>F</b></b></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>K</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><sup><span style="margin:0px; font-size:12pt">-1</span></sup></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>a</b></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">)</span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">/(1
+ k</span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>a</b></span><span style="margin:0px"><sup style="font-family:MJXc-TeX-main-B,MJXc-TeX-main-Bx,MJXc-TeX-main-Bw"><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><i>T</i></span></sup></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>K</b></span><sup style="font-family:Arial,Helvetica,sans-serif"><span style="margin:0px; font-size:12pt">-1</span></sup><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif"><b>a</b></span><span style="margin:0px; font-size:12pt; font-family:Arial,Helvetica,sans-serif">) </span></span></div>
<div style="margin:0px"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">=<span> </span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><b style="word-spacing:normal">u</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">_1
+ beta k<span> </span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><b style="word-spacing:normal">u</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">_2^T<span> </span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><b style="word-spacing:normal">F
u</b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">_2 (so I never need to invert<span> </span></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif"><b style="word-spacing:normal">K<span> </span></b></span><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal">either). </span><span style="font-family:Arial,Helvetica,sans-serif; word-spacing:normal">(To
Matt's point on gitlab, K is a symmetric sparse matrix arising from a bilinear form. ) Also, e<span style="font-family:Arial,Helvetica,sans-serif; background-color:rgb(255,255,255); display:inline">ventually, I want to have more than one low-rank updates to
K, but again, Sherman Morrisson Woodbury should still work. </span></span></div>
<div style="margin:0px"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif; word-spacing:normal"><br>
</span></div>
<div style="margin:0px"><span style="margin:0px; font-family:Arial,Helvetica,sans-serif">Being new to PETSc, I don't know if this algorithm directly translates into an efficient numerical solver. (I'm also not sure if Picard iteration would be useful here.)
What would it take to set up the KSP solver in SNES like I did below? <span style="font-family:Arial,Helvetica,sans-serif; background-color:rgb(255,255,255); display:inline">Is it possible<span> "out of the box"?</span></span></span><span style="font-family:Arial,Helvetica,sans-serif">
I looked at MatCreateLRC() - what would I pass this to? (A pointer to demo/tutorial would be very appreciated.) If there's a better way to go about all of this, I'm open to any and all ideas. My only limitation is that I use petsc4py exclusively since I/future
users of my code will not be comfortable with C.</span></div>
<div style="margin:0px"><span style="font-family:Arial,Helvetica,sans-serif"><br>
</span></div>
<div style="margin:0px"><span style="font-family:Arial,Helvetica,sans-serif">Thanks again for your help!</span></div>
<div style="margin:0px"><span style="font-family:Arial,Helvetica,sans-serif"><br>
</span></div>
</div>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<br>
</div>
<div id="x_gmail-m_-8841679181880611122Signature">
<div style="font-family:Tahoma; font-size:13px">
<div>
<div><font style="font-size:13px"><font><b style="color:rgb(34,34,34); font-size:small; background-color:rgb(255,255,255)">Alexander (Olek) Niewiarowski</b>
</font></font>
<div style="color:rgb(34,34,34); background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">PhD Candidate, Civil & Environmental Engineering</font></div>
<div style="color:rgb(34,34,34); background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">Princeton University, 2020</font></div>
<div style="color:rgb(34,34,34); background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">Cell: +1 (610) 393-2978</font></div>
</div>
</div>
</div>
</div>
<div id="x_gmail-m_-8841679181880611122appendonsend"></div>
<hr style="display:inline-block; width:98%">
<div id="x_gmail-m_-8841679181880611122divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" color="#000000" style="font-size:11pt"><b>From:</b> Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>><br>
<b>Sent:</b> Wednesday, February 5, 2020 15:46<br>
<b>To:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
<b>Cc:</b> Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>>;
<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b> Re: [petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?</font>
<div> </div>
</div>
<div><font size="2"><span style="font-size:11pt">
<div><br>
<a href="https://gitlab.com/petsc/petsc/issues/557" target="_blank">https://gitlab.com/petsc/petsc/issues/557</a><br>
<br>
<br>
> On Feb 5, 2020, at 7:35 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
> <br>
> Perhaps Barry is right that you want Picard, but suppose you really want Newton.<br>
> <br>
> "This problem can be solved efficiently using the Sherman-Morrison formula" Well, maybe. The main assumption here is that inverting K is cheap. I see two things you can do in a straightforward way:<br>
> <br>
> 1) Use MatCreateLRC() <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateLRC.html" target="_blank">
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateLRC.html</a> to create the Jacobian<br>
> and solve using an iterative method. If you pass just K was the preconditioning matrix, you can use common PCs.<br>
> <br>
> 2) We only implemented MatMult() for LRC, but you could stick your SMW code in for MatSolve_LRC if you really want to factor K. We would<br>
> of course help you do this.<br>
> <br>
> Thanks,<br>
> <br>
> Matt<br>
> <br>
> On Wed, Feb 5, 2020 at 1:36 AM Smith, Barry F. via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
> <br>
> I am not sure of everything in your email but it sounds like you want to use a "Picard" iteration to solve [K(u)−kaaT]Δu=−F(u). That is solve<br>
> <br>
> A(u^{n}) (u^{n+1} - u^{n}) = F(u^{n}) - A(u^{n})u^{n} where A(u) = K(u) - kaaT<br>
> <br>
> PETSc provides code to this with SNESSetPicard() (see the manual pages) I don't know if Petsc4py has bindings for this.<br>
> <br>
> Adding missing python bindings is not terribly difficult and you may be able to do it yourself if this is the approach you want.<br>
> <br>
> Barry<br>
> <br>
> <br>
> <br>
> > On Feb 4, 2020, at 5:07 PM, Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>> wrote:<br>
> > <br>
> > Hello, <br>
> > I am a FEniCS user but new to petsc4py. I am trying to modify the KSP solver through the SNES object to implement the Sherman-Morrison formula(e.g.
<a href="http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html" target="_blank">
http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html</a> ). I am solving a nonlinear system of the form [K(u)−kaaT]Δu=−F(u). Here the jacobian matrix K is modified by the term kaaT, where k is a scalar. Notably, K is sparse, while the term kaaT results
in a full matrix. This problem can be solved efficiently using the Sherman-Morrison formula :
<br>
> > <br>
> > [K−kaaT]-1 = K-1 - (kK-1 aaTK-1)/(1+kaTK-1a)<br>
> > I have managed to successfully implement this at the linear solve level (by modifying the KSP solver) inside a custom Newton solver in python by following an incomplete tutorial at
<a href="https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner" target="_blank">
https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner</a> :<br>
> > • while (norm(delU) > alpha): # while not converged<br>
> > • <br>
> > • self.update_F() # call to method to update r.h.s form<br>
> > • self.update_K() # call to update the jacobian form<br>
> > • K = assemble(self.K) # assemble the jacobian matrix<br>
> > • F = assemble(self.F) # assemble the r.h.s vector<br>
> > • a = assemble(self.a_form) # assemble the a_form (see Sherman Morrison formula)<br>
> > • <br>
> > • for bc in self.mem.bc: # apply boundary conditions<br>
> > • bc.apply(K, F) <br>
> > • bc.apply(K, a) <br>
> > • <br>
> > • B = PETSc.Mat().create() <br>
> > • <br>
> > • # Assemble the bilinear form that defines A and get the concrete
<br>
> > • # PETSc matrix <br>
> > • A = as_backend_type(K).mat() # get the PETSc objects for K and a<br>
> > • u = as_backend_type(a).vec() <br>
> > • <br>
> > • # Build the matrix "context" # see firedrake docs<br>
> > • Bctx = MatrixFreeB(A, u, u, self.k) <br>
> > • <br>
> > • # Set up B <br>
> > • # B is the same size as A <br>
> > • B.setSizes(*A.getSizes()) <br>
> > • B.setType(B.Type.PYTHON) <br>
> > • B.setPythonContext(Bctx) <br>
> > • B.setUp() <br>
> > • <br>
> > • <br>
> > • ksp = PETSc.KSP().create() # create the KSP linear solver object<br>
> > • ksp.setOperators(B) <br>
> > • ksp.setUp() <br>
> > • pc = ksp.pc <br>
> > • pc.setType(pc.Type.PYTHON) <br>
> > • pc.setPythonContext(MatrixFreePC()) <br>
> > • ksp.setFromOptions() <br>
> > • <br>
> > • solution = delU # the incremental displacement at this iteration<br>
> > • <br>
> > • b = as_backend_type(-F).vec() <br>
> > • delu = solution.vector().vec() <br>
> > • <br>
> > • ksp.solve(b, delu) <br>
> > <br>
> > • self.mem.u.vector().axpy(0.25, self.delU.vector()) # poor man's linesearch<br>
> > • counter += 1 <br>
> > Here is the corresponding petsc4py code adapted from the firedrake docs:<br>
> > <br>
> > • class MatrixFreePC(object): <br>
> > • <br>
> > • def setUp(self, pc): <br>
> > • B, P = pc.getOperators() <br>
> > • # extract the MatrixFreeB object from B <br>
> > • ctx = B.getPythonContext() <br>
> > • self.A = ctx.A <br>
> > • self.u = ctx.u <br>
> > • self.v = ctx.v <br>
> > • self.k = ctx.k <br>
> > • # Here we build the PC object that uses the concrete, <br>
> > • # assembled matrix A. We will use this to apply the action <br>
> > • # of A^{-1} <br>
> > • self.pc = PETSc.PC().create() <br>
> > • self.pc.setOptionsPrefix("mf_") <br>
> > • self.pc.setOperators(self.A) <br>
> > • self.pc.setFromOptions() <br>
> > • # Since u and v do not change, we can build the denominator <br>
> > • # and the action of A^{-1} on u only once, in the setup <br>
> > • # phase. <br>
> > • tmp = self.A.createVecLeft() <br>
> > • self.pc.apply(self.u, tmp) <br>
> > • self._Ainvu = tmp <br>
> > • self._denom = 1 + self.k*self.v.dot(self._Ainvu) <br>
> > • <br>
> > • def apply(self, pc, x, y): <br>
> > • # y <- A^{-1}x <br>
> > • self.pc.apply(x, y) <br>
> > • # alpha <- (v^T A^{-1} x) / (1 + v^T A^{-1} u) <br>
> > • alpha = (self.k*self.v.dot(y)) / self._denom <br>
> > • # y <- y - alpha * A^{-1}u <br>
> > • y.axpy(-alpha, self._Ainvu) <br>
> > • <br>
> > • <br>
> > • class MatrixFreeB(object): <br>
> > • <br>
> > • def __init__(self, A, u, v, k): <br>
> > • self.A = A <br>
> > • self.u = u <br>
> > • self.v = v <br>
> > • self.k = k <br>
> > • <br>
> > • def mult(self, mat, x, y): <br>
> > • # y <- A x <br>
> > • self.A.mult(x, y) <br>
> > • <br>
> > • # alpha <- v^T x <br>
> > • alpha = self.v.dot(x) <br>
> > • <br>
> > • # y <- y + alpha*u <br>
> > • y.axpy(alpha, self.u) <br>
> > However, this approach is not efficient as it requires many iterations due to the Newton step being fixed, so I would like to implement it using SNES and use line search. Unfortunately, I have not been able to find any documentation/tutorial on how to do
so. Provided I have the FEniCS forms for F, K, and a, I'd like to do something along the lines of:<br>
> > solver = PETScSNESSolver() # the FEniCS SNES wrapper<br>
> > snes = solver.snes() # the petsc4py SNES object<br>
> > ## ??<br>
> > ksp = snes.getKSP()<br>
> > # set ksp option similar to above<br>
> > solver.solve()<br>
> > <br>
> > I would be very grateful if anyone could could help or point me to a reference or demo that does something similar (or maybe a completely different way of solving the problem!).
<br>
> > Many thanks in advance!<br>
> > Alex<br>
> <br>
> <br>
> <br>
> -- <br>
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> -- Norbert Wiener<br>
> <br>
> <a href="https://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
<br>
</div>
</span></font></div>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr" class="x_gmail_signature">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener</div>
<div><br>
</div>
<div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>