<html xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.msonormal0, li.msonormal0, div.msonormal0
{mso-style-name:msonormal;
mso-margin-top-alt:auto;
margin-right:0in;
mso-margin-bottom-alt:auto;
margin-left:0in;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
span.EmailStyle20
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style>
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal">Ok, I do not pretend to be shallow water expert here. I only played with it maybe ten years ago.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">When you say 5 point stencil in 1d flow, that sounds so much like a finite volume method with TVD reconstruction.<o:p></o:p></p>
<p class="MsoNormal">I solved many 1d flow equation using PETSc (Euler eqn., single-phase flow eqn., two-phase flow eqn.), had a lot of ‘zero-flow’ conditions, so far PETSc works pretty well for me.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">-Ling<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b><span style="font-size:12.0pt;color:black">From: </span></b><span style="font-size:12.0pt;color:black">Alexander B Prescott <alexprescott@email.arizona.edu><br>
<b>Date: </b>Monday, April 6, 2020 at 8:44 PM<br>
<b>To: </b>Nathan Collier <nathaniel.collier@gmail.com><br>
<b>Cc: </b>"Zou, Ling" <lzou@anl.gov>, PETSc <petsc-users@mcs.anl.gov><br>
<b>Subject: </b>Re: [EXT]Re: [petsc-users] Discontinuities in the Jacobian matrix for nonlinear problem<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Hi Nathan, <o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Thanks for the thoughtful response. For the 1D toy version I force flow to occur at every node, but my ultimate objective is to apply this to a 2D floodplain in which a node may or may not have flow in the 'true' solution. Based on what
you point out, however, it seems that this approach may be doomed to fail. For one thing, using a 5 point stencil, a node with no flow and whose neighbor's have no flow will produce a row in the 'bent physics' Jacobian of all zero entries. Based on my experience
so far that will cause PETSc to return errors. Does that sound right to you?<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">I agree that the difficulty encountered with this sort of problem is counterintuitive to the relative simplification made on nature!<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Best,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal">Alexander<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal">On Mon, Apr 6, 2020 at 12:45 PM Nathan Collier <<a href="mailto:nathaniel.collier@gmail.com">nathaniel.collier@gmail.com</a>> wrote:<o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<p align="center" style="text-align:center"><strong><span style="font-family:"Calibri",sans-serif;color:red">External Email</span></strong><o:p></o:p></p>
<div>
<p class="MsoNormal">Alexander, <o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">I am not familiar with your specific model, but I do have experience working with the diffusive / kinematic wave approximation to shallow water equations. I always had a lot of trouble near ponding conditions (steady state), which is odd
because when nature says do nothing, you don't expect the equations to be hard to solve. Some kind soul pointed out to me that this was because the equations are derived using a power law model relating flow to slope/depth which breaks down when you should
have no flow. So as long as you have flow, the solver behaves well and things are fine (i.e. smooth terrain, academic test problems), but when you have local ponding due to real terrain effects, you are hosed. A better solution approach doesn't help if your
equations aren't meant for the problem you are trying to solve. <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">I think 'bending' the physics is the right idea, but I never had much luck myself with this. <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Nate<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal">On Mon, Apr 6, 2020 at 2:56 PM Zou, Ling via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<o:p></o:p></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">What about ‘bending’ the physics a bit?<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Q_i-1/2 proportional to sqrt(x[i-1] + z[i-1] - (x[i] + z[i])), if x[i-1]+z[i-1] > x[i]+z[i] + eps<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Q_i-1/2 proportional to -1.0*sqrt(x[i] + z[i] - (x[i-1] + z[i-1])), if x[i]+z[i] > x[i-1]+z[i-1] + eps<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Q_i-1/2 proportional to x[i-1] + z[i-1] - (x[i] + z[i]) in between<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">in which, eps is a very small positive number.<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">-Ling<o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><b><span style="font-size:12.0pt;color:black">From:
</span></b><span style="font-size:12.0pt;color:black">petsc-users <<a href="mailto:petsc-users-bounces@mcs.anl.gov" target="_blank">petsc-users-bounces@mcs.anl.gov</a>> on behalf of Alexander B Prescott <<a href="mailto:alexprescott@email.arizona.edu" target="_blank">alexprescott@email.arizona.edu</a>><br>
<b>Date: </b>Monday, April 6, 2020 at 1:06 PM<br>
<b>To: </b>PETSc <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject: </b>[petsc-users] Discontinuities in the Jacobian matrix for nonlinear problem</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Hello,
<o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">The non-linear boundary-value problem I am applying PETSc to is a relatively simple steady-state flow routing algorithm based on the continuity equation, such that Div(Q) = 0 everywhere
(Q=discharge). I use a finite volume approach to calculate flow between nodes, with Q calculated as a piecewise smooth function of the local flow depth and the water-surface slope. In 1D, the residual is calculated as R(x_i)=Q_i-1/2 - Q_i+1/2.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">For example, Q_i-1/2 at x[i]:<o:p></o:p></p>
</div>
<div>
<blockquote style="margin-left:30.0pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Q_i-1/2 proportional to sqrt(x[i-1] + z[i-1] - (x[i] + z[i])), if x[i-1]+z[i-1] > x[i]+z[i]<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Q_i-1/2 proportional to -1.0*sqrt(x[i] + z[i] - (x[i-1] + z[i-1])), if x[i]+z[i] > x[i-1]+z[i-1] <o:p></o:p></p>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Where z[i] is local topography and doesn't change over the iterations, and Q_i+1/2 is computed analogously. So the residual derivatives with respect to x[i-1], x[i] and x[i+1] are
not continuous when the water-surface slope = 0.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Are there intelligent ways to handle this problem? My 1D trial runs naively fix any zero-valued water-surface slopes to a small non-zero positive value (e.g. 1e-12). Solver convergence
has been mixed and highly dependent on the initial guess. So far, FAS with QN coarse solver has been the most robust.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Restricting x[i] to be non-negative is a separate issue, to which I have applied the SNES_VI solvers. They perform modestly but have been less robust.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Best,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">Alexander<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><br clear="all">
<o:p></o:p></p>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"> <o:p></o:p></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto">--
<o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">Alexander Prescott</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:9.5pt;font-family:"Arial",sans-serif"><a href="mailto:alexprescott@email.arizona.edu" target="_blank">alexprescott@email.arizona.edu</a></span><o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">PhD Candidate, The University of Arizona</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">Department of Geosciences</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">1040 E. 4th Street</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">Tucson, AZ, 85721</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</blockquote>
</div>
<p class="MsoNormal"><br clear="all">
<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal">-- <o:p></o:p></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">Alexander Prescott</span><span style="font-size:9.5pt"><o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Arial",sans-serif"><a href="mailto:alexprescott@email.arizona.edu" target="_blank">alexprescott@email.arizona.edu</a></span><span style="font-size:9.5pt"><o:p></o:p></span></p>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">PhD Candidate, The University of Arizona<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">Department of Geosciences<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">1040 E. 4th Street<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Arial",sans-serif">Tucson, AZ, 85721<o:p></o:p></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>