[petsc-users] PetscFE and TS

[Matthew Knepley]

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
>



-- 
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/e8a93580/attachment.html>