Reconstruct the interfaces states (into primitive variables formulation) by the requested order of accuracy.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| 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. |
pure subroutine reconstruct_interfaces_states(self, primitive, r_primitive)
!--------------------------------------------------------------------------------------------------------------------------------
!< Reconstruct the interfaces states (into primitive variables formulation) by the requested order of accuracy.
!--------------------------------------------------------------------------------------------------------------------------------
class(euler_1D), intent(IN) :: self !< Euler field.
real(R_P), intent(IN) :: primitive(1:self%Np, 1-self%Ng:self%Ni+self%Ng) !< Primitive variables.
real(R_P), intent(INOUT) :: r_primitive(1:self%Np, 1:2, 0:self%Ni+1) !< Reconstructed primitive variables.
real(R_P) :: C(1:2, 1-self%Ng:-1+self%Ng, 1:self%Ns+2) !< Pseudo characteristic variables.
real(R_P) :: CR(1:self%Ns+2, 1:2) !< Pseudo characteristic reconst. vars.
real(R_P) :: Pm(1:self%Np, 1:2) !< Mean of primitive variables.
real(R_P) :: LPm(1:self%Ns+2, 1:self%Ns+2, 1:2) !< Mean left eigenvectors matrix.
real(R_P) :: RPm(1:self%Ns+2, 1:self%Ns+2, 1:2) !< Mean right eigenvectors matrix.
integer(I_P) :: i !< Counter.
integer(I_P) :: j !< Counter.
integer(I_P) :: f !< Counter.
integer(I_P) :: v !< Counter.
!--------------------------------------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------------------------------------
associate(Ni=>self%Ni, Ng=>self%Ng, Ns=>self%Ns, Np=>self%Np, ord=>self%ord, cv0=>self%cv0, cp0=>self%cp0, weno=>self%weno)
select case(ord)
case(1) ! 1st order piecewise constant reconstruction
do i=0, Ni+1
r_primitive(:, 1, i) = primitive(:, i)
r_primitive(:, 2, i) = r_primitive(:, 1, i)
enddo
case(3, 5, 7) ! 3rd, 5th or 7th order WENO reconstruction
do i=0, Ni+1
! trasform primitive variables to pseudo charteristic ones
do f=1, 2
Pm(:,f) = 0.5_R_P * (primitive(:, i+f-2) + primitive(:, i+f-1))
enddo
do f=1, 2
LPm(:, :, f) = eigen_vect_L(Ns=Ns, Np=Np, primitive=Pm(:, f))
RPm(:, :, f) = eigen_vect_R(Ns=Ns, Np=Np, primitive=Pm(:, f))
enddo
do j=i+1-Ng, i-1+Ng
do f=1, 2
do v=1, Ns+2
C(f, j-i, v) = dot_product(LPm(v, 1:Ns+2, f), primitive(1:Ns+2, j))
enddo
enddo
enddo
! compute WENO reconstruction of pseudo charteristic variables
do v=1, Ns+2
call weno%interpolate(S=Ng, &
stencil=C(1:2, 1-Ng:-1+Ng, v), &
location='both', &
interpolation=CR(v, 1:2))
enddo
! trasform back reconstructed pseudo charteristic variables to primitive ones
do f=1, 2
do v=1, Ns+2
r_primitive(v, f, i) = dot_product(RPm(v, 1:Ns+2, f), CR(1:Ns+2, f))
enddo
r_primitive(Ns+3, f, i) = sum(r_primitive(1:Ns, f, i))
r_primitive(Ns+4, f, i) = dot_product(r_primitive(1:Ns, f, i) / r_primitive(Ns+3, f, i), cp0) / &
dot_product(r_primitive(1:Ns, f, i) / r_primitive(Ns+3, f, i), cv0)
enddo
enddo
endselect
endassociate
return
!--------------------------------------------------------------------------------------------------------------------------------
contains
pure function eigen_vect_L(Ns, Np, primitive) result(L)
!-------------------------------------------------------------------------------------------------------------------------------
!< Compute left eigenvectors from primitive variables.
!-------------------------------------------------------------------------------------------------------------------------------
integer(I_P), intent(IN) :: Ns !< Number of initial species.
integer(I_P), intent(IN) :: Np !< Number of primitive variables.
real(R_P), intent(IN) :: primitive(1:Np) !< Primitive variables.
real(R_P) :: L(1:Ns+2,1:Ns+2) !< Left eigenvectors matrix.
real(R_P) :: gp !< g*p.
real(R_P) :: gp_a !< g*p/a.
integer(I_P) :: i !< Counter.
integer(I_P) :: s !< Counter.
!-------------------------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------------------------
gp = primitive(Ns+4) * primitive(Ns+2)
gp_a = gp/a(p=primitive(Ns+2), r=primitive(Ns+3), g=primitive(Ns+4))
L = 0._R_P
L(1, Ns+1) = -gp_a ; L(1, Ns+2) = 1._R_P
do s=2, Ns+1
if (primitive(s-1)>0) L(s, s-1 ) = gp/primitive(s-1) ; L(s, Ns+2) = -1._R_P
enddo
L(Ns+2, Ns+1) = gp_a ; L(Ns+2, Ns+2) = 1._R_P
return
!-------------------------------------------------------------------------------------------------------------------------------
endfunction eigen_vect_L
pure function eigen_vect_R(Ns, Np, primitive) result(R)
!-------------------------------------------------------------------------------------------------------------------------------
!< Compute right eigenvectors from primitive variables.
!-------------------------------------------------------------------------------------------------------------------------------
integer(I_P), intent(IN) :: Ns !< Number of initial species.
integer(I_P), intent(IN) :: Np !< Number of primitive variables.
real(R_P), intent(IN) :: primitive(1:Np) !< Primitive variables.
real(R_P) :: R(1:Ns+2,1:Ns+2) !< Right eigenvectors matrix.
real(R_P) :: gp !< g*p.
real(R_P) :: ss !< Speed of sound, sqrt(g*p/r).
real(R_P) :: gp_inv !< 1/(g*p).
integer(I_P) :: i !< Counter.
integer(I_P) :: s !< Counter.
!-------------------------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------------------------
gp = primitive(Ns+4) * primitive(Ns+2)
ss = a(p=primitive(Ns+2), r=primitive(Ns+3), g=primitive(Ns+4))
gp_inv = 1._R_P/gp
R = 0._R_P
do s=1, Ns
R(s, 1) = 0.5_R_P*primitive(s) * gp_inv ; R(s, s+1) = primitive(s) * gp_inv ; R(s, Ns+2) = R(s, 1)
enddo
R(Ns+1, 1) = -0.5_R_P* ss *gp_inv ; R(Ns+1, Ns+2) = 0.5_R_P* ss * gp_inv
R(Ns+2, 1) = 0.5_R_P ; R(Ns+2, Ns+2) = 0.5_R_P
return
!-------------------------------------------------------------------------------------------------------------------------------
endfunction eigen_vect_R
endsubroutine reconstruct_interfaces_states