[petsc-users] PetscFE and TS

Matthew Knepley knepley at gmail.com
Thu Nov 10 06:15:37 CST 2016 

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

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


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

  Thanks,

     Matt


> Thanks again.
>
> Julian
>
> On Wed, Nov 9, 2016 at 9:25 PM, Matthew Knepley <knepley at gmail.com> wrote:
> > On Wed, Nov 9, 2016 at 8:14 AM, Matthew Knepley <knepley at gmail.com>
> wrote:
> >>
> >> On Wed, Nov 9, 2016 at 7:52 AM, Julian Andrej <juan at tf.uni-kiel.de>
> wrote:
> >>>
> >>> I tried the example with different signs for the forcing, so this
> >>> might be a mistake by me, yes. it doesn't change anything though. I
> >>> also have to correct myself, the time for the callback in
> >>> PetscDSAddBoundary changes correctly! Although u_t in all other
> >>> callbacks seems to remain zero.
> >>>
> >>> I am out of ideas on my side though :)
> >>
> >>
> >> The bug is that the boundary conditions are not being properly
> enforced. I
> >> am fixing it now.
> >> Thanks for the simple example demonstrating this. It made it much much
> >> easier to debug.
> >
> >
> > I think this works. Let me know if you agree.
> >
> > There is the business of u_t on the boundary, which is not handled right
> > now. I am putting that
> > in, but the homogeneous values for that is still alright here.
> >
> >   Thanks,
> >
> >      Matt
> >
> >>
> >>     Matt
> >>
> >>>
> >>> On Wed, Nov 9, 2016 at 2:33 PM, Matthew Knepley <knepley at gmail.com>
> >>> wrote:
> >>> > On Wed, Nov 9, 2016 at 12:32 AM, Julian Andrej <juan at tf.uni-kiel.de>
> >>> > wrote:
> >>> >>
> >>> >> I tested a few other things and it turned out that the function
> added
> >>> >> as a dirichlet boundary condition via PetscDSAddBoundary doesn't
> >>> >> receive the time values correctly (PetscReal time is always 0.0).
> >>> >
> >>> >
> >>> > I will check out your example and fix this. However, I wonder if
> there
> >>> > is a
> >>> > bug in your implementation:
> >>> >
> >>> > void f0_temp(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, const PetscReal x[], PetscScalar f0[])
> >>> > {
> >>> >   f0[0] = u_t[0] - 2.0;
> >>> > }
> >>> >
> >>> > I get
> >>> >
> >>> >   du/dt - \Delta u = -2
> >>> >
> >>> > so that if u = x^2 + y^2 + 2t, then
> >>> >
> >>> >   2 - (2 + 2) = -2
> >>> >
> >>> > so you would want u_t + 2.
> >>> >
> >>> > If you give the Laplacian the opposite sign, this is an unstable
> >>> > formulation.
> >>> >
> >>> >   Thanks,
> >>> >
> >>> >      Matt
> >>> >
> >>> >>
> >>> >> I also see no effect adding DMTSSetIJacobianLocal vs. not setting
> the
> >>> >> jacobian function. The analytic solution i used in my attached
> example
> >>> >> is the same as in the src/ts/examples/tutorials/ex32.c.
> >>> >>
> >>> >> Am i doing anything obviously wrong here? My next step would be to
> try
> >>> >> if assembling the matrices myself and writing a custom
> >>> >> IFunction/IJacobian which assembles the different parts of the
> matrix
> >>> >> like M and J for the stiff ODE with nontrivial mass matrix (see
> manual
> >>> >> section of TS -> F = M u' - f) but i think this should be obsolete
> >>> >> right?
> >>> >>
> >>> >> Appreciate your help.
> >>> >>
> >>> >> Regards
> >>> >> Julian
> >>> >>
> >>> >> On Mon, Nov 7, 2016 at 10:49 AM, Julian Andrej <juan at tf.uni-kiel.de
> >
> >>> >> wrote:
> >>> >> > Hello,
> >>> >> >
> >>> >> > i'm using PetscFE in combination with SNES and a hand written
> >>> >> > backward
> >>> >> > euler solver right now and like to hand that over to TS. I have
> been
> >>> >> > fiddling quite a while getting TS to work with PetscFE and have
> >>> >> > encountered a few unclarities.
> >>> >> >
> >>> >> > I have looked at examples which are outdated i guess
> >>> >> > (src/ts/examples/tutorials/ex32.c) and i am confused on how i
> have
> >>> >> > to
> >>> >> > formulate the discretization of the jacobian and the residual. It
> >>> >> > seems obvious to use the DMTSSetIFunctionLocal and
> >>> >> > DMTSSetIJacobianLocal functions because there are prepared
> wrappers
> >>> >> > from DMPlexTSComputeIFunctionFEM and DMPlexTSComputeIJacobianFEM.
> >>> >> >
> >>> >> > I need a Mass matrix for my problem, is that formed from the
> >>> >> > function
> >>> >> > space information or do i have to form it myself? Is there any
> >>> >> > working
> >>> >> > example which uses PetscFE and TS to form the
> DMTSSetIFunctionLocal
> >>> >> > and DMTSSetIJacobianLocal? (grepping the tree tells me "no"
> >>> >> > unfortunately).
> >>> >> >
> >>> >> > Regards
> >>> >> > Julian
> >>> >
> >>> >
> >>> >
> >>> >
> >>> > --
> >>> > 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
>



-- 
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/20161110/1c1fe01c/attachment.html>

More information about the petsc-users mailing list