euler_1D Derived Type

type, public, extends(integrand) :: euler_1D

type~~euler_1d~~InheritsGraph type~euler_1d euler_1D weno_interpolator_upwind weno_interpolator_upwind weno_interpolator_upwind->type~euler_1d weno integrand integrand integrand->type~euler_1d
Help


Euler 1D PDEs system field.

It is a FOODIE integrand class concrete extension.

1D Euler PDEs system

The 1D Euler PDEs system considered is a non linear, hyperbolic (inviscid) system of conservation laws for compressible gas dynamics, that reads as $$ \begin{matrix} U_t = R(U) \Leftrightarrow U_t = F(U)_x \\ U = \begin{bmatrix} \rho \\ \rho u \\ \rho E \end{bmatrix}\;\;\; F(U) = \begin{bmatrix} \rho u \\ \rho u^2 + p \\ \rho u H \end{bmatrix} \end{matrix} $$ where \(\rho\) is the density, \(u\) is the velocity, \(p\) the pressure, \(E\) the total internal specific energy and \(H\) the total specific enthalpy. The PDEs system must completed with the proper initial and boundary conditions. Moreover, an ideal (thermally and calorically perfect) gas is considered $$ \begin{matrix} R = c_p - c_v \\ \gamma = \frac{c_p}{c_v}\\ e = c_v T \\ h = c_p T \end{matrix} $$ where R is the gas constant, \(c_p\,c_v\) are the specific heats at constant pressure and volume (respectively), e is the internal energy, h is the internal enthalpy and T is the temperature. The following addition equations of state hold: $$ \begin{matrix} T = \frac{p}{\rho R} \\ E = \rho e + \frac{1}{2} \rho u^2 \\ H = \rho h + \frac{1}{2} \rho u^2 \\ a = \sqrt{\frac{\gamma p}{\rho}} \end{matrix} $$

Multi-fluid Euler PDEs system

An extension of the above Euler system is considered allowing the modelling of a multi-fluid mixture of different gas (with different physical characteristics). The well known Standard Thermodynamic Model is used to model the gas mixture replacing the density with the density fraction of each specie composing the mixture. This led to the following system: $$ \begin{matrix} U_t = R(U) \Leftrightarrow U_t = F(U)_x \\ U = \begin{bmatrix} \rho_s \\ \rho u \\ \rho E \end{bmatrix}\;\;\; F(U) = \begin{bmatrix} \rho_s u \\ \rho u^2 + p \\ \rho u H \end{bmatrix}\;\;\; for\; s=1,2,...N_s \\ \rho = \sum_{s=1}^{N_s}\rho_s \\ c_p = \sum_{s=1}^{N_S} \frac{\rho_s}{\rho} c_{p,s} \quad c_v = \sum_{s=1}^{N_S} \frac{\rho_s}{\rho} c_{v,s} \end{matrix} $$ where \(N_s\) is the number of initial species composing the gas mixture.

Numerical grid organization

The finite volume, Godunov's like approach is employed. The conservative variables (and the primitive ones) are co-located at the cell center. The cell and (inter)faces numeration is as follow.

                cell            (inter)faces
                 |                   |
                 v                   v
     |-------|-------|-.....-|-------|-------|-------|-------|-.....-|-------|-------|-------|-.....-|-------|-------|
     | 1-Ng  | 2-Ng  | ..... |  -1   |   0   |   1   |  2    | ..... |  Ni   | Ni+1  | Ni+1  | ..... |Ni+Ng-1| Ni+Ng |
     |-------|-------|-.....-|-------|-------|-------|-------|-.....-|-------|-------|-------|-.....-|-------|-------|
    0-Ng                             -1      0       1       2      Ni-1     Ni                                    Ni+Ng

Where Ni are the finite volumes (cells) used for discretizing the domain and Ng are the ghost cells used for imposing the left and right boundary conditions (for a total of 2Ng cells).

Primitive variables organization

