<div dir="ltr">Hi Jed,<div>Thank you very much for your detailed answer. It helps me a lot. </div><div>I'm going to use the "chain rule" method to solve my problem first, and may continue to try different ways later. </div><div>For your reply, I have the following questions: </div><div> 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?</div><div> f(\rho) = \rho - subroutine_compute_rho(P,E) </div><div><br></div><div> 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? <br></div><div><br></div><div> 3. For the new merge request you have mentioned, where should I get the information so as to get its new developments? </div><div><br></div><div>Thanks again for your help.<br></div><div>Yingjie</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>> 于2020年1月27日周一 下午9:54写道：<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Yingjie Wu <<a href="mailto:yjwu16@gmail.com" target="_blank">yjwu16@gmail.com</a>> writes:<br>
<br>
> Dear PETSc developers<br>
> Hi,<br>
> Recently, I am using PETSc to write a one-dimensional hydrodynamics solver.<br>
> At present, I use SNES object to complete the development of steady-state<br>
> (excluding time term) program, and the result is very good. I want to<br>
> continue with the dynamic solver, but have some problems. As shown in the<br>
> attachment, I solve three conservation equations (partial differential<br>
> equations), and use finite difference to separate one-dimensional meshes.<br>
> The three main variables I solved are velocity V, pressure P, and energy E.<br>
> therefore, I have the following question when writing transient programs<br>
> using TS objects:<br>
><br>
> 1. Three equations correspond to three variables. I used *TSSetIfunction*<br>
> to set the residual equation. In theory, I can use *Vec u_t* (time<br>
> derivative of state vector) to set the time term. But the time term in *Eq1*<br>
> is the time partial derivative of density *\ rho*. The density *\rho* is a<br>
> function of energy E and pressure P. How to set this time term which is not<br>
> a main variable, but an intermediate variable (*\rho (E, P)*)) which is<br>
> composed of two main variables?<br>
<br>
The standard answer to this problem is to use the chain rule to expand<br>
the time derivative in terms of the state variables.<br>
<br>
\rho_P P_t + \rho_E E_t + (\rho V)_z = 0.<br>
<br>
The trouble with this is that you will generally lose discrete<br>
conservation unless the function \rho(P,E) is sufficiently simple, e.g.,<br>
linear. Conservation error will converge to zero, but only at the order<br>
of time discretization, versus being exact (up to machine epsilon or<br>
algebraic solver tolerances).<br>
<br>
In the conservation law community, the standard answer to this is to<br>
insist on using conservative variables for the state. This requires<br>
that equations of state be invertible, e.g.,<br>
<br>
P(\rho, \rho e) = ...,<br>
<br>
perhaps by way of an implicit solve (implicitly-defined constitutive<br>
models are used in a number of domains). Sometimes this has negative<br>
effects on the conditioning of the algebraic systems; that can generally<br>
be repaired by preconditioning.<br>
<br>
A method that can be used with any time integrator and preconditioner is<br>
to expand the equations into a DAE, writing<br>
<br>
(conservative) - eq_of_state(primitive) = 0<br>
(conservative)_t - f(primitive) = 0<br>
<br>
This increases the number of degrees of freedom in your system and may<br>
add solver complexity (though it can typically be done at low overhead).<br>
<br>
As chance would have it, I'll soon be submitting a merge request that<br>
will let you set a callback to define "transient" variables (your<br>
conservative variables in this case) as a function of the state<br>
variables. Methods like BDF will then become conservative with a simple<br>
implementation (your IFunction will get time derivative of conserved<br>
variables as an input), though you'll still be expected to handle the<br>
chain rule correctly in IJacobian. This transformation won't work for<br>
explicit methods or Runge-Kutta methods.<br>
<br>
> (In the equations, the direction of the mesh is z. g represents the<br>
> acceleration of gravity. f_vis for resistance pressure drop, Q for heat<br>
> source)<br>
><br>
> I need some advices, if there are examples better.<br>
><br>
> Thanks,<br>
> Yingjie<br>
</blockquote></div>