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

[Zou, Ling]

Ok, I do not pretend to be shallow water expert here. I only played with it maybe ten years ago.

When you say 5 point stencil in 1d flow, that sounds so much like a finite volume method with TVD reconstruction.
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.

-Ling

From: Alexander B Prescott <alexprescott at email.arizona.edu>
Date: Monday, April 6, 2020 at 8:44 PM
To: Nathan Collier <nathaniel.collier at gmail.com>
Cc: "Zou, Ling" <lzou at anl.gov>, PETSc <petsc-users at mcs.anl.gov>
Subject: Re: [EXT]Re: [petsc-users] Discontinuities in the Jacobian matrix for nonlinear problem

Hi Nathan,

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?

I agree that the difficulty encountered with this sort of problem is counterintuitive to the relative simplification made on nature!

Best,
Alexander

On Mon, Apr 6, 2020 at 12:45 PM Nathan Collier <nathaniel.collier at gmail.com<mailto:nathaniel.collier at gmail.com>> wrote:

External Email
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<mailto: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<mailto:petsc-users-bounces at mcs.anl.gov>> on behalf of Alexander B Prescott <alexprescott at email.arizona.edu<mailto:alexprescott at email.arizona.edu>>
Date: Monday, April 6, 2020 at 1:06 PM
To: PETSc <petsc-users at mcs.anl.gov<mailto: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<mailto:alexprescott at email.arizona.edu>
PhD Candidate, The University of Arizona
Department of Geosciences
1040 E. 4th Street
Tucson, AZ, 85721


--
Alexander Prescott
alexprescott at email.arizona.edu<mailto: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/20200407/3196654b/attachment-0001.html>