[petsc-users] PetscFE and TS

[Matthew Knepley]

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...

  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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161114/5b2ea8bb/attachment-0001.html>