[petsc-users] Problem in using TS objects

[Jed Brown]

Yingjie Wu <yjwu16 at gmail.com> writes:

> Dear PETSc developers
> Hi,
> Recently, I am using PETSc to write a one-dimensional hydrodynamics solver.
> At present, I use SNES object to complete the development of steady-state
> (excluding time term) program, and the result is very good. I want to
> continue with the dynamic solver, but have some problems. As shown in the
> attachment, I solve three conservation equations (partial differential
> equations), and use finite difference to separate one-dimensional meshes.
> The three main variables I solved are velocity V, pressure P, and energy E.
> therefore, I have the following question when writing transient programs
> using TS objects:
>
> 1. Three equations correspond to three variables. I used *TSSetIfunction*
> to set the residual equation. In theory, I can use *Vec u_t* (time
> derivative of state vector) to set the time term. But the time term in *Eq1*
> is the time partial derivative of density *\ rho*. The density *\rho* is a
> function of energy E and pressure P. How to set this time term which is not
> a main variable, but an intermediate variable (*\rho (E, P)*))  which is
> composed of two main variables?

The standard answer to this problem is to use the chain rule to expand
the time derivative in terms of the state variables.

  \rho_P P_t + \rho_E E_t + (\rho V)_z = 0.

The trouble with this is that you will generally lose discrete
conservation unless the function \rho(P,E) is sufficiently simple, e.g.,
linear.  Conservation error will converge to zero, but only at the order
of time discretization, versus being exact (up to machine epsilon or
algebraic solver tolerances).

In the conservation law community, the standard answer to this is to
insist on using conservative variables for the state.  This requires
that equations of state be invertible, e.g.,

  P(\rho, \rho e) = ...,

perhaps by way of an implicit solve (implicitly-defined constitutive
models are used in a number of domains).  Sometimes this has negative
effects on the conditioning of the algebraic systems; that can generally
be repaired by preconditioning.

A method that can be used with any time integrator and preconditioner is
to expand the equations into a DAE, writing

  (conservative) - eq_of_state(primitive) = 0
  (conservative)_t - f(primitive) = 0

This increases the number of degrees of freedom in your system and may
add solver complexity (though it can typically be done at low overhead).

As chance would have it, I'll soon be submitting a merge request that
will let you set a callback to define "transient" variables (your
conservative variables in this case) as a function of the state
variables.  Methods like BDF will then become conservative with a simple
implementation (your IFunction will get time derivative of conserved
variables as an input), though you'll still be expected to handle the
chain rule correctly in IJacobian.  This transformation won't work for
explicit methods or Runge-Kutta methods.

> (In the equations, the direction of the mesh is z. g represents the
> acceleration of gravity. f_vis for resistance pressure drop, Q for heat
> source)
>
> I need some advices, if there are examples better.
>
> Thanks,
> Yingjie