Primitive variables are organized as an array of reals which the first index means:

  • 1 : density of species 1 (r1)
  • 2 : density of species 2 (r2)
  • … :
  • s : density of species s-th (rs)
  • … :
  • Ns : density of species Ns (rNs)
  • Ns+1 : velocity (u)
  • Ns+2 : pressure (p)
  • Ns+3 : density (r=sum(rs))
  • Ns+4 : specific heats ratio (g)

Conservative variables organization

Conservative variables are organized as an array (rank 2) of reals which the first index means:

  • 1 : mass conservation of species 1 (r1)
  • 2 : mass conservation of species 2 (r2)
  • … :
  • s : mass conservation of species s-th (rs)
  • … :
  • Ns : mass conservation of species Ns (rNs)
  • Ns+1 : momentum conservation (r*u)
  • Ns+2 : energy conservation (r*E)

Source Code


Components

TypeVisibility AttributesNameInitial
character(len=:), private, allocatable:: BC_L

Left boundary condition type.

character(len=:), private, allocatable:: BC_R

Right boundary condition type.

real(kind=R_P), private :: Dx =0._R_P

Space step.

integer(kind=I_P), private :: Nc =0

Number of conservative variables, Ns+2.

integer(kind=I_P), private :: Ng =0

Number of ghost cells for boundary conditions handling.

integer(kind=I_P), private :: Ni =0

Space dimension.

integer(kind=I_P), private :: Np =0

Number of primitive variables, Ns+4.

integer(kind=I_P), private :: Ns =0

Number of initial species.

real(kind=R_P), private, allocatable:: U(:,:)

Integrand (state) variables, whole physical domain [1:Nc,1-Ng:Ni+Ng].

real(kind=R_P), private, allocatable:: cp0(:)

Specific heat cp of initial species [1:Ns].

real(kind=R_P), private, allocatable:: cv0(:)

Specific heat cv of initial species [1:Ns].

integer(kind=I_P), private :: ord =0

Space accuracy formal order.

real(kind=R_P), private, allocatable:: previous(:,:,:)

Previous time steps states [1:Nc,1-Ng:Ni+Ng,1:steps].

integer(kind=I_P), private :: steps =0

Number of time steps stored.

type(weno_interpolator_upwind), private :: weno

WENO interpolator.


Type-Bound Procedures

procedure, public, pass(lhs) :: add => add_euler

Euler + Euler operator.

  • private function add_euler(lhs, rhs) result(opr)

    Add two Euler fields.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: lhs

    Left hand side.

    class(integrand), intent(in) :: rhs

    Right hand side.

    Return Value class(integrand), allocatable

    Operator result.

procedure, public, pass(lhs) :: assign_integrand => euler_assign_euler

Euler = Euler.

  • private subroutine euler_assign_euler(lhs, rhs)

    Assign one Euler field to another.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(inout) :: lhs

    Left hand side.

    class(integrand), intent(in) :: rhs

    Right hand side.

procedure, public, pass(lhs) :: assign_real => euler_assign_real

Euler = real.

  • private subroutine euler_assign_real(lhs, rhs)

    Assign one real to an Euler field.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(inout) :: lhs

    Left hand side.

    real(kind=R_P), intent(in) :: rhs

    Right hand side.

procedure, private, pass(self) :: conservative2primitive

Convert conservative variables to primitive ones.

  • private pure function conservative2primitive(self, conservative) result(primitive)

    Convert conservative variables to primitive variables.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    real(kind=R_P), intent(in) :: conservative(:)

    Conservative variables.

    Return Value real(kind=R_P) (1:self%Np)

    Primitive variables.

procedure, public, pass(self) :: destroy

Destroy field.

  • private pure subroutine destroy(self)

    Destroy field.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(inout) :: self

    Euler field.

procedure, public, pass(self) :: dt => compute_dt

