[petsc-users] TS question 1: how to stop explicit methods because they do not use SNES(VI)?

[Ed Bueler]

Jed --

> My point was that including sliding potentially adds another
> stiff/algebraic term so whatever interface we choose better be able to
> support at least two stiff terms.

Yes, totally agreed w.r.t. the interface design.  (The DAE that you are
solving in that case is stiff.)

However, *my* observation about the current interface design is actually
that it was never apparently tested with TS+SNESVI and is awkward in that
case.

IMHO the design for time-stepping on constrained problems should continue
to allow use of SNESVI (in implicit and imex cases) if an LHSfunction (or
similar) is provided, *and* have a mechanism where explicit time-solvers
can use projection onto the bounds---conceivably this could be a
trivialized SNESVI instance---if bounds are present.  This is because mere
evaluation of LHSFunction and RHSFunction will not otherwise generate the
"side effect" of enforcing the constraints which were provided by
FormBounds() in using the SNESVI.

In both the current and proposed designs in this thread, for implicit and
imex methods the constraint enforcement is inside SNESVI *but* for explicit
time-stepping it would be a user-written projection in TSPostStep().
Furthermore the LHSFunction might have to be rewritten to allow
non-feasible inputs just to support explicit multistage time-stepping.  Yuk.

Note that there are VIs in L^2 and they *are* "projection onto the bounds"
in the just-mentioned sense.  See section II.3 of Kinderlehrer and
Stampacchia.  (Elliptic VIs in W^{1,p} are what one usually sees.)  The
switch from VI in Sobolev space to VI in L^2 is exactly what you see when
you semi-discretize these obstacle-like problems in time only and then
compare implicit and explicit timestepping, respectively.

So ya'll are having a good argument ... but I am worried about what I asked
about anyway.

>> q_total = - Gamma H^{n+2} |grad s|^{n-1} grad s + V H
> Is V really a prescribed function or the solution of an equation that
> involves u?

