Yingjie Wu <[email protected]> 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
