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

[Nathan Collier]

Alexander,

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.

I think 'bending' the physics is the right idea, but I never had much luck
myself with this.

Nate






On Mon, Apr 6, 2020 at 2:56 PM Zou, Ling via petsc-users <
petsc-users at mcs.anl.gov> wrote:

> What about ‘bending’ the physics a bit?
>
>
>
> 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
>
> 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
>
> Q_i-1/2 proportional to x[i-1] + z[i-1] - (x[i] +
> z[i])                               in between
>
>
>
> in which, eps is a very small positive number.
>
>
>
> -Ling
>
>
>
>
>
> *From: *petsc-users <petsc-users-bounces at mcs.anl.gov> on behalf of
> Alexander B Prescott <alexprescott at email.arizona.edu>
> *Date: *Monday, April 6, 2020 at 1:06 PM
> *To: *PETSc <petsc-users at mcs.anl.gov>
> *Subject: *[petsc-users] Discontinuities in the Jacobian matrix for
> nonlinear problem
>
>
>
> 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/fdb45196/attachment.html>