I know that the problem you want me to solve is a big DAE including both
mass and momentum conservation.  (Do you really think that all this time I
haven't noticed?)

My interest is in good solvers for what you call the "local problem".  This
local problem is the one that limits PISM hybrid runs *and* I think it is
why Stokes is not really functional yet in anyone's time-stepping model.
 (If they treat the time-stepping explicitly then they have the SIA-type
time-step restriction kick in as a solver convergence problem, because the
SIA tells the truth in the majority of the area where the bed is sticky.)
 Thus my test case code does not waste time generating a real V (i.e.
solving simultaneously for H and V at each time step, for the DAE you want
me to solve ...).

> For the SIA thickness equation, I think it would be interesting to
> consider the conductivity tensor in the Newton linearization.
>   d_u [-div (u^k |grad u|^{p-2} grad u)] * w
>     = -div (k w^{k-1} |grad u|^{p-2} grad u)
>       -div (u^k [|grad u|^{p-2} I + (p-2) |grad u|^{p-4} grad u \otimes
grad u] grad w)
> The first part is transport of the perturbation w.  The next is
> diffusion of w with an anisotropic conductivity tensor.

Yes I have considered it.  See http://dx.doi.org/10.1017/jog.2015.3, both
the discussion of flux-splitting and the appendix on analytical jacobian.
What do *you* get by considering this conductivity tensor?

>> So now I'll emphasize what I said before: it is *degenerate* p-bratu.
> To my knowledge, "p-Bratu" is just something I invented to combine two
> nonlinearities.  You don't have an e^u nonlinearity so there is no
> Bratu-type nonlinearity.

There *is* a term I want to treat explicitly, and it is *not* related to
the stiffness etc., but to the coupling to the atmosphere.  (I said this
before.)  Namely, in my form u_t - f(t,u)=g(t,u) the RHS g(t,u) is an
elevation dependent mass balance term. Actually it can be (and might as
well be) bratu i.e. exponentialish in the thickness.  Mostly it is out of
my control because it is the altitude dependence of precipitation in the
bottom few km of the atmosphere arising from coupling to atmosphere
models.  My test case linearizes it a la a glaciologist, with an ELA and a
lapse rate.

By analogy, in other words, the time-dependent form of the
I-know-Jed-invented-it p-bratu problem is pretty representative of my
problem ... at least if you look in p>2 cases, add additional degeneracy to
the diffusivity, make the zero-order term signed, use SNESVI (because the
source term is signed and so a VI/NCP must be solved), and you add a
sliding flux.  I was trying to work by analogy ...

Ed


On Wed, Feb 15, 2017 at 9:37 AM, Jed Brown <jed at jedbrown.org> wrote:

> Ed Bueler <elbueler at alaska.edu> writes:
>
> > Jed --
> >
> >>> u_t - div (u^k |grad u|^{p-2} grad u) = g(t,u)
> >>Are you also interested in the case with sliding?
> >
> > You must be getting grumpy.  There is this 2009 paper about sliding ...
> > Yes I am interested in sliding but I was trying to keep things simple for
> > you.  ;-)
>
> My point was that including sliding potentially adds another
> stiff/algebraic term so whatever interface we choose better be able to
> support at least two stiff terms.
>
> > The code in question is an under-development example in my book that I am
> > trying to use to show off SNESVI functionality in an interesting case and
> > yet not get too nasty.  See the help string at the top of the file for
> what
> > it actually does:
> >
> > https://github.com/bueler/p4pdes/blob/master/c/ch11/ice.c
> >
> > In summary, the actual flux is
> >
> > q_total = - Gamma H^{n+2} |grad s|^{n-1} grad s + V H
>
> Is V really a prescribed function or the solution of an equation that
> involves u?
>
> > and it is all handled implicitly by preference, i.e. in the current API I
> > have IFunction = u_t - div q_total.
> >
> > The real interest of the problem is the failure of long-range
> communication
> > of low-(spatial)-frequency information when you hit the free boundary,
> and
> > these barriers can emerge anywhere in the domain.  (In a mountain range
> you
> > can have neighboring small glaciers completely out of phase; how should
> > multigrid ... or any other preconditioning one would hope for ... work
> for
> > all the glaciers in a mountain range?
>
> This already exists for variable coefficient problems -- the electric
> potential in a twisted wire pair is independent despite being nearby
> points in the solution of a Poisson equation.
>
> An important difference from ice sheets is that the Green's functions
> are highly nonlocal for the twisted wire pair, but very local for ice
> sheets with sticky beds.  With slippery beds (ice shelves are the
> limiting case), the velocity Green's functions become similar to
> Poisson.  That is the scenario where multigrid is needed.  Otherwise the
> solutions are quite local.
>
> For the SIA thickness equation, I think it would be interesting to
> consider the conductivity tensor in the Newton linearization.
>
>   d_u [-div (u^k |grad u|^{p-2} grad u)] * w
>     = -div (k w^{k-1} |grad u|^{p-2} grad u)
>       -div (u^k [|grad u|^{p-2} I + (p-2) |grad u|^{p-4} grad u \otimes
> grad u] grad w)
>
> The first part is transport of the perturbation w.  The next is
> diffusion of w with an anisotropic conductivity tensor.
>
> > So now I'll emphasize what I said before: it is *degenerate* p-bratu.
>
> To my knowledge, "p-Bratu" is just something I invented to combine two
> nonlinearities.  You don't have an e^u nonlinearity so there is no
> Bratu-type nonlinearity.
>
> > Doubly-degenerate because p>2, but the "H^{n+2}" aka "u^k" is the bad
> > thing that makes the abrupt free boundary.)  Sliding, in this regard,
> > is the easy case, so I am throwing it in by generating an artificial
> > velocity field; yes this is the coupling that should include the SSA
> > in the future ... we already have an SSA solver or two that behave
> > surprisingly well, you know, as long as the boundary conditions can be
> > figured-out.
>
> I think sliding is the only phenomena that produces a nonlocal response.
> Consider an ice shelf.  The thickness and slope at the grounding line
> instantaneously affect the seaward edge of the shelf.  Grounded SIA has
> no such dynamics.
>
> >> I would just have a PetscOptionsEList() for whether to evaluate the
> >> diffusive term on the left or right.  Shouldn't be more than a few lines
> >> of code.
> >
> > So I have something like?:
> >   -ice_diffusive_side left,right
> >
> > O.k.  What if user says
> >   -ts_type rk -ice_diffusive_side left
> > Then I get the current brokenness: code runs and produces nonsense
> because
> > RK does not call SNESVI (constraint never enforced) *and* diffusive stuff
> > is never evaluated (because RK does not call IFunction).
>
> Yes, I think TSRK should error if IFunction has been specified.
>
> > So maybe then I also go through possible -ts_types and check whether the
> > option combination made sense?  Or move some terms to the correct side
> > based on what that -ts_type should want?  Do I include -ts_type sundials
> > because somebody's config might include that third-party thing (even
> though
> > mine doesn't)?
>
> I am not suggesting switching on TS type.  I guess I felt that having
> dependent options was acceptable in exchange for it being easy/reliable
> to know which terms are being treated explicitly or implicitly.  If the
> method takes over that choice, it makes the interface "smarter" which is
> not necessarily a good thing.  It requires copious documentation and
> diagnostics so the user can tell what the hell is happening under the
> hood.  But I see the utility.
>
> > That's what I mean by fragile code.
> >
> > Such things build up!  Each time you shove cool stuff in my face I have
> to
> > get used to seeing ugly stuff turn into beautiful stuff as my aesthetic
> > sense follows my growing understanding.  But it takes time (learning
> curve)
> > and all the time it is fragile code (i.e. for me ... because I don't
> > appreciate all the motivations behind that part of the API ...)
> >
> > On the other hand ... I agree that if I am convinced this is really a
> good
> > API for something I want, then indeed I will suck it up and learn it and
> > write it.
> >
> > Recall my original request was for the ability to know if the TS type I
> > have chosen will actually use the SNES I have set ... so long ago it now
> > seems.
>
> I think this is just how design conversations work.
>



-- 
Ed Bueler
Dept of Math and Stat and Geophysical Institute
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
301C Chapman and 410D Elvey
907 474-7693 and 907 474-7199  (fax 907 474-5394)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170215/7a10c6f8/attachment.html>