<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Mon, Nov 14, 2016 at 7:35 AM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span class="">On Mon, Nov 14, 2016 at 3:24 AM, Julian Andrej <span dir="ltr"><<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">I think i'm understanding the concept now. Although i have another set<br>
of questions. If i am formulating a DAE (like in the linear stokes<br>
equation for example) it seems i am missing something to formulate the<br>
correct jacobian. The residual formulation with snes_mf or snes_fd<br>
works fine in terms of nonlinear convergence. If i use the hand coded<br>
jacobian the convergence is way slower, which is a reason for an<br>
incorrect jacobian formulation from my side (at least what i think<br>
right now). I should be able to test my jacobian with -snes_type test,<br>
even if i am using TS, right? This gives me a huge difference. I guess<br>
i am missing something in the formulation and thus in the<br>
implementation of the callbacks here.<br>
<br>
There is a small example attached.<br></blockquote><div><br></div></span><div>There are a few problems here, but first here are my options for the run </div><div><br></div><div>-vel_petscspace_order 2 -vel_petscspace_poly_tensor 1 -pres_petscspace_order 1 -pres_petscspace_poly_tensor 1<br></div><div><div>-ts_dt 0.0001 -ts_max_steps 5<br></div><div>-pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_<wbr>precondition full</div><div>  -fieldsplit_velocity_pc_type lu</div><div>  -fieldsplit_pressure_ksp_rtol 1.0e-9 -fieldsplit_pressure_pc_type svd</div></div><div>-ts_monitor -snes_monitor -snes_view -snes_converged_<wbr>reason -ksp_converged_reason -ksp_<wbr>monitor</div><div><br></div><div>1) Your g00 term was wrong since it did not handle all velocity components</div><div><br></div><div><div>void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,</div><span class=""><div><span class="m_-6317961402219212361gmail-Apple-tab-span" style="white-space:pre-wrap">       </span>   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],</div><div><span class="m_-6317961402219212361gmail-Apple-tab-span" style="white-space:pre-wrap">        </span>   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],</div></span><div><span class="m_-6317961402219212361gmail-Apple-tab-span" style="white-space:pre-wrap">   </span>   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])</div><div>{</div><div>  PetscInt d;</div><div>  for (d = 0; d < dim; ++d) g0[d*dim+d] = u_tShift*1.0;</div><div>}</div></div><div><br></div><div>2) You do not handle the null space for the pressure solve right now. The easy way is to use SVD for the PC, however this</div><div>     is not scalable. Eventually, you want to specify the constant null space for that solve. I do it in ex62.</div><div><br></div><div>3) I still do not know what is wrong with the exact solution...</div></div></div></div></blockquote><div><br></div><div>I found it. You are not making the initial condition. Uncomment the line above TSSolve().</div><div><br></div><div>  Thanks,</div><div><br></div><div>    Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>  Thanks,</div><div><br></div><div>     Matt</div><div><div class="h5"><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
I appreciate your time and help.<br>
<div class="m_-6317961402219212361gmail-HOEnZb"><div class="m_-6317961402219212361gmail-h5"><br>
On Thu, Nov 10, 2016 at 4:45 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
> On Thu, Nov 10, 2016 at 9:38 AM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
>><br>
>> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> writes:<br>
>><br>
>> > On Thu, Nov 10, 2016 at 1:31 AM, Julian Andrej <<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a>><br>
>> > wrote:<br>
>> ><br>
>> >> I'm getting the correct solution now, thanks! There are still a few<br>
>> >> open questions.<br>
>> >><br>
>> >> 1. The mass term is necessary _and_ using the u_tShift value if i<br>
>> >> provide the jacobian. After reading the manual countless times, i<br>
>> >> still don't get what the u_tShift value tells me in this context. Are<br>
>> >> there any resources which you could point me to?<br>
>> >><br>
>> ><br>
>> > The manual is always a lagging resource because we are trying things<br>
>> > out.<br>
>><br>
>> This has been explained in the manual similarly to your email for<br>
>> several years.  If anyone has a suggestion for how to make it better,<br>
>> we'd like to hear.  The typical syndrome is that once you learn it, the<br>
>> description suddenly makes sense and you wonder why you were ever<br>
>> confused.  It takes feedback from people just learning the material to<br>
>> make the explanation more clear.<br>
>><br>
>> > I learned about this by reading examples. The idea is the<br>
>> > following. We have an implicit definition of our timestepping method<br>
>> ><br>
>> >   F(u, grad u, u_t, x, t) = 0<br>
>><br>
>> It doesn't make sense to write grad u as a separate argument here.<br>
><br>
><br>
> Sorry, that is just how I think of it, not the actual interface.<br>
><br>
>><br>
>> Also, is `x` meant to be coordinates or some other statically prescribed<br>
><br>
><br>
> Coordinates, as distinct from u. Again this is my fault.<br>
><br>
>><br>
>> data?  It can't be actively changing if you expect a generic TS to be<br>
>> convergent.  (It's not an argument to the function.)  So really, you<br>
>> have:<br>
>><br>
>>   F(u, u_t, t) = 0<br>
>><br>
>> > which is not unlike a Lagrangian description of a mechanical system. The<br>
>> > Jacobian can be considered<br>
>> > formally to have two parts,<br>
>> ><br>
>> >   dF/du and dF/du_t<br>
>> ><br>
>> > just as in the Lagrangian setting. The u_tShift variable is the<br>
>> > multiplier<br>
>> > for dF/du_t so that you actually<br>
>> > form<br>
>> ><br>
>> >   dF/du + u_tShift dF/du_t<br>
>><br>
>> We're taking the total derivative of F with respect to u where u_t has<br>
>> been _defined_ as an affine function of u, i.e., shift*u+affine.<br>
>> u_tShift comes from the chain rule.  When computing the total<br>
>> derivative, we have<br>
>><br>
>>   dF/du + dF/du_t du_t/du.<br>
>><br>
>><br>
>> >> 2. Is DMProjectFunction necessary _before_ TSSolve? This acts like an<br>
>> >> initial condition if i'm not mistaken.<br>
>> >><br>
>> ><br>
>> > Yes, this is the initial condition.<br>
>> ><br>
>> ><br>
>> >> 3. A more in depth question i guess. Where is the "mass action"<br>
>> >> formed? As i could see by stepping through the debugger, there is just<br>
>> >> a formed action for the residual equation. It seems way over my<br>
>> >> understanding how this formulation works out. I'd also appreciate<br>
>> >> resources on that if thats possible.<br>
>> >><br>
>> ><br>
>> > Jed is the only one who understands this completely.<br>
>><br>
>> That's a stretch -- several other people have developed sophisticated TS<br>
>> implementations.<br>
><br>
><br>
> Matt does not understand this completely.<br>
><br>
>><br>
>><br>
>> > However, I guess by "mass action" you mean the dF/du_t term. In the<br>
>> > explicit methods, you just have u_t + ..., so<br>
>> ><br>
>> >   dF/du_t = M<br>
>> ><br>
>> > for finite element methods. So you are putting that there in<br>
>> > FormIJacobian(). In the simplest<br>
>> > case of backwards Euler then u_tShift would be 1/dt.<br>
>><br>
>> Specifically, for implicit methods, and many semi-implicit methods, we<br>
>> don't need M separate.<br>
><br>
><br>
> Yes.<br>
><br>
>    Matt<br>
><br>
> --<br>
> What most experimenters take for granted before they begin their experiments<br>
> is infinitely more interesting than any results to which their experiments<br>
> lead.<br>
> -- Norbert Wiener<br>
</div></div></blockquote></div></div></div><div><div class="h5"><br><br clear="all"><div><br></div>-- <br><div class="m_-6317961402219212361gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>