Compute the current time step, by means of CFL condition.

  • private pure function compute_dt(self, Nmax, Tmax, t, CFL) result(Dt)

    Compute the current time step by means of CFL condition.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    integer(kind=I_P), intent(in) :: Nmax

    Maximun number of iterates.

    real(kind=R_P), intent(in) :: Tmax

    Maximum time (ignored if Nmax>0).

    real(kind=R_P), intent(in) :: t

    Time.

    real(kind=R_P), intent(in) :: CFL

    CFL value.

    Return Value real(kind=R_P)

    Time step.

procedure, private, pass(self) :: impose_boundary_conditions

Impose boundary conditions.

  • private pure subroutine impose_boundary_conditions(self, primitive)

    Impose boundary conditions.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    real(kind=R_P), intent(inout) :: primitive(1:self%Np,1-self%Ng:self%Ni+self%Ng)

    Primitive variables [1:Np,1-Ng:Ni+Ng].

procedure, public, pass(self) :: init

Init field.

  • private subroutine init(self, Ni, Ns, Dx, BC_L, BC_R, initial_state, cp0, cv0, steps, ord)

    Init field.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(inout) :: self

    Euler field.

    integer(kind=I_P), intent(in) :: Ni

    Space dimension.

    integer(kind=I_P), intent(in) :: Ns

    Number of initial species.

    real(kind=R_P), intent(in) :: Dx

    Space step.

    character(len=*), intent(in) :: BC_L

    Left boundary condition type.

    character(len=*), intent(in) :: BC_R

    Right boundary condition type.

    real(kind=R_P), intent(in) :: initial_state(:,:)

    Initial state of primitive variables.

    real(kind=R_P), intent(in) :: cp0(:)

    Initial specific heat, constant pressure.

    real(kind=R_P), intent(in) :: cv0(:)

    Initial specific heat, constant volume.

    integer(kind=I_P), intent(in), optional :: steps

    Time steps stored.

    integer(kind=I_P), intent(in), optional :: ord

    Space accuracy formal order.

procedure, public, pass(lhs) :: integrand_multiply_integrand => euler_multiply_euler

Euler * Euler operator.

  • private function euler_multiply_euler(lhs, rhs) result(opr)

    Multiply an Euler field by another one.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: lhs

    Left hand side.

    class(integrand), intent(in) :: rhs

    Right hand side.

    Return Value class(integrand), allocatable

    Operator result.

procedure, public, pass(lhs) :: integrand_multiply_real => euler_multiply_real

Euler * real operator.

  • private function euler_multiply_real(lhs, rhs) result(opr)

    Multiply an Euler field by a real scalar.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: lhs

    Left hand side.

    real(kind=R_P), intent(in) :: rhs

    Right hand side.

    Return Value class(integrand), allocatable

    Operator result.

procedure, public, pass(lhs) :: local_error => euler_local_error

Local error.

  • private function euler_local_error(lhs, rhs) result(error)

    Estimate local truncation error between 2 euler approximations.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: lhs

    Left hand side.

    class(integrand), intent(in) :: rhs

    Right hand side.

    Return Value real(kind=R_P)

    Error estimation.

procedure, public, pass(self) :: output

Extract Euler field.

  • private pure function output(self) result(state)

    Output the Euler field state (primitive variables).

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    Return Value real(kind=R_P), dimension(:,:), allocatable

    Euler state vector.

procedure, public, pass(self) :: previous_step

Get a previous time step.

  • private function previous_step(self, n) result(previous)

    Extract previous time solution of Euler field.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    integer(kind=I_P), intent(in) :: n

    Time level.

    Return Value class(integrand), allocatable

    Previous time solution of Euler field.

procedure, private, pass(self) :: primitive2conservative

Convert primitive variables to conservative ones.

  • private pure function primitive2conservative(self, primitive) result(conservative)

    Convert primitive variables to conservative variables.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    real(kind=R_P), intent(in) :: primitive(:)

    Primitive variables.

    Return Value real(kind=R_P) (1:self%Nc)

    Conservative variables.

procedure, public, pass(rhs) :: real_multiply_integrand => real_multiply_euler

