[petsc-users] PetscFE and TS

[Julian Andrej]

You're absolutely right, i missed the y-direction component of the
momentum equation in the mass matrix integration.

I guess i have most of the parts i need from here to investigate the
remaining on my own!

Thank you so much!

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

> On Mon, Nov 14, 2016 at 7:35 AM, 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.
> >>  
> >
> > There are a few problems here, but first here are my options for
> > the run
> >
> > -vel_petscspace_order 2 -vel_petscspace_poly_tensor 1
> > -pres_petscspace_order 1 -pres_petscspace_poly_tensor 1
> > -ts_dt 0.0001 -ts_max_steps 5
> > -pc_type fieldsplit -pc_fieldsplit_type schur
> > -pc_fieldsplit_schur_precondition full
> >   -fieldsplit_velocity_pc_type lu
> >   -fieldsplit_pressure_ksp_rtol 1.0e-9 -fieldsplit_pressure_pc_type
> > svd -ts_monitor -snes_monitor -snes_view -snes_converged_reason
> > -ksp_converged_reason -ksp_monitor
> >
> > 1) Your g00 term was wrong since it did not handle all velocity
> > components
> >
> > void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
> >   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar
> > u[], const PetscScalar u_t[], const PetscScalar u_x[],
> >   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar
> > a[], const PetscScalar a_t[], const PetscScalar a_x[],
> >   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar
> > g0[]) {
> >   PetscInt d;
> >   for (d = 0; d < dim; ++d) g0[d*dim+d] = u_tShift*1.0;
> > }
> >
> > 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
> >      is not scalable. Eventually, you want to specify the constant
> > null space for that solve. I do it in ex62.
> >
> > 3) I still do not know what is wrong with the exact solution...
> >  
> 
> I found it. You are not making the initial condition. Uncomment the
> line above TSSolve().
> 
>   Thanks,
> 
>     Matt
> 
> 
> >   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  
> >>  
> >
> >
> >
> > --
> > 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
> >  
> 
> 
>