Hi Jed, Thank you very much for your detailed answer. It helps me a lot. I'm going to use the "chain rule" method to solve my problem first, and may continue to try different ways later. For your reply, I have the following questions: 1. The method of transforming equations into DAE, for my problem, does it mean taking \ rho as a variable and adding its constitutive equation as a residual equation? f(\rho) = \rho - subroutine_compute_rho(P,E)
2. I hope to use finite difference to construct Jacobian matrix first(this means that I may not provide the IJacobian function). Is the previous switch "-snes_fd" still available for transient calculation? 3. For the new merge request you have mentioned, where should I get the information so as to get its new developments? Thanks again for your help. Yingjie Jed Brown <[email protected]> 于2020年1月27日周一 下午9:54写道： > 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 >