Real * Euler operator.

  • private function real_multiply_euler(lhs, rhs) result(opr)

    Multiply a real scalar by an Euler field.

    Arguments

    Type IntentOptional AttributesName
    real(kind=R_P), intent(in) :: lhs

    Left hand side.

    class(euler_1D), intent(in) :: rhs

    Right hand side.

    Return Value class(integrand), allocatable

    Operator result.

procedure, private, pass(self) :: reconstruct_interfaces_states

Reconstruct interfaces states.

  • private pure subroutine reconstruct_interfaces_states(self, primitive, r_primitive)

    Reconstruct the interfaces states (into primitive variables formulation) by the requested order of accuracy.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    real(kind=R_P), intent(in) :: primitive(1:self%Np,1-self%Ng:self%Ni+self%Ng)

    Primitive variables.

    real(kind=R_P), intent(inout) :: r_primitive(1:self%Np,1:2,0:self%Ni+1)

    Reconstructed primitive variables.

procedure, private, pass(self) :: riemann_solver

Solve the Riemann Problem at cell interfaces.

  • private pure subroutine riemann_solver(self, p1, r1, u1, g1, p4, r4, u4, g4, F)

    Solve the Riemann problem between the state $1$ and $4$ using the (local) Lax Friedrichs (Rusanov) solver.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    real(kind=R_P), intent(in) :: p1

    Pressure of state 1.

    real(kind=R_P), intent(in) :: r1

    Density of state 1.

    real(kind=R_P), intent(in) :: u1

    Velocity of state 1.

    real(kind=R_P), intent(in) :: g1

    Specific heats ratio of state 1.

    real(kind=R_P), intent(in) :: p4

    Pressure of state 4.

    real(kind=R_P), intent(in) :: r4

    Density of state 4.

    real(kind=R_P), intent(in) :: u4

    Velocity of state 4.

    real(kind=R_P), intent(in) :: g4

    Specific heats ratio of state 4.

    real(kind=R_P), intent(out) :: F(1:self%Nc)

    Resulting fluxes.

procedure, public, pass(lhs) :: sub => sub_euler

Euler - Euler.

  • private function sub_euler(lhs, rhs) result(opr)

    Subtract two Euler fields.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: lhs

    Left hand side.

    class(integrand), intent(in) :: rhs

    Right hand side.

    Return Value class(integrand), allocatable

    Operator result.

procedure, public, pass(self) :: t => dEuler_dt

Time derivative, residuals function.

  • private function dEuler_dt(self, t) result(dState_dt)

    Time derivative of Euler field, the residuals function.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(in) :: self

    Euler field.

    real(kind=R_P), intent(in), optional :: t

    Time.

    Return Value class(integrand), allocatable

    Euler field time derivative.

procedure, public, pass(self) :: update_previous_steps

Update previous time steps.

  • private subroutine update_previous_steps(self, filter, weights)

    Update previous time steps.

    Arguments

    Type IntentOptional AttributesName
    class(euler_1D), intent(inout) :: self

    Euler field.

    class(integrand), intent(in), optional :: filter

    Filter field displacement.

    real(kind=R_P), intent(in), optional :: weights(:)

    Weights for filtering the steps.

Source Code

