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