[petsc-users] Questions about TS

[Blaise Bourdin]

Dear Jed,

I get it now... My confusion was due to being used to derived all scheme on the pde, then discretizing, whereas the documentation assumes that the equation is already discretized. I should have figured it out. 

Is it right to think of the division between ODE, DAE and IMEX in the documentation as Fully Explicit vs. Fully implicit vs. Semi-implicit?

Blaise

On Oct 19, 2011, at 12:55 PM, Jed Brown wrote:

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

Got it.


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

yep. Sorry, I meant finite differences, not structured mesh.


> 
> 

-- 
Department of Mathematics and Center for Computation & Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111019/0b155b4e/attachment-0001.htm>