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

[Alexander B Prescott]

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

-- 
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/a08780f3/attachment.html>