<div dir="ltr">Alexander,<div><br></div><div>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. </div><div><br></div><div>I think 'bending' the physics is the right idea, but I never had much luck myself with this. </div><div><br></div><div>Nate</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Apr 6, 2020 at 2:56 PM Zou, Ling via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</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 lang="EN-US">
<div class="gmail-m_-6667456765894780725WordSection1">
<p class="MsoNormal">What about ‘bending’ the physics a bit?<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">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<u></u><u></u></p>
<p class="MsoNormal">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<u></u><u></u></p>
<p class="MsoNormal">Q_i-1/2 proportional to x[i-1] + z[i-1] - (x[i] + z[i]) in between<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">in which, eps is a very small positive number.<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">-Ling<u></u><u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal"><u></u> <u></u></p>
<div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(181,196,223);padding:3pt 0in 0in">
<p class="MsoNormal"><b><span style="font-size:12pt;color:black">From: </span></b><span style="font-size:12pt;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<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">Hello, <u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">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.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">For example, Q_i-1/2 at x[i]:<u></u><u></u></p>
</div>
<div>
<blockquote style="margin-left:30pt;margin-right:0in">
<div>
<p class="MsoNormal">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]<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">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] <u></u><u></u></p>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">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.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">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.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">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.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">Best,<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">Alexander<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal"><br clear="all">
<u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<p class="MsoNormal">-- <u></u><u></u></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"><u></u><u></u></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"><u></u><u></u></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<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:Arial,sans-serif">Department of Geosciences<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:Arial,sans-serif">1040 E. 4th Street<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:9.5pt;font-family:Arial,sans-serif">Tucson, AZ, 85721<u></u><u></u></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote></div>