[petsc-users] Questions about TS

[Blaise Bourdin]

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?

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

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.

Thanks for the clarification.
Blaise 
-- 
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