[petsc-users] PetscFE and TS

Julian Andrej juan at tf.uni-kiel.de
Mon Nov 14 07:34:03 CST 2016 

Sorry, i'll make it a habit to include my runtime options ;)

./fluid -dm_refine 2 -vel_petscspace_poly_tensor -vel_petscspace_order 2
-pres_petscspace_poly_tensor -pres_petscspace_order 1 -ts_monitor
-ts_max_steps 10 -ts_dt 0.01 -ts_final_time 2 -snes_monitor
-ts_fd_color

0 TS dt 0.01 time 0.
    0 SNES Function norm 1.843877547310e+01 
    1 SNES Function norm 1.284864233819e-02 
    2 SNES Function norm 7.195271345724e-08 
1 TS dt 0.01 time 0.01 ...

works fine for example. Running with my hand coded jacobian with

./fluid -dm_refine 2 -vel_petscspace_poly_tensor -vel_petscspace_order 2
-pres_petscspace_poly_tensor -pres_petscspace_order 1 -ts_monitor
-ts_max_steps 10 -ts_dt 0.01 -ts_final_time 2 -snes_monitor

0 TS dt 0.01 time 0.
    0 SNES Function norm 1.843877547310e+01 
...
46 SNES Function norm 2.462286148582e-07 
   47 SNES Function norm 1.792848027835e-07 
1 TS dt 0.01 time 0.01 ...

And finally running

./fluid -dm_refine 2 -vel_petscspace_poly_tensor -vel_petscspace_order
2 -pres_petscspace_poly_tensor -pres_petscspace_order 1 -ts_monitor
-ts_max_steps 10 -ts_dt 0.01 -ts_final_time 2 -snes_type test

...
Norm of matrix ratio 0.0203667, difference 3.06481 (user-defined state)
Norm of matrix ratio 0.0203667, difference 3.06481 (constant state -1.0)
Norm of matrix ratio 0.0203667, difference 3.06481 (constant state 1.0)

Thanks again

On Mon, 14 Nov 2016 07:01:09 -0600 Matthew Knepley
<knepley at gmail.com> wrote:

> On Mon, Nov 14, 2016 at 3:24 AM, Julian Andrej <juan at tf.uni-kiel.de>
> wrote:
> 
> > I think i'm understanding the concept now. Although i have another
> > set of questions. If i am formulating a DAE (like in the linear
> > stokes equation for example) it seems i am missing something to
> > formulate the correct jacobian. The residual formulation with
> > snes_mf or snes_fd works fine in terms of nonlinear convergence. If
> > i use the hand coded jacobian the convergence is way slower, which
> > is a reason for an incorrect jacobian formulation from my side (at
> > least what i think right now). I should be able to test my jacobian
> > with -snes_type test, even if i am using TS, right? This gives me a
> > huge difference. I guess i am missing something in the formulation
> > and thus in the implementation of the callbacks here.
> >
> > There is a small example attached.
> >  
> 
> I will look at it. It would be really helpful for me if you would
> include the options you
> run with, and perhaps a little output with the report.
> 
>   Thanks,
> 
>      Matt
> 
> 
> > I appreciate your time and help.
> >
> > On Thu, Nov 10, 2016 at 4:45 PM, Matthew Knepley <knepley at gmail.com>
> > wrote:  
> > > On Thu, Nov 10, 2016 at 9:38 AM, Jed Brown <jed at jedbrown.org>
> > > wrote:  
> > >>
> > >> Matthew Knepley <knepley at gmail.com> writes:
> > >>  
> > >> > On Thu, Nov 10, 2016 at 1:31 AM, Julian Andrej
> > >> > <juan at tf.uni-kiel.de> wrote:
> > >> >  
> > >> >> I'm getting the correct solution now, thanks! There are still
> > >> >> a few open questions.
> > >> >>
> > >> >> 1. The mass term is necessary _and_ using the u_tShift value
> > >> >> if i provide the jacobian. After reading the manual countless
> > >> >> times, i still don't get what the u_tShift value tells me in
> > >> >> this context. Are there any resources which you could point
> > >> >> me to? 
> > >> >
> > >> > The manual is always a lagging resource because we are trying
> > >> > things out.  
> > >>
> > >> This has been explained in the manual similarly to your email for
> > >> several years.  If anyone has a suggestion for how to make it
> > >> better, we'd like to hear.  The typical syndrome is that once
> > >> you learn it, the description suddenly makes sense and you
> > >> wonder why you were ever confused.  It takes feedback from
> > >> people just learning the material to make the explanation more
> > >> clear. 
> > >> > I learned about this by reading examples. The idea is the
> > >> > following. We have an implicit definition of our timestepping
> > >> > method
> > >> >
> > >> >   F(u, grad u, u_t, x, t) = 0  
> > >>
> > >> It doesn't make sense to write grad u as a separate argument
> > >> here.  
> > >
> > >
> > > Sorry, that is just how I think of it, not the actual interface.
> > >  
> > >>
> > >> Also, is `x` meant to be coordinates or some other statically
> > >> prescribed  
> > >
> > >
> > > Coordinates, as distinct from u. Again this is my fault.
> > >  
> > >>
> > >> data?  It can't be actively changing if you expect a generic TS
> > >> to be convergent.  (It's not an argument to the function.)  So
> > >> really, you have:
> > >>
> > >>   F(u, u_t, t) = 0
> > >>  
> > >> > which is not unlike a Lagrangian description of a mechanical
> > >> > system.  
> > The  
> > >> > Jacobian can be considered
> > >> > formally to have two parts,
> > >> >
> > >> >   dF/du and dF/du_t
> > >> >
> > >> > just as in the Lagrangian setting. The u_tShift variable is the
> > >> > multiplier
> > >> > for dF/du_t so that you actually
> > >> > form
> > >> >
> > >> >   dF/du + u_tShift dF/du_t  
> > >>
> > >> We're taking the total derivative of F with respect to u where
> > >> u_t has been _defined_ as an affine function of u, i.e.,
> > >> shift*u+affine. u_tShift comes from the chain rule.  When
> > >> computing the total derivative, we have
> > >>
> > >>   dF/du + dF/du_t du_t/du.
> > >>
> > >>  
> > >> >> 2. Is DMProjectFunction necessary _before_ TSSolve? This acts
> > >> >> like an initial condition if i'm not mistaken.
> > >> >>  
> > >> >
> > >> > Yes, this is the initial condition.
> > >> >
> > >> >  
> > >> >> 3. A more in depth question i guess. Where is the "mass
> > >> >> action" formed? As i could see by stepping through the
> > >> >> debugger, there is  
> > just  
> > >> >> a formed action for the residual equation. It seems way over
> > >> >> my understanding how this formulation works out. I'd also
> > >> >> appreciate resources on that if thats possible.
> > >> >>  
> > >> >
> > >> > Jed is the only one who understands this completely.  
> > >>
> > >> That's a stretch -- several other people have developed
> > >> sophisticated TS implementations.  
> > >
> > >
> > > Matt does not understand this completely.
> > >  
> > >>
> > >>  
> > >> > However, I guess by "mass action" you mean the dF/du_t term.
> > >> > In the explicit methods, you just have u_t + ..., so
> > >> >
> > >> >   dF/du_t = M
> > >> >
> > >> > for finite element methods. So you are putting that there in
> > >> > FormIJacobian(). In the simplest
> > >> > case of backwards Euler then u_tShift would be 1/dt.  
> > >>
> > >> Specifically, for implicit methods, and many semi-implicit
> > >> methods, we don't need M separate.  
> > >
> > >
> > > Yes.
> > >
> > >    Matt
> > >
> > > --
> > > What most experimenters take for granted before they begin their  
> > experiments  
> > > is infinitely more interesting than any results to which their  
> > experiments  
> > > lead.
> > > -- Norbert Wiener  
> >  
> 
> 
>

More information about the petsc-users mailing list