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>
