[petsc-users] Discontinuities in the Jacobian matrix for nonlinear problem

[Alexander B Prescott]

Hello,

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.
For example, Q_i-1/2 at x[i]:

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]
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]


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.

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.

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.

Best,
Alexander



-- 
Alexander Prescott
alexprescott at email.arizona.edu
PhD Candidate, The University of Arizona
Department of Geosciences
1040 E. 4th Street
Tucson, AZ, 85721
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200406/e0a9d703/attachment.html>