type, extends(integrand) :: euler_1D
  !< Euler 1D PDEs system field.
  !<
  !< It is a FOODIE integrand class concrete extension.
  !<
  !<### 1D Euler PDEs system
  !< The 1D Euler PDEs system considered is a non linear, hyperbolic (inviscid) system of conservation laws for compressible gas
  !< dynamics, that reads as
  !<$$
  !<\begin{matrix}
  !<U_t = R(U)  \Leftrightarrow U_t = F(U)_x \\
  !<U = \begin{bmatrix}
  !<\rho \\
  !<\rho u \\
  !<\rho E
  !<\end{bmatrix}\;\;\;
  !<F(U) = \begin{bmatrix}
  !<\rho u \\
  !<\rho u^2 + p \\
  !<\rho u H
  !<\end{bmatrix}
  !<\end{matrix}
  !<$$
  !< where \(\rho\) is the density, \(u\) is the velocity, \(p\) the pressure, \(E\) the total internal specific energy and \(H\)
  !< the total specific enthalpy. The PDEs system must completed with the proper initial and boundary conditions. Moreover, an ideal
  !< (thermally and calorically perfect) gas is considered
  !<$$
  !<\begin{matrix}
  !<R = c_p - c_v \\
  !<\gamma = \frac{c_p}{c_v}\\
  !<e = c_v T \\
  !<h = c_p T
  !<\end{matrix}
  !<$$
  !< where *R* is the gas constant, \(c_p\,c_v\) are the specific heats at constant pressure and volume (respectively), *e* is the
  !< internal energy, *h* is the internal enthalpy and *T* is the temperature. The following addition equations of state hold:
  !<$$
  !<\begin{matrix}
  !<T = \frac{p}{\rho R} \\
  !<E = \rho e + \frac{1}{2} \rho u^2 \\
  !<H = \rho h + \frac{1}{2} \rho u^2 \\
  !<a = \sqrt{\frac{\gamma p}{\rho}}
  !<\end{matrix}
  !<$$
  !<
  !<### Multi-fluid Euler PDEs system
  !< An extension of the above Euler system is considered allowing the modelling of a multi-fluid mixture of different gas (with
  !< different physical characteristics). The well known Standard Thermodynamic Model is used to model the gas mixture replacing the
  !< density with the density fraction of each specie composing the mixture. This led to the following system:
  !<$$
  !<\begin{matrix}
  !<U_t = R(U)  \Leftrightarrow U_t = F(U)_x \\
  !<U = \begin{bmatrix}
  !<\rho_s \\
  !<\rho u \\
  !<\rho E
  !<\end{bmatrix}\;\;\;
  !<F(U) = \begin{bmatrix}
  !<\rho_s u \\
  !<\rho u^2 + p \\
  !<\rho u H
  !<\end{bmatrix}\;\;\; for\; s=1,2,...N_s \\
  !<\rho = \sum_{s=1}^{N_s}\rho_s \\
  !<c_p = \sum_{s=1}^{N_S} \frac{\rho_s}{\rho} c_{p,s} \quad  c_v = \sum_{s=1}^{N_S} \frac{\rho_s}{\rho} c_{v,s}
  !<\end{matrix}
  !<$$
  !< where \(N_s\) is the number of initial species composing the gas mixture.
  !<
  !<#### Numerical grid organization
  !< The finite volume, Godunov's like approach is employed. The conservative variables (and the primitive ones) are co-located at
  !< the cell center. The cell and (inter)faces numeration is as follow.
  !<```
  !<                cell            (inter)faces
  !<                 |                   |
  !<                 v                   v
  !<     |-------|-------|-.....-|-------|-------|-------|-------|-.....-|-------|-------|-------|-.....-|-------|-------|
  !<     | 1-Ng  | 2-Ng  | ..... |  -1   |   0   |   1   |  2    | ..... |  Ni   | Ni+1  | Ni+1  | ..... |Ni+Ng-1| Ni+Ng |
  !<     |-------|-------|-.....-|-------|-------|-------|-------|-.....-|-------|-------|-------|-.....-|-------|-------|
  !<    0-Ng                             -1      0       1       2      Ni-1     Ni                                    Ni+Ng
  !<```
  !< Where *Ni* are the finite volumes (cells) used for discretizing the domain and *Ng* are the ghost cells used for imposing the
  !< left and right boundary conditions (for a total of *2Ng* cells).
  !<
  !<#### Primitive variables organization
  !< Primitive variables are organized as an array of reals which the first index means:
  !<
  !< + 1    : density of species 1    (r1)
  !< + 2    : density of species 2    (r2)
  !< + ...  :
  !< + s    : density of species s-th (rs)
  !< + ...  :
  !< + Ns   : density of species Ns   (rNs)
  !< + Ns+1 : velocity                (u)
  !< + Ns+2 : pressure                (p)
  !< + Ns+3 : density                 (r=sum(rs))
  !< + Ns+4 : specific heats ratio    (g)
  !<
  !<#### Conservative variables organization
  !< Conservative variables are organized as an array (rank 2) of reals which the first index means:
  !<
  !< + 1    : mass conservation of species 1    (r1)
  !< + 2    : mass conservation of species 2    (r2)
  !< + ...  :
  !< + s    : mass conservation of species s-th (rs)
  !< + ...  :
  !< + Ns   : mass conservation of species Ns   (rNs)
  !< + Ns+1 : momentum conservation             (r*u)
  !< + Ns+2 : energy conservation               (r*E)
  private
  integer(I_P)                   :: steps=0         !< Number of time steps stored.
  integer(I_P)                   :: ord=0           !< Space accuracy formal order.
  integer(I_P)                   :: Ni=0            !< Space dimension.
  integer(I_P)                   :: Ng=0            !< Number of ghost cells for boundary conditions handling.
  integer(I_P)                   :: Ns=0            !< Number of initial species.
  integer(I_P)                   :: Nc=0            !< Number of conservative variables, Ns+2.
  integer(I_P)                   :: Np=0            !< Number of primitive variables, Ns+4.
  real(R_P)                      :: Dx=0._R_P       !< Space step.
  type(weno_interpolator_upwind) :: weno            !< WENO interpolator.
  real(R_P),    allocatable      :: U(:,:)          !< Integrand (state) variables, whole physical domain [1:Nc,1-Ng:Ni+Ng].
  real(R_P),    allocatable      :: previous(:,:,:) !< Previous time steps states [1:Nc,1-Ng:Ni+Ng,1:steps].
  real(R_P),    allocatable      :: cp0(:)          !< Specific heat cp of initial species [1:Ns].
  real(R_P),    allocatable      :: cv0(:)          !< Specific heat cv of initial species [1:Ns].
  character(:), allocatable      :: BC_L            !< Left boundary condition type.
  character(:), allocatable      :: BC_R            !< Right boundary condition type.
  contains
    ! auxiliary methods
    procedure, pass(self), public :: init             !< Init field.
    procedure, pass(self), public :: destroy          !< Destroy field.
    procedure, pass(self), public :: output           !< Extract Euler field.
    procedure, pass(self), public :: dt => compute_dt !< Compute the current time step, by means of CFL condition.
    ! ADT integrand deferred methods
    procedure, pass(self), public :: t => dEuler_dt                                       !< Time derivative, residuals function.
    procedure, pass(lhs),  public :: local_error => euler_local_error                     !< Local error.
    procedure, pass(self), public :: update_previous_steps                                !< Update previous time steps.
    procedure, pass(self), public :: previous_step                                        !< Get a previous time step.
    procedure, pass(lhs),  public :: integrand_multiply_integrand => euler_multiply_euler !< Euler * Euler operator.
    procedure, pass(lhs),  public :: integrand_multiply_real => euler_multiply_real       !< Euler * real operator.
    procedure, pass(rhs),  public :: real_multiply_integrand => real_multiply_euler       !< Real * Euler operator.
    procedure, pass(lhs),  public :: add => add_euler                                     !< Euler + Euler operator.
    procedure, pass(lhs),  public :: sub => sub_euler                                     !< Euler - Euler.
    procedure, pass(lhs),  public :: assign_integrand => euler_assign_euler               !< Euler = Euler.
    procedure, pass(lhs),  public :: assign_real => euler_assign_real                     !< Euler = real.
    ! private methods
    procedure, pass(self), private :: primitive2conservative        !< Convert primitive variables to conservative ones.
    procedure, pass(self), private :: conservative2primitive        !< Convert conservative variables to primitive ones.
    procedure, pass(self), private :: impose_boundary_conditions    !< Impose boundary conditions.
    procedure, pass(self), private :: reconstruct_interfaces_states !< Reconstruct interfaces states.
    procedure, pass(self), private :: riemann_solver                !< Solve the Riemann Problem at cell interfaces.
endtype euler_1D