[petsc-users] Questions about TS

[Jed Brown]

On Wed, Oct 19, 2011 at 12:32, Blaise Bourdin <bourdin at lsu.edu> wrote:

> Hi,
>
>
> > On Wed, Oct 19, 2011 at 10:54, Blaise Bourdin <bourdin at math.lsu.edu>
> wrote:
> > Hi,
> >
> > I am trying to use TS to solve a simple transient problem in an
> unstructured finite element f90 code.
> >
> > 1. Section 6.1.1 of the manual refers to a TSSetMatrices function that
> can be used to set the RHS and LHS matrices, but I can;t find it. Is this
> section outdated?
> >
> > Yes, I must have missed this section when updating the documentation.
>
> OK, that makes more sense now.
> >
> > 2. Since we are using unstructured finite elements, the LHS matrix is not
> the identity. As far as I understand, we have two possible choices:
> >   - Use a mass lumping approximation of the variational identity matrix
> (mass matrix), M, and use M^{-1}K for the RHS matrix instead of K.
> >
> > You can also write this as a special case of the choice below where you
> use -ksp_type preonly -pc_type jacobi.
> I am not sure I am following you. Can you elaborate?
>

To solve

M*Xdot = F(X)

IFunction(X,Xdot) = M*Xdot
IJacobian(X,Xdot) = M
RHSFunction(X) = F(X)

Now if you lump M, then you can solve it exactly with one iteration of
Jacobi. My options were just turn off any extra norms that would normally be
computed by an implicit solve. If M is the consistent mass matrix, then you
will need some iterations.


>
> >   - Use an IMEX method where the implicit matrix is the variational
> identity M.
> >  Is this right? What is the recommended way?
> >
> > I would just do this because it's the most flexible. See the user's
> manual section on IMEX methods. If you are interested in adaptive error
> control, then you should also check out -ts_type rosw in petsc-dev. In any
> case, you can write your mass matrix as well as any stiff terms that you
> want to treat implicitly into TSSetIFunction(), provide an (approximate)
> Jacobian with TSSetIJacobian(), and put the rest in TSSetRHSFunction(). You
> can be even more sloppy about it with TSROSW.
> >
> > Look at src/ts/examples/tutorials/ex22.c (has a Fortran twin, ex22f.F) or
> ex25.c in petsc-dev.
> >
>
> How about recasting the problem as a DAE? The documentation seems to imply
> that this is feasible. "For ODE with nontrivial mass matrices such as arise
> in FEM, the implicit/DAE interface significantly reduces overhead to prepare
> the system for algebraic solvers (SNES/KSP) by having the user assemble the
> correctly shifted matrix."
>
> Following ex15, solving
>        u_t = \Delta u
> can be recast as solving
>        F(t,u,\dot u) = 0
> with
>        F(t,u,v) = v-\Delta u.
>
> In this case, the IJacobian would be
>        K+aM,
> where K is the stiffness matrix (K_{i,j} = \int_\Omega \nabla \phi_i \cdot
> \nabla \phi_j\,dx) and M the mass matrix (M_{i,j} = \int_\Omega \phi_i
>  \phi_j\,dx) .
>
> At the continuous level, the IFunction would be v-\Delta u, which cannot be
> evaluated directly in the finite element basis by either solving
>        M F = M\dot u + Ku
> or using mass lumping.
>
> Am I expected to do this in my IFunction or is there a way to pass the mass
> matrix to the TS?
>

You should have some sense for what K is and whether it makes sense to
integrate your system implicitly or not. If you know that you want to use
fully implicit methods, then just dump everything into the IFunction. If
part of your system is non-stiff (and especially if it has less than C^1
regularity, or you need special properties from a certain explicit method),
then you can use the IMEX form G(X,Xdot) = F(X) with

G(X,Xdot) := M*Xdot

You can pass TSComputeIFunctionLinear and/or TSComputeIJacobianConstant to
TSSetIFunction() and TSSetIJacobian respectively if you have a linear
constant mass matrix.


>
> As far as I understand, using IMPEX will lead to the same issue regardless
> of the way I write my problem,  i.e. wether I write G(t,u,v) =  v-\Delta u
> and F(t,u) = 0, or G(t,u,v) =  v and F(t,u) = \Delta u.
>
> As far as I can see, all examples and tests use structured meshes where the
> mass matrix is the identity matrix.
>

Identity as a mass matrix has nothing to do with structured vs.
unstructured, it's a matter of continuous finite elements versus finite
difference/finite volume/certain Petrov-Galerkin methods.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111019/bcc9a8e8/attachment.htm>