[petsc-users] Problem in using TS objects

[Yingjie Wu]

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 <jed at jedbrown.org> 于2020年1月27日周一 下午9:54写道：

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200128/398dad3b/attachment-0001.html>