<div class="gmail_quote">On Wed, Oct 19, 2011 at 12:32, Blaise Bourdin <span dir="ltr"><<a href="mailto:bourdin@lsu.edu">bourdin@lsu.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Hi,<br>
<div class="im"><br>
<br>
> On Wed, Oct 19, 2011 at 10:54, Blaise Bourdin <<a href="mailto:bourdin@math.lsu.edu">bourdin@math.lsu.edu</a>> wrote:<br>
> Hi,<br>
><br>
> I am trying to use TS to solve a simple transient problem in an unstructured finite element f90 code.<br>
><br>
> 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?<br>
><br>
> Yes, I must have missed this section when updating the documentation.<br>
<br>
</div>OK, that makes more sense now.<br>
<div class="im">><br>
> 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:<br>
> - 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.<br>
><br>
> You can also write this as a special case of the choice below where you use -ksp_type preonly -pc_type jacobi.<br>
</div>I am not sure I am following you. Can you elaborate?<br></blockquote><div><br></div><div>To solve</div><div><br></div><div>M*Xdot = F(X)</div><div><br></div><div>IFunction(X,Xdot) = M*Xdot</div><div>IJacobian(X,Xdot) = M</div>
<div>RHSFunction(X) = F(X)</div><div><br></div><div>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.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div class="im"><br>
> - Use an IMEX method where the implicit matrix is the variational identity M.<br>
> Is this right? What is the recommended way?<br>
><br>
> 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.<br>
><br>
> Look at src/ts/examples/tutorials/ex22.c (has a Fortran twin, ex22f.F) or ex25.c in petsc-dev.<br>
><br>
<br>
</div>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."<br>
<br>
Following ex15, solving<br>
u_t = \Delta u<br>
can be recast as solving<br>
F(t,u,\dot u) = 0<br>
with<br>
F(t,u,v) = v-\Delta u.<br>
<br>
In this case, the IJacobian would be<br>
K+aM,<br>
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) .<br>
<br>
At the continuous level, the IFunction would be v-\Delta u, which cannot be evaluated directly in the finite element basis by either solving<br>
M F = M\dot u + Ku<br>
or using mass lumping.<br>
<br>
Am I expected to do this in my IFunction or is there a way to pass the mass matrix to the TS?<br></blockquote><div><br></div><div>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</div>
<div><br></div><div>G(X,Xdot) := M*Xdot</div><div><br></div><div>You can pass TSComputeIFunctionLinear and/or TSComputeIJacobianConstant to TSSetIFunction() and TSSetIJacobian respectively if you have a linear constant mass matrix.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<br>
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.<br>
<br>
As far as I can see, all examples and tests use structured meshes where the mass matrix is the identity matrix.<br></blockquote><div><br></div><div>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.</div>
<div><br></div><div><br></div></div>