diff --git a/Solver/Makefile b/Solver/Makefile index 39691e9..4637cc1 100755 --- a/Solver/Makefile +++ b/Solver/Makefile @@ -48,10 +48,10 @@ ProblemFile: SMConstants Setup InitialConditions: SMConstants Physics Setup: SMConstants Paramfile DGWeakIntegrals: SMConstants Physics QuadElement MatrixOperations -DGSpatialDiscretizationMethods: SMConstants Physics QuadMesh DGViscousMethods MatrixOperations DGInviscidMethods DGWeakIntegrals +DGSpatialDiscretizationMethods: SMConstants Physics QuadMesh DGViscousMethods MatrixOperations DGInviscidMethods DGWeakIntegrals DGArtificialDissipation +DGArtificialDissipation: SMConstants QuadMesh QuadElement Setup DGViscousMethods: SMConstants Physics Setup QuadElement MatrixOperations NodesAndWeightsClass QuadMesh Headers DGWeakIntegrals DGInviscidMethods: SMConstants Physics Setup QuadElement MatrixOperations NodesAndWeightsClass -DGInviscid_StandardDG: DGInviscidMethods DGViscous_BR1: DGViscousMethods DGViscous_IP: DGViscousMethods DGBoundaryConditions: SMConstants Setup Physics QuadElement Paramfile @@ -132,7 +132,7 @@ INCLUDE= -I$(INC_DIR) -I$(NETCDF_INC) $(F90_INCLUDE) all: def_navier_stokes_macro $(EXECUTABLE) runcyl: - (cd ../TestCases/Benchmarktests/FlowOverCircle && ./$(EXECUTABLE) ./SETUP/Cylinder.HiOCase -check) + (cd ../TestCases/Benchmarktests/FlowOverCircle && ../../../Solver/bin/$(EXECUTABLE) ./SETUP/Cylinder.HiOCase -check) debugcyl: (cd ../TestCases/Benchmarktests/FlowOverCircle && gdb ./$(EXECUTABLE)) diff --git a/Solver/src/Approximation/DGArtificialDissipation.f90 b/Solver/src/Approximation/DGArtificialDissipation.f90 new file mode 100644 index 0000000..5038bdf --- /dev/null +++ b/Solver/src/Approximation/DGArtificialDissipation.f90 @@ -0,0 +1,340 @@ +#ifdef NAVIER_STOKES + +module DGArtificialDissipation + use SMConstants + use QuadMeshClass + use QuadElementClass + use Physics + implicit none + +#include "Defines.h" + + + private + public ArtificialDissipation_t , ArtificialDissipation_Initialization +! +!//////////////////////////////////////////////////////////////////////////// +! +! ARTIFICIAL DISSIPATION CLASSES +! ------------------------------ +!//////////////////////////////////////////////////////////////////////////// +! + type ArtificialDissipation_t + real(kind=RP) :: Ceps + procedure(ElementViscosityFCN), private, nopass, pointer :: ElementViscosity => NULL() + procedure(EdgeViscosityFCN), private, nopass, pointer :: EdgeViscosity => NULL() + contains + procedure :: ComputeVolumeFluxes => ArtificialDissipation_ComputeVolumeFluxes + procedure :: ComputeFaceFluxes => ArtificialDissipation_ComputeFaceFluxes + procedure :: ComputeElementViscosity => ArtificialDissipation_ComputeElementViscosity + procedure :: ComputeEdgeViscosity => ArtificialDissipation_ComputeEdgeViscosity + end type ArtificialDissipation_t + + type, extends(ArtificialDissipation_t) :: LaplaceDissipation_t + contains + procedure :: ComputeVolumeFluxes => LaplaceDissipation_ComputeVolumeFluxes + procedure :: ComputeFaceFluxes => LaplaceDissipation_ComputeFaceFluxes + end type LaplaceDissipation_t + + type, extends(ArtificialDissipation_t) :: PhysicalDissipation_t + contains + procedure :: ComputeElementViscosity => PhysicalDissipation_ComputeElementViscosity + procedure :: ComputeEdgeViscosity => PhysicalDissipation_ComputeEdgeViscosity + end type PhysicalDissipation_t + + abstract interface + pure function ElementViscosityFCN( self , e ) result ( mu ) + use SMConstants + use QuadElementClass + import ArtificialDissipation_t + implicit none + class(ArtificialDissipation_t), intent(in) :: self + class(QuadElement_t) , intent(in) :: e + real(kind=RP) :: mu + end function ElementViscosityFCN + + pure function EdgeViscosityFCN( self , edge ) result ( mu ) + use SMConstants + use QuadElementClass + import ArtificialDissipation_t + implicit none + class(ArtificialDissipation_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge + real(kind=RP) :: mu(NDIM) + end function EdgeViscosityFCN + end interface +! + +! ======== + contains +! ======== +! +! +!////////////////////////////////////////////////////////////////////////////////////// +! +! INITIALIZATION +! -------------- +!////////////////////////////////////////////////////////////////////////////////////// +! + function ArtificialDissipation_Initialization() result( ArtificialDissipation ) + use Setup_class + implicit none + class(ArtificialDissipation_t), pointer :: ArtificialDissipation + + + if ( Setup % artificialDissipation ) then + + if ( trim( Setup % artificialDissipationType) .eq. "Laplacian" ) then + allocate ( LaplaceDissipation_t :: ArtificialDissipation ) + + elseif ( trim( Setup % artificialDissipationType) .eq. "Physical" ) then + allocate ( PhysicalDissipation_t :: ArtificialDissipation ) + + else + allocate ( ArtificialDissipation_t :: ArtificialDissipation ) + + end if + + if ( trim(Setup % artificialDissipationIndicator) .eq. "Residuals-based" ) then +! ArtificialDissipation % ElementViscosity => ResidualsBasedElementViscosity +! ArtificialDissipation % EdgeViscosity => ResidualsBasedEdgeViscosity + + elseif ( trim(Setup % artificialDissipationIndicator) .eq. "Jumps-based" ) then + ArtificialDissipation % ElementViscosity => JumpsBasedElementViscosity + ArtificialDissipation % EdgeViscosity => JumpsBasedEdgeViscosity + + end if + + else + allocate ( ArtificialDissipation_t :: ArtificialDissipation ) + + end if + + ArtificialDissipation % Ceps = Setup % artificialDissipationIntensity + + + end function ArtificialDissipation_Initialization +! +!////////////////////////////////////////////////////////////////////////////////////// +! +! ARTIFICIAL DISSIPATION PROCEDURES +! --------------------------------- +!////////////////////////////////////////////////////////////////////////////////////// +! +!TODO: pure + function ArtificialDissipation_ComputeVolumeFluxes( self , e ) result ( F ) + implicit none + class(ArtificialDissipation_t), intent (in) :: self + class(QuadElement_t) , intent (in) :: e + real(kind=RP) :: F( 0 : e % spA % N , 0 : e % spA % N , 1 : NCONS , 1:NDIM) +! +! --------------------------- +! The base class does nothing. +! --------------------------- +! + F = 0.0_RP + + end function ArtificialDissipation_ComputeVolumeFluxes + + function ArtificialDissipation_ComputeFaceFluxes( self , edge , QL , QR , dQL , dQR , normal ) result ( F ) + implicit none + class(ArtificialDissipation_t), intent (in) :: self + class(Edge_t) , intent (in) :: edge + real(kind=RP) , intent (in) :: QL(0 : edge % spA % N , 1:NCONS) + real(kind=RP) , intent (in) :: QR(0 : edge % spA % N , 1:NCONS) + real(kind=RP) , intent (in) :: dQL(0 : edge % spA % N , 1:NDIM , 1:NCONS) + real(kind=RP) , intent (in) :: dQR(0 : edge % spA % N , 1:NDIM , 1:NCONS) + real(kind=RP) , intent (in) :: normal(1:NDIM , 0 : edge % spA % N ) + real(kind=RP) :: F( 0 : edge % spA % N , 1 : NCONS) +! +! --------------------------- +! The base class does nothing. +! --------------------------- +! + F = 0.0_RP + + end function ArtificialDissipation_ComputeFaceFluxes + + function ArtificialDissipation_ComputeElementViscosity( self , e ) result ( mu ) + implicit none + class(ArtificialDissipation_t), intent (in) :: self + class(QuadElement_t) , intent (in) :: e + real(kind=RP) :: mu +! +! --------------------------- +! The base class does nothing. +! --------------------------- +! + mu = 0.0_RP + + end function ArtificialDissipation_ComputeElementViscosity + + function ArtificialDissipation_ComputeEdgeViscosity( self , ed ) result ( mu ) + implicit none + class(ArtificialDissipation_t), intent (in) :: self + class(Edge_t) , intent (in) :: ed + real(kind=RP) :: mu +! +! --------------------------- +! The base class does nothing. +! --------------------------- +! + mu = 0.0_RP + + end function ArtificialDissipation_ComputeEdgeViscosity + + function LaplaceDissipation_ComputeVolumeFluxes ( self , e ) result ( F ) + implicit none + class(LaplaceDissipation_t) , intent(in) :: self + class(QuadElement_t) , intent(in) :: e + real(kind=RP) :: F( 0 : e % spA % N , 0 : e % spA % N , 1 : NCONS , 1:NDIM) +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: dQ( 0 : e % spA % N , 0 : e % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP) :: mu + + dQ = e % ComputeInteriorGradient() + mu = self % ElementViscosity(self,e) + + F(:,:,:,IX) = mu * dQ(:,:,IX,:) + F(:,:,:,IY) = mu * dQ(:,:,IY,:) + + end function LaplaceDissipation_ComputeVolumeFluxes + + function LaplaceDissipation_ComputeFaceFluxes( self, edge , QL , QR , dQL , dQR , normal) result ( F ) + implicit none + class(LaplaceDissipation_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge + real(kind=RP) , intent(in) :: QL(0 : edge % spA % N , 1:NCONS) + real(kind=RP) , intent(in) :: QR(0 : edge % spA % N , 1:NCONS) + real(kind=RP) , intent(in) :: dQL(0 : edge % spA % N , 1:NDIM , 1:NCONS) + real(kind=RP) , intent(in) :: dQR(0 : edge % spA % N , 1:NDIM , 1:NCONS) + real(kind=RP) , intent(in) :: normal(1:NDIM , 0 : edge % spA % N) + real(kind=RP) :: F( 0 : edge % spA % N , 1 : NCONS ) +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: mu + + mu = 0.5_RP * sum(self % EdgeViscosity(self,edge)) + F = - mu * edge % Area * ( QL - QR ) + + end function LaplaceDissipation_ComputeFaceFluxes + + pure function PhysicalDissipation_ComputeElementViscosity( self , e ) result ( mu ) + implicit none + class(PhysicalDissipation_t) , intent(in) :: self + class(QuadElement_t) , intent(in) :: e + real(kind=RP) :: mu + + mu = self % ElementViscosity(self,e) + + end function PhysicalDissipation_ComputeElementViscosity + + pure function PhysicalDissipation_ComputeEdgeViscosity( self , ed ) result ( mu ) + implicit none + class(PhysicalDissipation_t), intent(in) :: self + class(Edge_t) , intent(in) :: ed + real(kind=RP) :: mu + + mu = maxval(self % EdgeViscosity(self,ed)) + + end function PhysicalDissipation_ComputeEdgeViscosity +! +!///////////////////////////////////////////////////////////////////////////////////////////////// +! +! ARTIFICIAL VISCOSITY ESTIMATORS +! ------------------------------- +!///////////////////////////////////////////////////////////////////////////////////////////////// +! + pure function JumpsBasedElementViscosity(self , e ) result ( mu ) + implicit none + class(ArtificialDissipation_t), intent(in) :: self + class(QuadElement_t) , intent(in) :: e + real(kind=RP) :: mu +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: gk , gkRHO , gkRHOE + real(kind=RP) :: capitalGk + integer :: ed + real(kind=RP) :: xiL + real(kind=RP) :: xiR + real(kind=RP) :: maxVisc + real(kind=RP) :: s0 + real(kind=RP), parameter :: k = 5.0_RP +! +! The jumps indicator is computed squaring all interface jumps +! ------------------------------------------------------------ + gkRHO = e % edges(EBOTTOM) % f % computeJumps(IRHO) & + + e % edges(ERIGHT ) % f % computeJumps(IRHO) & + + e % edges(ETOP ) % f % computeJumps(IRHO) & + + e % edges(ELEFT ) % f % computeJumps(IRHO) + + gkRHOE = e % edges(EBOTTOM) % f % computeJumps(IRHOE) & + + e % edges(ERIGHT ) % f % computeJumps(IRHOE) & + + e % edges(ETOP ) % f % computeJumps(IRHOE) & + + e % edges(ELEFT ) % f % computeJumps(IRHOE) + + gk = max(gkRHO,gkRHOE) + + gk = gk / ( e % edges(EBOTTOM) % f % Area + e % edges(ERIGHT) % f % Area & + +e % edges(ETOP) % f % Area + e % edges(ELEFT) % f % Area ) + + if ( gk .gt. 1.0e-8 ) then + gk = log10(gk) + else + gk = -HUGE(1.0_RP) + end if + +! +! Compute the smooth discrete jumps indicator +! ------------------------------------------- + maxVisc = sqrt(e % Volume) / ( e % spA % N + 1 ) + s0 = log10( 1.0_RP / ( e % spA % N + 1) ** 4.0_RP ) + 1 + + if ( gk .lt. (s0-k) ) then + capitalGk = 0.0_RP + + elseif ( gk .lt. (s0+k) ) then + capitalGk = 0.5_RP * sin( 0.5_RP * PI * ( gk - s0 )/(k) ) + 0.5_RP + + else + capitalGk = 1.0_RP + + end if + + mu = self % Ceps * maxVisc * capitalGk + + end function JumpsBasedElementViscosity + + pure function JumpsBasedEdgeViscosity( self , ed ) result ( mu ) + implicit none + class(ArtificialDissipation_t), intent(in) :: self + class(Edge_t) , intent(in) :: ed + real(kind=RP) :: mu(NDIM) + + if ( size( ed % storage ) .eq. 1 ) then + mu = self % ElementViscosity( self , ed % quads(1) % e ) + + else + mu(LEFT) = self % ElementViscosity( self , ed % quads(LEFT) % e ) + mu(RIGHT) = self % ElementViscosity( self , ed % quads(RIGHT) % e ) + + end if + + end function JumpsBasedEdgeViscosity + +! + + +end module DGArtificialDissipation + +#endif diff --git a/Solver/src/Approximation/DGInviscidMethods.f90 b/Solver/src/Approximation/DGInviscidMethods.f90 index cb92d72..789be81 100755 --- a/Solver/src/Approximation/DGInviscidMethods.f90 +++ b/Solver/src/Approximation/DGInviscidMethods.f90 @@ -23,10 +23,6 @@ module DGInviscidMethods integer :: formulation procedure(RiemannSolverFunction), pointer, nopass :: RiemannSolver => NULL() contains - generic, public :: ComputeRiemannFluxes => ComputeRiemannFluxes_Interior , ComputeRiemannFluxes_StraightBdry , ComputeRiemannFluxes_CurvedBdry - procedure, private :: ComputeRiemannFluxes_Interior => StdDG_ComputeRiemannFluxes_Interior - procedure, private :: ComputeRiemannFluxes_StraightBdry => StdDG_ComputeRiemannFluxes_StraightBdry - procedure, private :: ComputeRiemannFluxes_CurvedBdry => StdDG_ComputeRiemannFluxes_CurvedBdry procedure :: ComputeInnerFluxes => StdDG_ComputeInnerFluxes procedure, non_overridable :: Describe => InviscidMethod_describe end type InviscidMethod_t @@ -47,37 +43,7 @@ module DGInviscidMethods type, extends(InviscidMethod_t) :: SplitDG_t real(kind=RP) :: alpha end type SplitDG_t - - interface - module subroutine StdDG_ComputeRiemannFluxes_Interior( self , ed , FL , FR ) - use MatrixOperations - use QuadElementClass - implicit none - class(InviscidMethod_t) :: self - type(Edge_t) :: ed - real ( kind=RP ) :: FL ( 0 : ed % storage ( LEFT ) % spA % N , 1 : NCONS ) - real ( kind=RP ) :: FR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS ) - end subroutine StdDG_ComputeRiemannFluxes_Interior - - module subroutine StdDG_ComputeRiemannFluxes_StraightBdry( self , ed , F ) - use MatrixOperations - use QuadElementClass - implicit none - class(InviscidMethod_t) :: self - type(StraightBdryEdge_t) :: ed - real ( kind=RP ) :: F ( 0 : ed % spA % N , 1 : NCONS ) - end subroutine StdDG_ComputeRiemannFluxes_StraightBdry - - module subroutine StdDG_ComputeRiemannFluxes_CurvedBdry( self , ed , F ) - use MatrixOperations - use QuadElementClass - implicit none - class(InviscidMethod_t) :: self - type(CurvedBdryEdge_t) :: ed - real ( kind=RP ) :: F ( 0 : ed % spA % N , 1 : NCONS ) - end subroutine StdDG_ComputeRiemannFluxes_CurvedBdry - end interface - +! ! ! ======== contains @@ -196,10 +162,7 @@ end function InviscidMethod_Initialization pure function StdDG_ComputeInnerFluxes( self , e ) result (F) ! ! ********************************************************************** -! This subroutine computes the contravariant fluxes of the element -! The fluxes read: -! F <- F * Ja(1,1) + G * Ja(2,1) -! G <- F * Ja(1,2) + G * Ja(2,2) +! This subroutine computes the cartesian fluxes of the element ! ********************************************************************** ! use Physics @@ -207,25 +170,8 @@ pure function StdDG_ComputeInnerFluxes( self , e ) result (F) class(InviscidMethod_t), intent (in) :: self class(QuadElement_t), intent (in) :: e real(kind=RP) :: F(0:e % spA % N , 0:e % spA % N , 1:NCONS , 1:NDIM) -! ------------------------------------------------------------- - real(kind=RP) :: F_cartesian(0:e % spA % N,0:e % spA % N,1:NCONS,1:NDIM) - integer :: eq - integer :: N - - N = e % spA % N - F_cartesian = InviscidFlux( e % spA % N , e % Q ) - - do eq = 1 , NCONS -! -! F flux (contravariant) -! ---------------------- - F(0:N,0:N,eq,IX) = F_cartesian(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,1) + F_cartesian(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,1) -! -! G flux (contravariant) -! ---------------------- - F(0:N,0:N,eq,IY) = F_cartesian(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,2) + F_cartesian(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,2) - end do + F = InviscidFlux( e % spA % N , e % Q ) end function StdDG_ComputeInnerFluxes ! diff --git a/Solver/src/Approximation/DGInviscid_StandardDG.f90 b/Solver/src/Approximation/DGInviscid_StandardDG.f90 deleted file mode 100644 index a7f9677..0000000 --- a/Solver/src/Approximation/DGInviscid_StandardDG.f90 +++ /dev/null @@ -1,278 +0,0 @@ -submodule (DGInviscidMethods) DGInviscid_StandardDG - use Physics - -#include "Defines.h" - - contains - module subroutine StdDG_ComputeRiemannFluxes_Interior( self , ed , FL , FR ) - use MatrixOperations - use QuadElementClass - implicit none - class(InviscidMethod_t) :: self - type(Edge_t) :: ed - real ( kind=RP ) :: FL ( 0 : ed % storage ( LEFT ) % spA % N , 1 : NCONS ) - real ( kind=RP ) :: FR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS ) -! --------------------------------------------------------------------------------------------- - real ( kind=RP ), target :: QL ( 0 : ed % spA % N , 1 : NCONS ) - real ( kind=RP ), target :: QR ( 0 : ed % spA % N , 1 : NCONS ) - real ( kind=RP ) :: QL1D(1:NCONS) , QR1D(1:NCONS) - integer :: iXi - - if ( ed % transform(LEFT) .and. ed % inverse ) then -! -! Transform the LEFT edge -! ----------------------- - call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) -! -! Invert the RIGHT edge -! --------------------- - QR = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , 1:NCONS ) -! -! Compute the Riemann solver -! -------------------------- - do iXi = 0 , ed % spA % N - QL1D = QL(iXi,:) - QR1D = QR(iXi,:) - FR(iXi,1:NCONS) = self % RiemannSolver(QL1D , QR1D , ed % n(IX:IY,0) ) * ed % dS(0) - end do -! -! Transform the LEFT -! ----------------------- - call Mat_x_Mat( ed % T_backward , FR , FL ) -! -! Invert the RIGHT -! --------------------- - FR(0:ed % spA % N , 1:NCONS) = FR(ed % spA % N : 0 : -1 , 1:NCONS) - - elseif ( ed % transform(LEFT) ) then -! -! Transform the LEFT -! ----------------------- - call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) -! -! Get the RIGHT state -! ------------------------ - QR = ed % storage(RIGHT) % Q -! -! Compute the Riemann solver -! -------------------------- - do iXi = 0 , ed % spA % N - QL1D = QL(iXi,:) - QR1D = QR(iXi,:) - FR(iXi,1:NCONS) = self % RiemannSolver(QL1D , QR1D , ed % n(IX:IY,0) ) * ed % dS(0) - end do -! -! Transform the LEFT -! ----------------------- - call Mat_x_Mat( ed % T_backward , FR , FL ) - - - elseif ( ed % transform(RIGHT) .and. ed % inverse ) then -! -! Get the LEFT -! ----------------- - QL = ed % storage(LEFT) % Q -! -! Transform the RIGHT -! ------------------------ - call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) -! -! Invert the RIGHT -! --------------------- - QR(0:ed % spA % N,1:NCONS) = QR(ed % spA % N : 0 : -1 , 1:NCONS) -! -! Compute the Riemann solver -! -------------------------- - do iXi = 0 , ed % spA % N - QL1D = QL(iXi,:) - QR1D = QR(iXi,:) - FL(iXi,1:NCONS) = self % RiemannSolver(QL1D , QR1D , ed % n(IX:IY,0) ) * ed % dS(0) - end do -! -! Undo the transformation for the RIGHT -! ------------------------------------------ - call Mat_x_Mat( ed % T_backward , FL , FR ) -! -! Invert the -! --------------- - FR(0:ed % NLow,1:NCONS) = FR(ed % NLow:0:-1 , 1:NCONS) - - elseif ( ed % transform(RIGHT) ) then -! -! Get the LEFT -! ----------------- - QL = ed % storage(LEFT) % Q -! -! Transform the RIGHT -! ------------------------ - call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) -! -! Compute the Riemann solver -! -------------------------- - do iXi = 0 , ed % spA % N - QL1D = QL(iXi,:) - QR1D = QR(iXi,:) - FL(iXi,1:NCONS) = self % RiemannSolver(QL1D , QR1D , ed % n(IX:IY,0) ) * ed % dS(0) - end do -! -! -! Undo the transformation for the RIGHT -! ------------------------------------------ - call Mat_x_Mat( ed % T_backward , FL , FR ) - - elseif ( ed % inverse ) then -! -! Get the LEFT -! ----------------- - QL = ed % storage(LEFT) % Q -! -! Invert the RIGHT -! --------------------- - QR(0:ed % spA % N , 1:NCONS ) = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , 1:NCONS) -! -! Compute the Riemann solver -! -------------------------- - do iXi = 0 , ed % spA % N - QL1D = QL(iXi,:) - QR1D = QR(iXi,:) - FL(iXi,1:NCONS) = self % RiemannSolver(QL1D , QR1D , ed % n(IX:IY,0) ) * ed % dS(0) - end do -! -! Invert the RIGHT -! --------------------- - FR( 0 : ed % spA % N , 1:NCONS ) = FL ( ed % spA % N : 0 : -1 , 1:NCONS ) - - else -! -! Get the LEFT -! ----------------- - QL = ed % storage(LEFT) % Q -! -! Get the RIGHT -! ------------------ - QR = ed % storage(RIGHT) % Q -! -! Compute the Riemann solver -! -------------------------- - do iXi = 0 , ed % spA % N - QL1D = QL(iXi,:) - QR1D = QR(iXi,:) - FL(iXi,1:NCONS) = self % RiemannSolver(QL1D , QR1D , ed % n(IX:IY,0) ) * ed % dS(0) - FR(iXi,1:NCONS) = FL(iXi,1:NCONS) - end do - - end if -! -! Change the sign of FR so that it points towards the element outside -! ------------------------------------------------------------------- - FR = -1.0_RP * FR - - end subroutine StdDG_ComputeRiemannFluxes_Interior - - module subroutine StdDG_ComputeRiemannFluxes_StraightBdry( self , ed , F ) - use QuadElementClass - implicit none - class(InviscidMethod_t) :: self - type(StraightBdryEdge_t) :: ed - real(kind=RP) :: F( 0 : ed % spA % N , 1 : NCONS ) - real(kind=RP) :: QL(1:NCONS) , QR(1:NCONS) - integer, pointer :: N - integer :: iXi - procedure(RiemannSolverFunction), pointer :: RiemannSolver - - N => ed % spA % N -! -! Select boundary condition type -! ------------------------------ - select case ( ed % BCWeakType ) - - case ( WEAK_PRESCRIBED ) -! -! Prescribed boundary conditions: just grab the value -! --------------------------------------------------- - F = ed % FB - - case ( WEAK_RIEMANN ) -! -! Weak boundary conditions -! ------------------------ - if ( associated ( ed % RiemannSolver ) ) then - RiemannSolver => ed % RiemannSolver - - else - RiemannSolver => self % RiemannSolver - - end if -! -! Select LEFT and RIGHT states -! ---------------------------- - if ( .not. ed % inverse ) then ! This just occurs in periodic BCs - do iXi = 0 , N - QL = ed % storage(1) % Q(iXi , 1:NCONS) - QR = ed % uB(iXi , 1:NCONS) - - F(iXi , 1:NCONS) = RiemannSolver( QL , QR , ed % n(IX:IY,0) ) * ed % dS(0) - end do - - else - do iXi = 0 , N - QL = ed % storage(1) % Q(iXi , 1:NCONS) - QR = ed % uB(N - iXi , 1:NCONS) - - F(iXi , 1:NCONS) = RiemannSolver( QL , QR , ed % n(IX:IY,0) ) * ed % dS(0) - end do - - end if - end select - - end subroutine StdDG_ComputeRiemannFluxes_StraightBdry - - module subroutine StdDG_ComputeRiemannFluxes_CurvedBdry( self , ed , F) - use QuadElementClass - implicit none - class(InviscidMethod_t) :: self - type(CurvedBdryEdge_t) :: ed - real(kind=RP) :: F( 0 : ed % spA % N , 1 : NCONS ) - real(kind=RP) :: QL(1:NCONS) , QR(1:NCONS) - integer, pointer :: N - integer :: iXi - procedure(RiemannSolverFunction), pointer :: RiemannSolver - - N => ed % spA % N -! -! Select boundary condition type -! ------------------------------ - select case ( ed % BCWeakType ) - - case ( WEAK_PRESCRIBED ) -! -! Prescribed boundary conditions: just grab the value -! --------------------------------------------------- - F = ed % FB - - case ( WEAK_RIEMANN ) -! -! Weak boundary conditions -! ------------------------ - if ( associated ( ed % RiemannSolver ) ) then - RiemannSolver => ed % RiemannSolver - - else - RiemannSolver => self % RiemannSolver - - end if -! -! Select LEFT and RIGHT states -! ---------------------------- - do iXi = 0 , N - QL = ed % storage(1) % Q(iXi , 1:NCONS) - QR = ed % uB(iXi , 1:NCONS) - - F(iXi , 1:NCONS) = RiemannSolver( QL , QR , ed % n(IX:IY,iXi) ) * ed % dS(iXi) - end do - - end select - - end subroutine StdDG_ComputeRiemannFluxes_CurvedBdry - -end submodule DGInviscid_StandardDG diff --git a/Solver/src/Approximation/DGSpatialDiscretizationMethods.f90 b/Solver/src/Approximation/DGSpatialDiscretizationMethods.f90 index 905567f..c454ac2 100755 --- a/Solver/src/Approximation/DGSpatialDiscretizationMethods.f90 +++ b/Solver/src/Approximation/DGSpatialDiscretizationMethods.f90 @@ -5,6 +5,7 @@ module DGSpatialDiscretizationMethods use DGInviscidMethods use DGWeakIntegrals #ifdef NAVIER_STOKES + use DGArtificialDissipation use DGViscousMethods #endif implicit none @@ -21,7 +22,8 @@ module DGSpatialDiscretizationMethods ! class(InviscidMethod_t), pointer :: InviscidMethod #ifdef NAVIER_STOKES - class(ViscousMethod_t), pointer :: ViscousMethod + class(ViscousMethod_t), pointer :: ViscousMethod + class(ArtificialDissipation_t), pointer :: ArtificialDissipation #endif interface DGSpatial_QDotFaceLoop @@ -53,7 +55,8 @@ subroutine DGSpatial_Initialization() ! Initialize Viscous method ! ------------------------- #ifdef NAVIER_STOKES - ViscousMethod => ViscousMethod_Initialization() + ViscousMethod => ViscousMethod_Initialization() + ArtificialDissipation => ArtificialDissipation_Initialization() #endif end subroutine DGSpatial_Initialization @@ -85,7 +88,7 @@ subroutine DGSpatial_computeTimeDerivative( mesh ) end subroutine DGSpatial_computeTimeDerivative - subroutine DGSpatial_newTimeStep( mesh ) + subroutine DGSpatial_newTimeStep( mesh ) ! ! ************************************************************* ! This subroutine prepares the mesh struct @@ -101,10 +104,6 @@ subroutine DGSpatial_newTimeStep( mesh ) class(QuadMesh_t) :: mesh integer :: zoneID ! -! Reset QDot -! ---------- - call DGSpatial_resetQDot( mesh ) -! ! Interpolate solution to boundaries ! ---------------------------------- call DGSpatial_interpolateSolutionToBoundaries( mesh ) @@ -127,18 +126,6 @@ subroutine DGSpatial_newTimeStep( mesh ) #endif end subroutine DGSpatial_newTimeStep -! - subroutine DGSpatial_resetQDot( mesh ) - implicit none - class(QuadMesh_t) :: mesh -! -------------------------------------- - integer :: eID - - do eID = 1 , mesh % no_of_elements - mesh % elements(eID) % QDot = 0.0_RP - end do - - end subroutine DGSpatial_resetQDot subroutine DGSpatial_computeQDot( mesh ) ! @@ -211,20 +198,24 @@ subroutine DGSpatial_QDotVolumeLoop( e ) class(QuadElement_t) :: e real ( kind=RP ) :: Fi ( 0:e % spA % N , 0:e % spA % N , 1:NCONS , 1:NDIM ) real ( kind=RP ) :: Fv ( 0:e % spA % N , 0:e % spA % N , 1:NCONS , 1:NDIM ) + real ( kind=RP ) :: Fa ( 0:e % spA % N , 0:e % spA % N , 1:NCONS , 1:NDIM ) real ( kind=RP ) :: F ( 0:e % spA % N , 0:e % spA % N , 1:NCONS , 1:NDIM ) real(kind=RP), pointer :: QDot(:,:,:) Fi = InviscidMethod % ComputeInnerFluxes( e ) #ifdef NAVIER_STOKES + e % mu_a = ArtificialDissipation % ComputeElementViscosity ( e ) + Fv = ViscousMethod % ComputeInnerFluxes( e ) + Fa = ArtificialDissipation % ComputeVolumeFluxes( e ) - F = Fi - Fv + F = Fi - Fv - Fa #else F = Fi #endif QDot(0:,0:,1:) => e % QDot - QDot = QDot + ScalarWeakIntegrals % StdVolumeGreen(e , F) + QDot = ScalarWeakIntegrals % StdVolumeGreen(e , F) end subroutine DGSpatial_QDotVolumeLoop @@ -236,22 +227,15 @@ subroutine DGSpatial_QDotFaceLoop_Interior( ed ) real ( kind=RP ) :: FiR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS ) real ( kind=RP ) :: FvL ( 0 : ed % storage ( LEFT ) % spA % N , 1 : NCONS ) real ( kind=RP ) :: FvR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS ) - real ( kind=RP ) :: GvL ( 0 : ed % storage ( LEFT ) % spA % N , 1 : NCONS , 1 : NDIM) - real ( kind=RP ) :: GvR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS , 1 : NDIM) + real ( kind=RP ) :: GL ( 0 : ed % storage ( LEFT ) % spA % N , 1 : NCONS , 1 : NDIM) + real ( kind=RP ) :: GR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS , 1 : NDIM) real ( kind=RP ) :: FL ( 0 : ed % storage ( LEFT ) % spA % N , 1 : NCONS ) real ( kind=RP ) :: FR ( 0 : ed % storage ( RIGHT ) % spA % N , 1 : NCONS ) real ( kind=RP ), pointer :: QDot(:,:,:) - - call InviscidMethod % ComputeRiemannFluxes( ed , FiL , FiR ) -#ifdef NAVIER_STOKES - call ViscousMethod % ComputeRiemannFluxes( ed , FvL , FvR , GvL , GvR ) - FL = FiL - FvL - FR = FiR - FvR -#else - FL = FiL - FR = FiR -#endif - +! +! Compute the Riemann Solver (FL,FR) and the Gradient Riemann Solver (GL,GR) +! -------------------------------------------------------------------------- + call ComputeRiemannSolver_InteriorEdge( ed , FL , FR , GL , GR) ! ! Add the contribution to the LEFT element ! ---------------------------------------- @@ -264,17 +248,20 @@ subroutine DGSpatial_QDotFaceLoop_Interior( ed ) QDot = QDot - ScalarWeakIntegrals % StdFace ( ed , RIGHT , FR ) #ifdef NAVIER_STOKES +! +! Add the gradient fluxes term +! ---------------------------- if ( ViscousMethod % computeRiemannGradientFluxes ) then ! ! Add the contribution to the LEFT element ! ---------------------------------------- QDot(0:,0:,1:) => ed % quads(LEFT) % e % QDot - QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , LEFT , GvL ) + QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , LEFT , GL ) ! ! Add the contribution to the RIGHT element ! ----------------------------------------- QDot(0:,0:,1:) => ed % quads(RIGHT) % e % QDot - QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , RIGHT , GvR ) + QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , RIGHT , GR ) ! end if #endif @@ -286,26 +273,18 @@ subroutine DGSpatial_QDotFaceLoop_StraightBdry( ed ) use QuadElementClass implicit none type(StraightBdryEdge_t) :: ed - real(kind=RP) :: Fi ( 0 : ed % spA % N , 1:NCONS) - real(kind=RP) :: Fv ( 0 : ed % spA % N , 1:NCONS) - real(kind=RP) :: Gv ( 0 : ed % spA % N , 1:NCONS , 1:NDIM ) + real(kind=RP) :: G ( 0 : ed % spA % N , 1:NCONS , 1:NDIM ) real(kind=RP) :: F ( 0 : ed % spA % N , 1:NCONS) real ( kind=RP ), pointer :: QDot(:,:,:) - call InviscidMethod % ComputeRiemannFluxes( ed , Fi ) -#ifdef NAVIER_STOKES - call ViscousMethod % ComputeRiemannFluxes( ed , Fv , Gv) - F = Fi - Fv -#else - F = Fi -#endif + call ComputeRiemannSolver_StraightBdryEdge( ed , F , G ) QDot(0:,0:,1:) => ed % quads(1) % e % QDot QDot = QDot - ScalarWeakIntegrals % StdFace( ed , 1 , F ) #ifdef NAVIER_STOKES if ( ViscousMethod % computeRiemannGradientFluxes ) then - QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , 1 , Gv ) + QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , 1 , G ) end if #endif @@ -316,26 +295,18 @@ subroutine DGSpatial_QDotFaceLoop_CurvedBdry( ed ) use QuadElementClass implicit none type(CurvedBdryEdge_t) :: ed - real(kind=RP) :: Fi ( 0 : ed % spA % N , 1:NCONS) - real(kind=RP) :: Fv ( 0 : ed % spA % N , 1:NCONS) - real(kind=RP) :: Gv ( 0 : ed % spA % N , 1:NCONS , 1:NDIM) + real(kind=RP) :: G ( 0 : ed % spA % N , 1:NCONS , 1:NDIM) real(kind=RP) :: F ( 0 : ed % spA % N , 1:NCONS) real ( kind=RP ), pointer :: QDot(:,:,:) - call InviscidMethod % ComputeRiemannFluxes( ed , Fi ) -#ifdef NAVIER_STOKES - call ViscousMethod % ComputeRiemannFluxes( ed , Fv , Gv ) - F = Fi - Fv -#else - F = Fi -#endif + call ComputeRiemannSolver_CurvedBdryEdge( ed , F , G ) QDot(0:,0:,1:) => ed % quads(1) % e % QDot QDot = QDot - ScalarWeakIntegrals % StdFace( ed , 1 , F ) #ifdef NAVIER_STOKES if ( ViscousMethod % computeRiemannGradientFluxes ) then - QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , 1 , Gv ) + QDot = QDot - ScalarWeakIntegrals % StdGradientFace ( ed , 1 , G ) end if #endif @@ -344,6 +315,788 @@ end subroutine DGSpatial_QDotFaceLoop_CurvedBdry ! !///////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! COMPUTE INTERIOR FLUXES SUBROUTINES +! ----------------------------------- +!///////////////////////////////////////////////////////////////////////////////////////////////////////// +! + subroutine ComputeRiemannSolver_InteriorEdge( ed , FL , FR , GL , GR) +! +! ***************************************************************************** +! This routine computes the Viscous Riemann problem in the "ed" edge. +! -> The result is FStarL and FStarR + +! >> It considers several possibilities: +! +! 1/ LEFT edge needs a p-Transformation and RIGHT edge is reversed. +! 2/ LEFT edge needs a p-Transformation +! 3/ RIGHT edge needs a p-Transformation and RIGHT edge is reversed. +! 4/ RIGHT edge needs a p-Transformation +! 5/ RIGHT edge is reversed. +! 6/ No transformations are needed. +! +! The Riemann solver is computed with the +!> self % RiemannSolver( N , QL , QR , dQL , dQR , normal ) +! procedure. +! +! ***************************************************************************** +! + use QuadElementClass + use MatrixOperations + implicit none + type(Edge_t) :: ed + real(kind=RP) :: FL( 0 : ed % storage(LEFT ) % spA % N , NCONS ) + real(kind=RP) :: FR( 0 : ed % storage(RIGHT) % spA % N , NCONS ) + real(kind=RP), intent (out) :: GL( 0 : ed % storage(LEFT ) % spA % N , 1:NCONS , 1:NDIM) + real(kind=RP), intent (out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1:NCONS , 1:NDIM) +! +! --------------- +! Local variables +! --------------- +! + real ( kind=RP) :: Fi( 0 : ed % spA % N , 1 : NCONS ) + real ( kind=RP) :: Fstar( 0 : ed % spA % N , 1 : NCONS ) + real ( kind=RP ), target :: QL ( 0 : ed % spA % N , 1 : NCONS ) + real ( kind=RP ), target :: QR ( 0 : ed % spA % N , 1 : NCONS ) + real ( kind=RP ), target :: dQL ( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real ( kind=RP ), target :: dQR ( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) +#ifdef NAVIER_STOKES + real ( kind=RP) :: Fv( 0 : ed % spA % N , 1 : NCONS ) + real ( kind=RP) :: Fa( 0 : ed % spA % N , 1 : NCONS ) + real ( kind=RP ) :: GauxL( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM) + real ( kind=RP ) :: GauxR( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM) +#endif + real ( kind=RP ) :: normal(NDIM , 0 : ed % spA % N ) + integer :: eq , iDim +! +! Compute the normal +! ------------------ + normal = spread( ed % n(IX:IY,0) , ncopies = ed % spA % N + 1 , dim = 2 ) +! +! Compute the edge artificial dissipation +! --------------------------------------- +#ifdef NAVIER_STOKES + ed % mu_a = ArtificialDissipation % ComputeEdgeViscosity( ed ) +#endif +! +! Get the solution projection onto the edge +! ----------------------------------------- + call ed % ProjectSolution(ed , QL , QR , dQL , dQR ) +! +! Compute the inviscid Riemann solver +! ----------------------------------- + Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +#ifdef NAVIER_STOKES +! +! Compute the viscous Riemann solver +! ---------------------------------- + Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! +! Compute the artificial dissipation Riemann solver +! ------------------------------------------------- + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! +! The resulting flux is: FStar = ( Inviscid - Viscous - ArtificialDissipation ) dS +! -------------------------------------------------------------------------------- + FStar = ( Fi - Fv - Fa) * ed % dS(0) +#else + FStar = Fi * ed % dS(0) +#endif +! +! Return the resulting Riemann flux to each element frame +! ------------------------------------------------------- + call ed % ProjectFluxes( ed , FStar , FL , FR ) + +#ifdef NAVIER_STOKES + if ( ViscousMethod % computeRiemannGradientFluxes ) then + call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GauxL , GauxR ) + + call ed % ProjectGradientFluxes( ed , GauxL , GauxR , GL , GR ) + + GL = GL * ed % dS(0) + GR = GR * ed % dS(0) + + end if +#endif +! +! +! if ( ed % transform(LEFT) .and. ed % inverse ) then +!! +!! ------------------------------------------------------- +!!> First case: LEFT needs p-Adaption and RIGHT is reversed +!! ------------------------------------------------------- +!! +!! Transform the LEFT edge +!! ----------------------- +! call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) +!#ifdef NAVIER_STOKES +! do eq = 1 , NCONS +! call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % dQ(:,:,eq) , dQL(:,:,eq) ) +! end do +!#endif +!! +!! Invert the RIGHT edge +!! --------------------- +! QR = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , 1:NCONS ) +!#ifdef NAVIER_STOKES +! dQR = ed % storage(RIGHT) % dQ(ed % spA % N : 0 : -1 , 1:NDIM , 1:NCONS ) +!#endif +!! +!! Compute the Riemann solver +!! -------------------------- +! Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +!#ifdef NAVIER_STOKES +! Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! +! FR = ( Fi - Fv - Fa) * ed % dS(0) +!#else +! FR = Fi * ed % dS(0) +!#endif +!! +!! Transform the LEFT edge +!! ----------------------- +! call Mat_x_Mat( ed % T_backward , FR , FL ) +!! +!! Invert the RIGHT edge +!! --------------------- +! FR(0:ed % spA % N , 1:NCONS) = FR(ed % spA % N : 0 : -1 , 1:NCONS) +! +!#ifdef NAVIER_STOKES +!! +!! Compute the Gradient Riemann solver if proceeds +!! ----------------------------------------------- +! if ( ViscousMethod % computeRiemannGradientFluxes ) then +! call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GauxL , GauxR ) +! +! do iDim = 1 , NDIM +! call Mat_x_Mat( ed % T_backward , GauxL(:,:,iDim) , GL(:,:,iDim) ) +! end do +! +! GL = GL * ed % dS(0) +! GR = GauxR(ed % spA % N : 0 : -1 , 1:NCONS , 1:NDIM) * ed % dS(0) +! +! end if +!#endif +! +! elseif ( ed % transform(LEFT) ) then +!! +!! ---------------------------------- +!!> Second case: LEFT needs p-Adaption. +!! ---------------------------------- +!! +!! Transform the LEFT edge +!! ----------------------- +! call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) +!#ifdef NAVIER_STOKES +! do eq = 1 , NCONS +! call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % dQ(:,:,eq) , dQL(:,:,eq) ) +! end do +!#endif +!! +!! Get the RIGHT edge state +!! ------------------------ +! QR = ed % storage(RIGHT) % Q +!#ifdef NAVIER_STOKES +! dQR = ed % storage(RIGHT) % dQ +!#endif +!! +!! Compute the Riemann solver +!! -------------------------- +! Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +!#ifdef NAVIER_STOKES +! Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! FR = ( Fi - Fv - Fa) * ed % dS(0) +!#else +! FR = Fi * ed % dS(0) +!#endif +! +!! +!! Transform the LEFT edge +!! ----------------------- +! call Mat_x_Mat( ed % T_backward , FR , FL ) +! +!#ifdef NAVIER_STOKES +!! +!! Compute the Gradient Riemann solver if proceeds +!! ----------------------------------------------- +! if ( ViscousMethod % computeRiemannGradientFluxes ) then +! call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GauxL , GR ) +! +! do iDim = 1 , NDIM +! call Mat_x_Mat( ed % T_backward , GauxL(:,:,iDim) , GL(:,:,iDim) ) +! end do +! +! GL = GL * ed % dS(0) +! GR = GR * ed % dS(0) +! end if +! +!#endif +! +! elseif ( ed % transform(RIGHT) .and. ed % inverse ) then +!! +!! ---------------------------------------------------- +!!> Third case: RIGHT needs p-Adaption, and its reversed. +!! ---------------------------------------------------- +!! +!! Get the LEFT edge +!! ----------------- +! QL = ed % storage(LEFT) % Q +!#ifdef NAVIER_STOKES +! dQL = ed % storage(LEFT) % dQ +!#endif +!! +!! Transform the RIGHT edge +!! ------------------------ +! call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) +!#ifdef NAVIER_STOKES +! do eq = 1 , NCONS +! call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % dQ(:,:,eq) , dQR(:,:,eq) ) +! end do +!#endif +!! +!! Invert the RIGHT edge +!! --------------------- +! QR(0:ed % spA % N,1:NCONS) = QR(ed % spA % N : 0 : -1 , 1:NCONS) +!#ifdef NAVIER_STOKES +! dQR = dQR(ed % spA % N : 0 : -1 , 1:NDIM , 1:NCONS ) +!#endif +!! +!! Compute the Riemann solver +!! -------------------------- +! Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +!#ifdef NAVIER_STOKES +! Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! FL = ( Fi - Fv - Fa) * ed % dS(0) +!#else +! FL = Fi * ed % dS(0) +!#endif +!! +!! Undo the transformation for the RIGHT edge +!! ------------------------------------------ +! call Mat_x_Mat( ed % T_backward , FL , FR ) +!! +!! Invert the edge +!! --------------- +! FR(0:ed % NLow,1:NCONS) = FR(ed % NLow:0:-1 , 1:NCONS) +! +!#ifdef NAVIER_STOKES +!! +!! Compute the Gradient Riemann solver if proceeds +!! ----------------------------------------------- +! if ( ViscousMethod % computeRiemannGradientFluxes ) then +! call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GL , GauxR ) +! +! do iDim = 1 , NDIM +! call Mat_x_Mat( ed % T_backward , GauxR(:,:,iDim) , GR(:,:,iDim) ) +! end do +! +! GL = GL * ed % dS(0) +! GR = GR(ed % spA % N : 0 : -1 , 1:NCONS, 1:NDIM) * ed % dS(0) +! end if +!#endif +! +! elseif ( ed % transform(RIGHT) ) then +!! +!! ----------------------------------- +!!> Fourth case: RIGHT needs p-Adaption. +!! ----------------------------------- +!! +!! Get the LEFT edge +!! ----------------- +! QL = ed % storage(LEFT) % Q +!#ifdef NAVIER_STOKES +! dQL = ed % storage(LEFT) % dQ +!#endif +!! +!! Transform the RIGHT edge +!! ------------------------ +! call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) +!#ifdef NAVIER_STOKES +! do eq = 1 , NCONS +! call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % dQ(:,:,eq) , dQR(:,:,eq) ) +! end do +!#endif +!! +!! Compute the Riemann solver +!! -------------------------- +! Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +!#ifdef NAVIER_STOKES +! Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! FL = ( Fi - Fv - Fa) * ed % dS(0) +!#else +! FL = Fi * ed % dS(0) +!#endif +!! +!! Undo the transformation for the RIGHT +!! ------------------------------------------ +! call Mat_x_Mat( ed % T_backward , FL , FR ) +! +!#ifdef NAVIER_STOKES +!! +!! Compute the Gradient Riemann solver if proceeds +!! ----------------------------------------------- +! if ( ViscousMethod % computeRiemannGradientFluxes ) then +! call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GL , GauxR ) +! +! do iDim = 1 , NDIM +! call Mat_x_Mat( ed % T_backward , GauxR(:,:,iDim) , GR(:,:,iDim) ) +! end do +! +! GL = GL * ed % dS(0) +! GR = GR * ed % dS(0) +! end if +!#endif +! +! elseif ( ed % inverse ) then +!! +!! ----------------------------- +!!> Fifth case: RIGHT is reversed. +!! ----------------------------- +!! +!! Get the LEFT edge +!! ----------------- +! QL = ed % storage(LEFT) % Q +!#ifdef NAVIER_STOKES +! dQL = ed % storage(LEFT) % dQ +!#endif +!! +!! Invert the RIGHT edge +!! --------------------- +! QR(0:ed % spA % N , 1:NCONS ) = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , 1:NCONS) +!#ifdef NAVIER_STOKES +! dQR = ed % storage(RIGHT) % dQ(ed % spA % N : 0 : -1 , 1:NDIM , 1:NCONS ) +!#endif +!! +!! Compute the Riemann solver +!! -------------------------- +! Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +!#ifdef NAVIER_STOKES +! Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! +! FL = ( Fi - Fv - Fa) * ed % dS(0) +!#else +! FL = Fi * ed % dS(0) +!#endif +!! +!! Invert the RIGHT edge +!! --------------------- +! FR( 0 : ed % spA % N , 1:NCONS ) = FL ( ed % spA % N : 0 : -1 , 1:NCONS ) +! +!#ifdef NAVIER_STOKES +!! +!! Compute the Gradient Riemann solver if proceeds +!! ----------------------------------------------- +! if ( ViscousMethod % computeRiemannGradientFluxes ) then +! call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GL , GR ) +! +! GL = GL * ed % dS(0) +! GR = GR(ed % spA % N : 0 : -1 , 1:NCONS , 1:NDIM) * ed % dS(0) +! end if +!#endif +! +! else +!! +!! ----------------------------------------------------------------------- +!!> Sixth case: Default case: neither p-Adaption nor inversion are required. +!! ----------------------------------------------------------------------- +!! +!! Get the LEFT edge +!! ----------------- +! QL = ed % storage(LEFT) % Q +!#ifdef NAVIER_STOKES +! dQL = ed % storage(LEFT) % dQ +!#endif +!! +!! Get the RIGHT edge +!! ------------------ +! QR = ed % storage(RIGHT) % Q +!#ifdef NAVIER_STOKES +! dQR = ed % storage(RIGHT) % dQ +!#endif +!! +!! Compute the Riemann solver +!! -------------------------- +! Fi = InviscidMethod % RiemannSolver( ed % spA % N , QL , QR , normal ) +!#ifdef NAVIER_STOKES +! Fv = ViscousMethod % RiemannSolver( ed , ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) +! Fa = ArtificialDissipation % ComputeFaceFluxes( ed , QL , QR , dQL , dQR , normal ) +! +! FL = ( Fi - Fv - Fa ) * ed % dS(0) +! FR = FL +!#else +! FL = Fi * ed % dS(0) +! FR = FL +!#endif +! +!#ifdef NAVIER_STOKES +!! +!! Compute the Gradient Riemann solver if proceeds +!! ----------------------------------------------- +! if ( ViscousMethod % computeRiemannGradientFluxes ) then +! call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , QL , QR , normal , GL , GR ) +! +! GL = GL * ed % dS(0) +! GR = GR * ed % dS(0) +! end if +!#endif +! +! end if +!! +!! Change the sign of FR so that it points towards the element outside +!! ------------------------------------------------------------------- +! FR = -1.0_RP * FR + + end subroutine ComputeRiemannSolver_InteriorEdge + + subroutine ComputeRiemannSolver_StraightBdryEdge( ed , F , G ) + use QuadElementClass + implicit none + type(StraightBdryEdge_t) :: ed + real(kind=RP) :: F( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: G( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: Fi( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: Q (0 : ed % spA % N , 1:NCONS) + real(kind=RP) :: Qb(0 : ed % spA % N , 1:NCONS) +#ifdef NAVIER_STOKES + real(kind=RP) :: Fv( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: Fa( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: Gv( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP) :: Gaux( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP) :: dQ (0 : ed % spA % N , 1 : NDIM , 1 : NCONS) + real(kind=RP) :: dQb(0 : ed % spA % N , 1 : NDIM , 1 : NCONS) + real(kind=RP) :: Q1D ( 1 : NCONS ) + real(kind=RP) :: Qb1D ( 1 : NCONS ) + real(kind=RP) :: dQ1D ( 1 : NDIM , 1 : NCONS ) + real(kind=RP) :: dQb1D( 1 : NDIM , 1 : NCONS ) +#endif + real(kind=RP) :: normal( NDIM , 0 : ed % spA % N ) + integer, pointer :: N + integer :: iXi + procedure(RiemannSolverFunction), pointer :: RiemannSolver + + N => ed % spA % N + + normal = spread( ed % n(IX:IY,0) , ncopies = ed % spA % N + 1 , dim = 2 ) +! +! Compute the edge artificial dissipation +! --------------------------------------- +#ifdef NAVIER_STOKES + ed % mu_a = ArtificialDissipation % ComputeEdgeViscosity( ed ) +#endif + +! +! =============== +!> INVISCID FLUXES +! =============== +! +! Select boundary condition type +! ------------------------------ + select case ( ed % inviscidBCType ) + + case ( WEAK_PRESCRIBED ) +! +! Prescribed boundary conditions: just grab the value +! --------------------------------------------------- + Fi = ed % FB + + case ( WEAK_RIEMANN ) +! +! Weak boundary conditions +! ------------------------ + if ( associated ( ed % RiemannSolver ) ) then + RiemannSolver => ed % RiemannSolver + + else + RiemannSolver => InviscidMethod % RiemannSolver + + end if +! +! Select LEFT and RIGHT states +! ---------------------------- + if ( .not. ed % inverse ) then ! This just occurs in periodic BCs + Q = ed % storage(1) % Q + Qb = ed % uB + Fi = RiemannSolver ( ed % spA % N , Q , Qb , normal ) * ed % dS(0) + + else + Q = ed % storage(1) % Q + Qb = ed % uB( ed % spA % N : 0 : -1 , 1:NCONS ) + Fi = RiemannSolver ( ed % spA % N , Q , Qb , normal ) * ed % dS(0) + + end if + end select + +#ifdef NAVIER_STOKES +! +! ============== +!> VISCOUS FLUXES +! ============== +! +!> Select the appropriate boundary condition +! ----------------------------------------- +! + + if ( ed % viscousBCType(0) .eq. PERIODIC ) then +! +! Periodic boundary conditions +! ---------------------------- + if ( ed % inverse ) then + Q = ed % storage(1) % Q + dQ = ed % storage(1) % dQ + + Qb = ed % uB(N:0:-1,1:NCONS) + dQb = ed % gB(N:0:-1,1:NDIM,1:NCONS) + + normal = spread( ed % n(IX:IY,0) , ncopies = N+1 , dim = 2 ) + Fv = ViscousMethod % RiemannSolver( ed , N , ed % invh , Q , Qb , dQ , dQb , normal ) * ed % dS(0) + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , Q , Qb , dQ , dQb , normal ) * ed % dS(0) + + if ( ViscousMethod % computeRiemannGradientFluxes ) then + call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , Q , Qb , normal , Gv , Gaux ) + Gv = Gv * ed % dS(0) + + end if + else + + Q = ed % storage(1) % Q + dQ = ed % storage(1) % dQ + + Qb = ed % uB + dQb = ed % gB + + normal = spread( ed % n(IX:IY,0) , ncopies = N+1 , dim = 2 ) + Fv = ViscousMethod % RiemannSolver( ed , N , ed % invh , Q , Qb , dQ , dQb , normal ) * ed % dS(0) + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , Q , Qb , dQ , dQb , normal ) * ed % dS(0) + + if ( ViscousMethod % computeRiemannGradientFluxes ) then + call ViscousMethod % GradientRiemannSolver ( ed , ed % spA % N , Q , Qb , normal , Gv , Gaux ) + Gv = Gv * ed % dS(0) + + end if + + end if + + elseif ( ed % viscousBCType(0) .eq. ADIABATIC ) then +! +! Adiabatic Dirichlet boundary conditions +! --------------------------------------- + Q = ed % storage(1) % Q + dQ = ed % storage(1) % dQ + Qb = ed % uSB + + normal = spread( ed % n(IX:IY,0) , ncopies = N+1 , dim = 2 ) + + Fv = ViscousMethod % RiemannSolver_Adiabatic( ed , N , ed % invh , Q , dQ , Qb , normal ) * ed % dS(0) + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , Q , Qb , dQ , dQ , normal ) * ed % dS(0) + + if ( ViscousMethod % computeRiemannGradientFluxes ) then + Gv = ViscousMethod % GradientRiemannSolver_Adiabatic( ed , N , Q , Qb , normal ) * ed % dS(0) + + end if + + else +! +! Dirichlet/Neumann boundary conditions +! ------------------------------------- + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , ed % storage(1) % Q , ed % uSB , ed % storage(1) % dQ , ed % storage(1) % dQ , normal ) * ed % dS(0) + do iXi = 0 , ed % spA % N + select case ( ed % viscousBCType(iXi) ) + + case ( DIRICHLET ) + + Q1D = ed % storage(1) % Q(iXi,:) + dQ1D = ed % storage(1) % dQ(iXi,:,:) + Qb1D = ed % uSB(iXi,:) + Fv(iXi,:) = ViscousMethod % RiemannSolver_Dirichlet ( ed , ed % spA % N , ed % invh , Q1D , dQ1D , Qb1D , ed % n(IX:IY,0) ) * ed % dS(0) + + if ( ViscousMethod % computeRiemannGradientFluxes ) then + Gv(iXi,:,:) = ViscousMethod % GradientRiemannSolver_BoundaryCondition( ed , Q1D , Qb1D , ed % n(IX:IY,0) ) * ed % dS(0) + end if + + case ( NEUMANN ) + + Fv(iXi,:) = 0.0_RP + Fa(iXi,:) = 0.0_RP + Gv(iXi,:,:) = 0.0_RP + + end select + end do + + + end if +#endif + +#ifdef NAVIER_STOKES + !Fa = 0.0_RP + F = Fi - Fv - Fa + G = Gv +#else + F = Fi + G = 0.0_RP +#endif + + end subroutine ComputeRiemannSolver_StraightBdryEdge + + subroutine ComputeRiemannSolver_CurvedBdryEdge( ed , F , G ) + use QuadElementClass + implicit none + type(CurvedBdryEdge_t) :: ed + real(kind=RP) :: F( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: G( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: Fi( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: Q(0 : ed % spA % N , 1:NCONS) + real(kind=RP) :: Qb(0 : ed % spA % N , 1:NCONS) +#ifdef NAVIER_STOKES + real(kind=RP) :: Fv( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: Fa( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP) :: Gv( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP) :: Gaux( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP) :: dQ (0 : ed % spA % N , 1 : NDIM , 1 : NCONS) + real(kind=RP) :: dQb(0 : ed % spA % N , 1 : NDIM , 1 : NCONS) + real(kind=RP) :: Q1D ( 1 : NCONS ) + real(kind=RP) :: Qb1D ( 1 : NCONS ) + real(kind=RP) :: dQ1D ( 1 : NDIM , 1 : NCONS ) + real(kind=RP) :: dQb1D( 1 : NDIM , 1 : NCONS ) +#endif + integer, pointer :: N + integer :: iXi + procedure(RiemannSolverFunction), pointer :: RiemannSolver + + N => ed % spA % N +! +! Compute the edge artificial dissipation +! --------------------------------------- +#ifdef NAVIER_STOKES + ed % mu_a = ArtificialDissipation % ComputeEdgeViscosity( ed ) +#endif +! +! =============== +!> INVISCID FLUXES +! =============== +! +! Select boundary condition type +! ------------------------------ + select case ( ed % inviscidBCType ) + + case ( WEAK_PRESCRIBED ) +! +! Prescribed boundary conditions: just grab the value +! --------------------------------------------------- + Fi = ed % FB + + case ( WEAK_RIEMANN ) +! +! Weak boundary conditions +! ------------------------ + if ( associated ( ed % RiemannSolver ) ) then + RiemannSolver => ed % RiemannSolver + + else + RiemannSolver => InviscidMethod % RiemannSolver + + end if +! +! Select LEFT and RIGHT states +! ---------------------------- + if ( .not. ed % inverse ) then ! This just occurs in periodic BCs + Q = ed % storage(1) % Q + Qb = ed % uB + Fi = RiemannSolver ( ed % spA % N , Q , Qb , ed % n ) * spread( ed % dS , ncopies = NCONS , dim = 2 ) + + else + Q = ed % storage(1) % Q + Qb = ed % uB( ed % spA % N : 0 : -1 , 1:NCONS ) + Fi = RiemannSolver ( ed % spA % N , Q , Qb , ed % n ) * spread( ed % dS , ncopies = NCONS , dim = 2 ) + + end if + end select + +#ifdef NAVIER_STOKES +! +! ============== +!> VISCOUS FLUXES +! ============== +! +!> Select the appropriate boundary condition +! ----------------------------------------- +! + if ( ed % viscousBCType(0) .eq. ADIABATIC ) then +! +! Adiabatic Dirichlet boundary conditions +! --------------------------------------- + Q = ed % storage(1) % Q + dQ = ed % storage(1) % dQ + Qb = ed % uSB + + Fv = ViscousMethod % RiemannSolver_Adiabatic( ed , N , ed % invh , Q , dQ , Qb , ed % n ) * spread ( ed % dS , ncopies = NCONS , dim = 2 ) + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , Q , Qb , dQ , dQ , ed % n ) * spread( ed % dS , ncopies = NCONS , dim = 2 ) + + if ( ViscousMethod % computeRiemannGradientFluxes ) then + Gv = ViscousMethod % GradientRiemannSolver_Adiabatic( ed , N , Q , Qb , ed % n ) * spread( spread ( ed % dS , ncopies = NCONS , dim = 2 ) , ncopies = NDIM , dim = 3) + + end if + + else +! +! Dirichlet/Neumann boundary conditions +! ------------------------------------- + Fa = ArtificialDissipation % ComputeFaceFluxes( ed , ed % storage(1) % Q , ed % uSB , ed % storage(1) % dQ , ed % storage(1) % dQ , ed % n ) * spread( ed % dS , ncopies = NCONS , dim = 2 ) + do iXi = 0 , ed % spA % N + select case ( ed % viscousBCType(iXi) ) + + case ( DIRICHLET ) + + Q1D = ed % storage(1) % Q(iXi,:) + dQ1D = ed % storage(1) % dQ(iXi,:,:) + Qb1D = ed % uSB(iXi,:) + Fv(iXi,:) = ViscousMethod % RiemannSolver_Dirichlet ( ed , ed % spA % N , ed % invh , Q1D , dQ1D , Qb1D , ed % n(IX:IY,iXi) ) * ed % dS(iXi) + + if ( ViscousMethod % computeRiemannGradientFluxes ) then + Gv(iXi,:,:) = ViscousMethod % GradientRiemannSolver_BoundaryCondition( ed , Q1D , Qb1D , ed % n(IX:IY,iXi) ) * ed % dS(iXi) + end if + + case ( NEUMANN ) + + Fv(iXi,:) = 0.0_RP + Fa(iXi,:) = 0.0_RP + Gv(iXi,:,:) = 0.0_RP + + end select + + + end do + + end if +#endif + +#ifdef NAVIER_STOKES + !Fa = 0.0_RP + F = Fi - Fv - Fa + G = Gv +#else + F = Fi + G = 0.0_RP +#endif + + end subroutine ComputeRiemannSolver_CurvedBdryEdge + +! +!///////////////////////////////////////////////////////////////////////////////////////////////////////// +! ! INTERPOLATE TO BOUNDARIES SUBROUTINES ! ------------------------------------- !///////////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/Solver/src/Approximation/DGTimeIntegrator.f90 b/Solver/src/Approximation/DGTimeIntegrator.f90 index f948c66..df9501b 100755 --- a/Solver/src/Approximation/DGTimeIntegrator.f90 +++ b/Solver/src/Approximation/DGTimeIntegrator.f90 @@ -129,7 +129,6 @@ function TimeIntegrator_NewTimeIntegrator(mesh) result (Integrator) Monitors = ConstructMonitors( mesh ) end if - end function TimeIntegrator_NewTimeIntegrator subroutine TimeIntegrator_Integrate( self , mesh , Storage ) diff --git a/Solver/src/Approximation/DGViscousMethods.f90 b/Solver/src/Approximation/DGViscousMethods.f90 index 8266de0..4a1ea48 100755 --- a/Solver/src/Approximation/DGViscousMethods.f90 +++ b/Solver/src/Approximation/DGViscousMethods.f90 @@ -37,8 +37,8 @@ module DGViscousMethods ! --------------------------------------------------------------- ! type ViscousMethod_t - logical :: computeRiemannGradientFluxes - character(len=STR_LEN_VISCOUS) :: method + logical :: computeRiemannGradientFluxes + character(len=STR_LEN_VISCOUS) :: method contains ! ! *************************************** @@ -69,23 +69,17 @@ module DGViscousMethods ! RIEMANN SOLVER PROCEDURE ! *************************************** ! - procedure, private :: RiemannSolver => BaseClass_RiemannSolver - procedure, private :: RiemannSolver_Dirichlet => BaseClass_RiemannSolver_Dirichlet - procedure, private :: RiemannSolver_Adiabatic => BaseClass_RiemannSolver_Adiabatic - generic, public :: ComputeRiemannFluxes => ComputeRiemannFluxes_Interior , & - ComputeRiemannFluxes_StraightBdry , & - ComputeRiemannFluxes_CurvedBdry - procedure, private :: ComputeRiemannFluxes_Interior => BaseClass_ComputeRiemannFluxes_Interior - procedure, private :: ComputeRiemannFluxes_StraightBdry => BaseClass_ComputeRiemannFluxes_StraightBdry - procedure, private :: ComputeRiemannFluxes_CurvedBdry => BaseClass_ComputeRiemannFluxes_CurvedBdry + procedure :: RiemannSolver => BaseClass_RiemannSolver + procedure :: RiemannSolver_Dirichlet => BaseClass_RiemannSolver_Dirichlet + procedure :: RiemannSolver_Adiabatic => BaseClass_RiemannSolver_Adiabatic ! ! *************************************** ! GRADIENT RIEMANN SOLVER PROCEDURE ! *************************************** ! - procedure, private :: GradientRiemannSolver => BaseClass_GradientRiemannSolver - procedure, private :: GradientRiemannSolver_BoundaryCondition => BaseClass_GradientRiemannSolver_BoundaryCondition - procedure, private :: GradientRiemannSolver_Adiabatic => BaseClass_GradientRiemannSolver_Adiabatic + procedure :: GradientRiemannSolver => BaseClass_GradientRiemannSolver + procedure :: GradientRiemannSolver_BoundaryCondition => BaseClass_GradientRiemannSolver_BoundaryCondition + procedure :: GradientRiemannSolver_Adiabatic => BaseClass_GradientRiemannSolver_Adiabatic procedure :: Describe => ViscousMethod_describe end type ViscousMethod_t ! @@ -118,9 +112,9 @@ module DGViscousMethods procedure :: ComputeGradient => BR1_ComputeGradient procedure :: ComputeInnerFluxes => BR1_ComputeInnerFluxes procedure, private :: SolutionRiemannSolver => BR1_SolutionRiemannSolver - procedure, private :: RiemannSolver => BR1_RiemannSolver - procedure, private :: RiemannSolver_Dirichlet => BR1_RiemannSolver_Dirichlet - procedure, private :: RiemannSolver_Adiabatic => BR1_RiemannSolver_Adiabatic + procedure :: RiemannSolver => BR1_RiemannSolver + procedure :: RiemannSolver_Dirichlet => BR1_RiemannSolver_Dirichlet + procedure :: RiemannSolver_Adiabatic => BR1_RiemannSolver_Adiabatic end type BR1Method_t ! ! --------------------------------------------------------------- @@ -162,9 +156,10 @@ module pure function BR1_SolutionRiemannSolver( self , N , UL , UR ) result ( uS real(kind=RP) :: uStar(0:N,1:NCONS) end function BR1_SolutionRiemannSolver - module pure function BR1_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) + module pure function BR1_RiemannSolver( self , edge , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) implicit none class(BR1Method_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) @@ -175,9 +170,10 @@ module pure function BR1_RiemannSolver( self , N , invh_edge , UL , UR , dUL , d real(kind=RP) :: Fstar(0:N , 1:NCONS) end function BR1_RiemannSolver - module pure function BR1_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function BR1_RiemannSolver_Dirichlet( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) implicit none class(BR1Method_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer , intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(NCONS) @@ -187,10 +183,11 @@ module pure function BR1_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , real(kind=RP) :: Fstar(1:NCONS) end function BR1_RiemannSolver_Dirichlet - module pure function BR1_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function BR1_RiemannSolver_Adiabatic( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) implicit none class(BR1Method_t), intent(in) :: self integer , intent(in) :: N + class(Edge_t) , intent(in) :: edge real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(0:N , NCONS) real(kind=RP), intent(in) :: g(0:N , NDIM,NCONS) @@ -220,9 +217,10 @@ module pure function IP_ComputeInnerFluxes( self , e ) result (Fv) real(kind=RP) :: Fv(0 : e % spA % N , 0 : e % spA % N , 1:NCONS , 1:NDIM) end function IP_ComputeInnerFluxes - module pure function IP_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) + module pure function IP_RiemannSolver( self , edge , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) @@ -233,9 +231,10 @@ module pure function IP_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dU real(kind=RP) :: Fstar(0:N , 1:NCONS) end function IP_RiemannSolver - module pure function IP_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function IP_RiemannSolver_Dirichlet( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer , intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(NCONS) @@ -245,9 +244,10 @@ module pure function IP_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , real(kind=RP) :: Fstar(1:NCONS) end function IP_RiemannSolver_Dirichlet - module pure function IP_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function IP_RiemannSolver_Adiabatic( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) implicit none class(IPMethod_t), intent (in) :: self + class(Edge_t) , intent(in) :: edge integer , intent (in) :: N real(kind=RP), intent (in) :: invh_edge real(kind=RP), intent (in) :: u ( 0:N , NCONS ) @@ -257,9 +257,10 @@ module pure function IP_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , real(kind=RP) :: Fstar ( 0:N , 1:NCONS ) end function IP_RiemannSolver_Adiabatic - module pure subroutine IP_GradientRiemannSolver( self , N , UL , UR , normal , GstarL , GstarR ) + module pure subroutine IP_GradientRiemannSolver( self , edge , N , UL , UR , normal , GstarL , GstarR ) implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) real(kind=RP), intent(in) :: uR(0:N , 1:NCONS) @@ -268,18 +269,20 @@ module pure subroutine IP_GradientRiemannSolver( self , N , UL , UR , normal , G real(kind=RP), intent(out) :: GstarR(0:N , 1:NCONS , 1:NDIM) end subroutine IP_GradientRiemannSolver - module pure function IP_GradientRiemannSolver_BoundaryCondition( self , u , uB , normal ) result ( Gstar ) + module pure function IP_GradientRiemannSolver_BoundaryCondition( self , edge , u , uB , normal ) result ( Gstar ) implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge real(kind=RP), intent(in) :: u(1:NCONS) real(kind=RP), intent(in) :: uB(1:NCONS) real(kind=RP), intent(in) :: normal(IX:IY) real(kind=RP) :: Gstar(1:NCONS , 1:NDIM) end function IP_GradientRiemannSolver_BoundaryCondition - module pure function IP_GradientRiemannSolver_Adiabatic( self , N , u , uB , normal ) result ( Gstar ) + module pure function IP_GradientRiemannSolver_Adiabatic( self , edge , N , u , uB , normal ) result ( Gstar ) implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: u (0:N , 1:NCONS) real(kind=RP), intent(in) :: uB(0:N , 1:NCONS) @@ -412,9 +415,10 @@ pure function BaseClass_SolutionRiemannSolver( self , N , uL , uR ) result ( uSt ! end function BaseClass_SolutionRiemannSolver - pure function BaseClass_RiemannSolver( self , N , invh_edge , uL , uR , duL , duR , normal ) result ( FStar ) + pure function BaseClass_RiemannSolver( self , edge , N , invh_edge , uL , uR , duL , duR , normal ) result ( FStar ) implicit none class(ViscousMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) @@ -430,9 +434,10 @@ pure function BaseClass_RiemannSolver( self , N , invh_edge , uL , uR , duL , du ! end function BaseClass_RiemannSolver - pure function BaseClass_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + pure function BaseClass_RiemannSolver_Dirichlet( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) implicit none class(ViscousMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer , intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(NCONS) @@ -447,9 +452,10 @@ pure function BaseClass_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , ! end function BaseClass_RiemannSolver_Dirichlet - pure function BaseClass_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + pure function BaseClass_RiemannSolver_Adiabatic( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) implicit none class(ViscousMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer , intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(0:N , NCONS) @@ -464,13 +470,14 @@ pure function BaseClass_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , ! end function BaseClass_RiemannSolver_Adiabatic - pure subroutine BaseClass_GradientRiemannSolver( self , N , UL , UR , normal , GstarL , GstarR ) + pure subroutine BaseClass_GradientRiemannSolver( self , edge , N , UL , UR , normal , GstarL , GstarR ) ! ! ***************************************************************************************** ! ***************************************************************************************** ! implicit none class(ViscousMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) real(kind=RP), intent(in) :: uR(0:N , 1:NCONS) @@ -479,18 +486,20 @@ pure subroutine BaseClass_GradientRiemannSolver( self , N , UL , UR , normal , G real(kind=RP), intent(out) :: GstarR(0:N , 1:NCONS , 1:NDIM) end subroutine BaseClass_GradientRiemannSolver - pure function BaseClass_GradientRiemannSolver_BoundaryCondition( self , u , uB , normal ) result ( Gstar ) + pure function BaseClass_GradientRiemannSolver_BoundaryCondition( self , edge , u , uB , normal ) result ( Gstar ) implicit none class(ViscousMethod_t), intent(in) :: self + class(Edge_t), intent(in) :: edge real(kind=RP), intent(in) :: u(1:NCONS) real(kind=RP), intent(in) :: uB(1:NCONS) real(kind=RP), intent(in) :: normal(1:NDIM) real(kind=RP) :: Gstar(1:NCONS , 1:NDIM) end function BaseClass_GradientRiemannSolver_BoundaryCondition - pure function BaseClass_GradientRiemannSolver_Adiabatic( self , N , u , uB , normal ) result ( Gstar ) + pure function BaseClass_GradientRiemannSolver_Adiabatic( self , edge , N , u , uB , normal ) result ( Gstar ) implicit none class(ViscousMethod_t), intent(in) :: self + class(Edge_t), intent(in) :: edge integer , intent(in) :: N real(kind=RP), intent(in) :: u(0:N,1:NCONS) real(kind=RP), intent(in) :: uB(0:N,1:NCONS) @@ -801,612 +810,6 @@ subroutine BaseClass_ComputeSolutionRiemann_CurvedBdry( self , ed , uStar) end subroutine BaseClass_ComputeSolutionRiemann_CurvedBdry ! -!///////////////////////////////////////////////////////////////////////////////////////////////////////// -! -! VISCOUS FLUXES SUBROUTINES -! -------------------------- -!///////////////////////////////////////////////////////////////////////////////////////////////////////// -! -!TODO Pure - subroutine BaseClass_ComputeRiemannFluxes_Interior( self , ed , FStarL , FStarR , GStarL , GStarR ) -! -! ***************************************************************************** -! This routine computes the Viscous Riemann problem in the "ed" edge. -! -> The result is FStarL and FStarR -! >> Recall that viscous fluxes are defined just for RHOU,RHOV, and RHOE -! states. -! -! This procedure is common to all ViscousMethods, since it just provides -! a framework to deal with the data sharing across the interface. Then, -! each method considers its own Riemann Solver procedure -! ------------------ -! >> Therefore, it considers several possibilities: -! -! 1/ LEFT edge needs a p-Transformation and RIGHT edge is reversed. -! 2/ LEFT edge needs a p-Transformation -! 3/ RIGHT edge needs a p-Transformation and RIGHT edge is reversed. -! 4/ RIGHT edge needs a p-Transformation -! 5/ RIGHT edge is reversed. -! 6/ No transformations are needed. -! -! The Riemann solver is computed with the -!> self % RiemannSolver( N , QL , QR , dQL , dQR , normal ) -! procedure. -! -! ***************************************************************************** -! - - use QuadElementClass - use MatrixOperations - implicit none - class(ViscousMethod_t), intent(in) :: self - type(Edge_t), intent(in) :: ed - real(kind=RP), intent(out) :: FStarL( 0 : ed % storage(LEFT ) % spA % N , 1:NCONS ) - real(kind=RP), intent(out) :: FStarR( 0 : ed % storage(RIGHT) % spA % N , 1:NCONS ) - real(kind=RP), intent(out) :: GStarL( 0 : ed % storage(LEFT ) % spA % N , 1:NCONS , 1:NDIM) - real(kind=RP), intent(out) :: GStarR( 0 : ed % storage(RIGHT) % spA % N , 1:NCONS , 1:NDIM) -! -! --------------- -! Local variables -! --------------- -! - real ( kind=RP ), target :: QL ( 0 : ed % spA % N , 1 : NCONS ) - real ( kind=RP ), target :: QR ( 0 : ed % spA % N , 1 : NCONS ) - real ( kind=RP ), target :: dQL ( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) - real ( kind=RP ), target :: dQR ( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) - real ( kind=RP ) :: GStarL_aux( 0 : ed % spA % N , 1 : NCONS , 1:NDIM) - real ( kind=RP ) :: GStarR_aux( 0 : ed % spA % N , 1 : NCONS , 1:NDIM) - real ( kind=RP ) :: normal(1:NDIM , 0 : ed % spA % N ) - integer :: iXi , eq , iDim -! -! ---------------------------------------------------------------------------------------------- -!> Straight boundaries are considered, so replicate the normal in the set of interpolation points -! ---------------------------------------------------------------------------------------------- -! - do iXi = 0 , ed % spA % N - normal(IX:IY,iXi) = ed % n(IX:IY,0) !spread( ed % n(IX:IY,0) , ncopies = ed % spA % N+1 , dim = 2 ) - end do - - if ( ed % transform(LEFT) .and. ed % inverse ) then -! -! ------------------------------------------------------- -!> First case: LEFT needs p-Adaption and RIGHT is reversed -! ------------------------------------------------------- -! -! Transform the LEFT edge -! ----------------------- - call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) - do eq = 1 , NCONS - call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % dQ(:,:,eq) , dQL(:,:,eq) ) - end do -! -! Invert the RIGHT edge -! --------------------- - QR = ed % storage(RIGHT) % Q (ed % spA % N : 0 : -1 , 1:NCONS ) - dQR = ed % storage(RIGHT) % dQ(ed % spA % N : 0 : -1 , 1:NDIM , 1:NCONS ) -! -! Compute the Riemann solver -! -------------------------- - FStarR = self % RiemannSolver( ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) * ed % dS(0) - -! -! Transform the LEFT edge -! ----------------------- - call Mat_x_Mat( ed % T_backward , FStarR , FStarL ) -! -! Invert the RIGHT edge -! --------------------- - FStarR(0:ed % spA % N , 1:NCONS) = FStarR(ed % spA % N : 0 : -1 , 1:NCONS) -! -! Compute the Gradient Riemann solver if proceeds -! ----------------------------------------------- - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , QL , QR , normal , GstarL_aux , GstarR_aux ) - - do iDim = 1 , NDIM - call Mat_x_Mat( ed % T_backward , GstarL_aux(:,:,iDim) , GstarL(:,:,iDim) ) - end do - GstarL = GstarL * ed % dS(0) - GStarR = GStarR_aux(ed % spA % N : 0 : -1 , 1:NCONS , 1:NDIM) * ed % dS(0) - end if - - elseif ( ed % transform(LEFT) ) then -! -! ---------------------------------- -!> Second case: LEFT needs p-Adaption. -! ---------------------------------- -! -! Transform the LEFT edge -! ----------------------- - call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) - do eq = 1 , NCONS - call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % dQ(:,:,eq) , dQL(:,:,eq) ) - end do -! -! Get the RIGHT edge state -! ------------------------ - QR = ed % storage(RIGHT) % Q - dQR = ed % storage(RIGHT) % dQ -! -! Compute the Riemann solver -! -------------------------- - FStarR = self % RiemannSolver( ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) * ed % dS(0) -! -! Transform the LEFT edge -! ----------------------- - call Mat_x_Mat( ed % T_backward , FStarR , FStarL ) -! -! Compute the Gradient Riemann solver if proceeds -! ----------------------------------------------- - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , QL , QR , normal , GstarL_aux , GstarR ) - - do iDim = 1 , NDIM - call Mat_x_Mat( ed % T_backward , GstarL_aux(:,:,iDim) , GstarL(:,:,iDim) ) - end do - - GstarL = GstarL * ed % dS(0) - GstarR = GstarR * ed % dS(0) - end if - - - elseif ( ed % transform(RIGHT) .and. ed % inverse ) then -! -! ---------------------------------------------------- -!> Third case: RIGHT needs p-Adaption, and its reversed. -! ---------------------------------------------------- -! -! Get the LEFT edge -! ----------------- - QL = ed % storage(LEFT) % Q - dQL = ed % storage(LEFT) % dQ -! -! Transform the RIGHT edge -! ------------------------ - call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) - do eq = 1 , NCONS - call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % dQ(:,:,eq) , dQR(:,:,eq) ) - end do -! -! Invert the RIGHT edge -! --------------------- - QR(0:ed % spA % N,1:NCONS) = QR(ed % spA % N : 0 : -1 , 1:NCONS) - dQR = dQR(ed % spA % N : 0 : -1 , 1:NDIM , 1:NCONS ) -! -! Compute the Riemann solver -! -------------------------- - FStarL = self % RiemannSolver( ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) * ed % dS(0) -! -! Undo the transformation for the RIGHT edge -! ------------------------------------------ - call Mat_x_Mat( ed % T_backward , FStarL , FStarR ) -! -! Invert the RIGHT edge -! --------------------- - FStarR(0:ed % NLow,1:NCONS) = FStarR(ed % NLow:0:-1 , 1:NCONS) -! -! Compute the Gradient Riemann solver if proceeds -! ----------------------------------------------- - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , QL , QR , normal , GstarL , GstarR_aux ) - - do iDim = 1 , NDIM - call Mat_x_Mat( ed % T_backward , GstarR_aux(:,:,iDim) , GstarR(:,:,iDim) ) - end do - GstarL = GstarL * ed % dS(0) - GStarR = GStarR(ed % spA % N : 0 : -1 , 1:NCONS, 1:NDIM) * ed % dS(0) - end if - - - elseif ( ed % transform(RIGHT) ) then -! -! ----------------------------------- -!> Fourth case: RIGHT needs p-Adaption. -! ----------------------------------- -! -! Get the LEFT edge -! ----------------- - QL = ed % storage(LEFT) % Q - dQL = ed % storage(LEFT) % dQ -! -! Transform the RIGHT edge -! ------------------------ - call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) - do eq = 1 , NCONS - call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % dQ(:,:,eq) , dQR(:,:,eq) ) - end do -! -! Compute the Riemann solver -! -------------------------- - FStarL = self % RiemannSolver( ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) * ed % dS(0) -! -! Undo the transformation for the RIGHT edge -! ------------------------------------------ - call Mat_x_Mat( ed % T_backward , FStarL , FStarR ) -! -! Compute the Gradient Riemann solver if proceeds -! ----------------------------------------------- - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , QL , QR , normal , GstarL , GstarR_aux ) - - do iDim = 1 , NDIM - call Mat_x_Mat( ed % T_backward , GstarR_aux(:,:,iDim) , GstarR(:,:,iDim) ) - end do - - GstarL = GstarL * ed % dS(0) - GStarR = GStarR * ed % dS(0) - end if - - - elseif ( ed % inverse ) then -! -! ----------------------------- -!> Fifth case: RIGHT is reversed. -! ----------------------------- -! -! Get the LEFT edge -! ----------------- - QL = ed % storage(LEFT) % Q - dQL = ed % storage(LEFT) % dQ -! -! Invert the RIGHT edge -! --------------------- - QR = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , 1:NCONS) - dQR = ed % storage(RIGHT) % dQ(ed % spA % N : 0 : -1 , 1:NDIM , 1:NCONS ) -! -! Compute the Riemann solver -! -------------------------- - FStarL = self % RiemannSolver( ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) * ed % dS(0) -! -! Invert the RIGHT edge -! --------------------- - FStarR( 0 : ed % spA % N , 1:NCONS ) = FStarL ( ed % spA % N : 0 : -1 , 1:NCONS ) -! -! Compute the Gradient Riemann solver if proceeds -! ----------------------------------------------- - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , QL , QR , normal , GstarL , GstarR ) - - GstarL = GstarL * ed % dS(0) - GStarR = GStarR(ed % spA % N : 0 : -1 , 1:NCONS , 1:NDIM) * ed % dS(0) - end if - - else -! -! ----------------------------------------------------------------------- -!> Sixth case: Default case: neither p-Adaption nor inversion are required. -! ----------------------------------------------------------------------- -! -! Get the LEFT edge -! ----------------- - QL = ed % storage(LEFT) % Q - dQL = ed % storage(LEFT) % dQ -! -! Get the RIGHT edge -! ------------------ - QR = ed % storage(RIGHT) % Q - dQR = ed % storage(RIGHT) % dQ -! -! Compute the Riemann solver -! -------------------------- - FStarL = self % RiemannSolver( ed % spA % N , ed % invh , QL , QR , dQL , dQR , normal ) * ed % dS(0) - FStarR = FStarL -! -! Compute the Gradient Riemann solver if proceeds -! ----------------------------------------------- - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , QL , QR , normal , GstarL , GstarR ) - - GstarL = GstarL * ed % dS(0) - GStarR = GStarR * ed % dS(0) - end if - - end if -! -! ------------------------------------------------------------------- -!> Change the sign of FR so that it points towards the element outside -! ------------------------------------------------------------------- -! - FStarR = -1.0_RP * FStarR - - end subroutine BaseClass_ComputeRiemannFluxes_Interior - - subroutine BaseClass_ComputeRiemannFluxes_StraightBdry( self , ed , FStar , GStar) -! -! ***************************************************************************** -! This routine computes the Viscous Riemann problem in the "ed" edge, which -! is a Straight boundary. -! -> The result is FStar - -! >> Recall that viscous fluxes are defined just for RHOU,RHOV, and RHOE -! states. -! -! >> The Riemann solver depends on the Boundary condition type: -! -! 1/ Periodic BCs: self % RiemannSolver (Same than interior edges) -! 2/ Dirichlet BCs: self % RiemannSolver_Dirichlet -! 3/ Neumann BCs: FStar = 0.0_RP -! 4/ Adiabatic BCs: self % RiemannSolver_Adiabatic -! -! ***************************************************************************** -! - - use QuadElementClass - implicit none - class(ViscousMethod_t) , intent(in) :: self - type(StraightBdryEdge_t), intent(in) :: ed - real(kind=RP) , intent(out) :: FStar( 0 : ed % spA % N , 1:NCONS ) - real(kind=RP) , intent(out) :: GStar( 0 : ed % spA % N , 1:NCONS , 1:NDIM) -! -! --------------- -! Local variables -! --------------- -! - real ( kind=RP ) :: Q ( 0 : ed % spA % N , 1:NCONS ) - real ( kind=RP ) :: dQ ( 0 : ed % spA % N , 1:NDIM , 1:NCONS ) - real ( kind=RP ) :: Qb ( 0 : ed % spa % N , 1:NCONS ) - real ( kind=RP ) :: dQb ( 0 : ed % spA % N , 1:NDIM,1:NCONS ) - real ( kind=RP ) :: normal ( 1 : NDIM , 0 : ed % spA % N ) - real ( kind=RP ) :: Q1D ( 1 : NCONS ) - real ( kind=RP ) :: dQ1D ( 1 : NDIM , 1 : NCONS ) - real ( kind=RP ) :: Qb1D ( 1 : NCONS ) - real(kind=RP) :: Gaux( 0 : ed % spA % N , 1:NCONS , 1:NDIM) - integer :: N - integer :: iXi - - N = ed % spA % N -! -! ----------------------------------------- -!> Select the appropriate boundary condition -! ----------------------------------------- -! - if ( ed % viscousBCType(0) .eq. PERIODIC ) then -! -! Periodic boundary conditions -! ---------------------------- - if ( ed % inverse ) then - Q = ed % storage(1) % Q - dQ = ed % storage(1) % dQ - - Qb = ed % uB(N:0:-1,1:NCONS) - dQb = ed % gB(N:0:-1,1:NDIM,1:NCONS) - - normal = spread( ed % n(IX:IY,0) , ncopies = N+1 , dim = 2 ) - Fstar = self % RiemannSolver( N , ed % invh , Q , Qb , dQ , dQb , normal ) * ed % dS(0) - - if ( self % computeRiemannGradientFluxes ) then - call self % GradientRiemannSolver ( ed % spA % N , Q , Qb , normal , GStar , Gaux ) - Gstar = GStar * ed % dS(0) - -! do iXi = 0 , ed % spA % N -! Q1D = Q(iXi,:) -! Qb1D = Qb(iXi,:) -! Gstar(iXi,:,:) = 0.5_RP * self % GradientRiemannSolver_BoundaryCondition( Q1D , Qb1D , ed % n(IX:IY,0) ) * ed % dS(0) -! -! end do - end if - else - - Q = ed % storage(1) % Q - dQ = ed % storage(1) % dQ - - Qb = ed % uB - dQb = ed % gB - - normal = spread( ed % n(IX:IY,0) , ncopies = N+1 , dim = 2 ) - Fstar = self % RiemannSolver( N , ed % invh , Q , Qb , dQ , dQb , normal ) * ed % dS(0) - - if ( self % computeRiemannGradientFluxes ) then - - call self % GradientRiemannSolver ( ed % spA % N , Q , Qb , normal , GStar , Gaux ) - Gstar = GStar * ed % dS(0) -! do iXi = 0 , ed % spA % N -! Q1D = Q(iXi,:) -! Qb1D = Qb(iXi,:) -! Gstar(iXi,:,:) = 0.5_RP * self % GradientRiemannSolver_BoundaryCondition( Q1D , Qb1D , ed % n(IX:IY,0) ) * ed % dS(0) -! -! end do - end if - - end if - - elseif ( ed % viscousBCType(0) .eq. ADIABATIC ) then -! -! Adiabatic Dirichlet boundary conditions -! --------------------------------------- - Q = ed % storage(1) % Q - dQ = ed % storage(1) % dQ - Qb = ed % uSB - - normal = spread( ed % n(IX:IY,0) , ncopies = N+1 , dim = 2 ) - - Fstar = self % RiemannSolver_Adiabatic( N , ed % invh , Q , dQ , Qb , normal ) * ed % dS(0) - - if ( self % computeRiemannGradientFluxes ) then - Gstar = self % GradientRiemannSolver_Adiabatic( N , Q , Qb , normal ) * ed % dS(0) - end if - - else -! -! Dirichlet/Neumann boundary conditions -! ------------------------------------- - do iXi = 0 , ed % spA % N - select case ( ed % viscousBCType(iXi) ) - - case ( DIRICHLET ) - - Q1D = ed % storage(1) % Q(iXi,:) - dQ1D = ed % storage(1) % dQ(iXi,:,:) - Qb1D = ed % uSB(iXi,:) - Fstar(iXi,:) = self % RiemannSolver_Dirichlet ( ed % spA % N , ed % invh , Q1D , dQ1D , Qb1D , ed % n(IX:IY,0) ) * ed % dS(0) - - if ( self % computeRiemannGradientFluxes ) then - Gstar(iXi,:,:) = self % GradientRiemannSolver_BoundaryCondition( Q1D , Qb1D , ed % n(IX:IY,0) ) * ed % dS(0) - end if - - case ( NEUMANN ) - - Fstar(iXi,:) = 0.0_RP - Gstar(iXi,:,:) = 0.0_RP - - end select - - - end do - - end if - - end subroutine BaseClass_ComputeRiemannFluxes_StraightBdry - - pure subroutine BaseClass_ComputeRiemannFluxes_CurvedBdry( self , ed , Fstar , Gstar) -! -! ***************************************************************************** -! This routine computes the Viscous Riemann problem in the "ed" edge, which -! is a Curved boundary. -! -> The result is FStar - -! >> Recall that viscous fluxes are defined just for RHOU,RHOV, and RHOE -! states. -! -! >> The Riemann solver depends on the Boundary condition type: -! -! 1/ Periodic BCs: self % RiemannSolver (Same than interior edges) -! 2/ Dirichlet BCs: self % RiemannSolver_Dirichlet -! 3/ Neumann BCs: FStar = 0.0_RP -! 4/ Adiabatic BCs: self % RiemannSolver_Adiabatic -! -! ***************************************************************************** -! - use QuadElementClass - implicit none - class(ViscousMethod_t) , intent (in) :: self - type(CurvedBdryEdge_t), intent (in) :: ed - real(kind=RP) , intent (out) :: FStar( 0 : ed % spA % N , 1:NCONS ) - real(kind=RP) , intent (out) :: GStar( 0 : ed % spA % N , 1:NCONS , 1:NDIM) - -! -! --------------- -! Local variables -! --------------- -! - real ( kind=RP ) :: Q ( 0 : ed % spA % N , 1:NCONS ) - real ( kind=RP ) :: dQ ( 0 : ed % spA % N , 1:NDIM , 1:NCONS ) - real ( kind=RP ) :: Qb ( 0 : ed % spa % N , 1:NCONS ) - real ( kind=RP ) :: dQb ( 0 : ed % spA % N , 1:NDIM,1:NCONS ) - real ( kind=RP ) :: Q1D ( 1 : NCONS ) - real ( kind=RP ) :: dQ1D ( 1 : NDIM , 1 : NCONS ) - real ( kind=RP ) :: Qb1D ( 1 : NCONS ) - integer :: N - integer :: iXi - integer(kind=1) :: eq - - N = ed % spA % N -! -! ----------------------------------------- -!> Select the appropriate boundary condition -! ----------------------------------------- -! - if ( ed % viscousBCType(0) .eq. PERIODIC ) then -! -! Periodic boundary conditions -! ---------------------------- - if ( ed % inverse ) then - Q = ed % storage(1) % Q - dQ = ed % storage(1) % dQ - - Qb = ed % uSB(N:0:-1,1:NCONS) - dQb = ed % gB(N:0:-1,1:NDIM,1:NCONS) - - Fstar = self % RiemannSolver( N , ed % invh , Q , Qb , dQ , dQb , ed % n ) - - do eq = 1 , NCONS - Fstar(:,eq) = Fstar(:,eq) * ed % dS - end do - - if ( self % computeRiemannGradientFluxes ) then - do iXi = 0 , ed % spA % N - Q1D = Q(iXi,:) - Qb1D = Qb(iXi,:) - Gstar(iXi,:,:) = self % GradientRiemannSolver_BoundaryCondition( Q1D , Qb1D , ed % n(IX:IY,iXi) ) * ed % dS(iXi) - - end do - end if - - else - Q = ed % storage(1) % Q - dQ = ed % storage(1) % dQ - - Qb = ed % uSB - dQb = ed % gB - - Fstar = self % RiemannSolver( N , ed % invh , Q , Qb , dQ , dQb , ed % n ) - - do eq = 1 , NCONS - Fstar(:,eq) = Fstar(:,eq) * ed % dS - end do - - if ( self % computeRiemannGradientFluxes ) then - do iXi = 0 , ed % spA % N - Q1D = Q(iXi,:) - Qb1D = Qb(iXi,:) - Gstar(iXi,:,:) = self % GradientRiemannSolver_BoundaryCondition( Q1D , Qb1D , ed % n(IX:IY,iXi) ) * ed % dS(iXi) - - end do - end if - - - end if - - elseif ( ed % viscousBCType(0) .eq. ADIABATIC ) then -! -! Adiabatic Dirichlet boundary conditions -! --------------------------------------- - Q = ed % storage(1) % Q - dQ = ed % storage(1) % dQ - Qb = ed % uSB - - Fstar = self % RiemannSolver_Adiabatic( N , ed % invh , Q , dQ , Qb , ed % n ) - - do eq = 1 , NCONS - Fstar(:,eq) = Fstar(:,eq) * ed % dS - end do - - if ( self % computeRiemannGradientFluxes ) then - Gstar = self % GradientRiemannSolver_Adiabatic( N , Q , Qb , ed % n ) - do iXi = 0 , N - Gstar(iXi,:,:) = Gstar(iXi,:,:) * ed % dS(iXi) - end do - end if - - else -! -! Dirichlet/Neumann boundary conditions -! ------------------------------------- - do iXi = 0 , ed % spA % N - select case ( ed % viscousBCType(iXi) ) - - case ( DIRICHLET ) - - Q1D = ed % storage(1) % Q(iXi,:) - dQ1D = ed % storage(1) % dQ(iXi,:,:) - Qb1D = ed % uSB(iXi,:) - Fstar(iXi,:) = self % RiemannSolver_Dirichlet ( ed % spA % N , ed % invh , Q1D , dQ1D , Qb1D , ed % n(IX:IY,iXi) ) * ed % dS(iXi) - - if ( self % computeRiemannGradientFluxes ) then - Gstar(iXi,:,:) = self % GradientRiemannSolver_BoundaryCondition( Q1D , Qb1D , ed % n(IX:IY,iXi) ) * ed % dS(iXi) - end if - - case ( NEUMANN ) - - Fstar(iXi,:) = 0.0_RP - Gstar(iXi,:,:) = 0.0_RP - - end select - end do - end if - - end subroutine BaseClass_ComputeRiemannFluxes_CurvedBdry -! !////////////////////////////////////////////////////////////////////////////////////////////////// ! ! AUXILIAR PROCEDURES diff --git a/Solver/src/Approximation/DGViscous_BR1.f90 b/Solver/src/Approximation/DGViscous_BR1.f90 index eda4364..5d9390e 100644 --- a/Solver/src/Approximation/DGViscous_BR1.f90 +++ b/Solver/src/Approximation/DGViscous_BR1.f90 @@ -127,7 +127,13 @@ subroutine BR1_dQFaceLoop_CurvedBdry( ViscousMethod , ed ) dQ = dQ + VectorWeakIntegrals % StdFace( ed , 1 , ed % uSB ) ! end subroutine BR1_dQFaceLoop_CurvedBdry - +! +!////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +! +! TIME DERIVATIVE PROCEDURES +! -------------------------- +!////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +! module pure function BR1_ComputeInnerFluxes( self , e ) result (Fv) ! ! *************************************************************************************** @@ -145,31 +151,9 @@ module pure function BR1_ComputeInnerFluxes( self , e ) result (Fv) class(QuadElement_t), intent(in) :: e real(kind=RP) :: Fv(0 : e % spA % N , 0 : e % spA % N , 1:NCONS , 1:NDIM) ! -! --------------- -! Local variables -! --------------- -! - real(kind=RP) :: F_cartesian(0:e % spA % N,0:e % spA % N,1:NCONS,1:NDIM) - integer :: eq - integer :: N - - N = e % spA % N -! ! Compute the cartesian flux ! -------------------------- - F_cartesian = ViscousFlux( e % spA % N , e % Q , e % dQ) - - do eq = 1 , NCONS -! -! F flux (contravariant) -! ---------------------- - Fv(0:N,0:N,eq,IX) = F_cartesian(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,1) + F_cartesian(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,1) -! -! G flux (contravariant) -! ---------------------- - Fv(0:N,0:N,eq,IY) = F_cartesian(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,2) + F_cartesian(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,2) - end do - + Fv = (dimensionless % mu + e % mu_a) * ViscousFlux( e % spA % N , e % Q , e % dQ) end function BR1_ComputeInnerFluxes @@ -192,7 +176,7 @@ module pure function BR1_SolutionRiemannSolver( self , N , UL , UR ) result ( uS end function BR1_SolutionRiemannSolver - module pure function BR1_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) + module pure function BR1_RiemannSolver( self , edge , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) ! ! ***************************************************************************************** ! The BR1 Viscous Riemann Solver averages the fluxes obtained from both sides. @@ -203,6 +187,7 @@ module pure function BR1_RiemannSolver( self , N , invh_edge , UL , UR , dUL , d ! implicit none class(BR1Method_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) @@ -222,8 +207,8 @@ module pure function BR1_RiemannSolver( self , N , invh_edge , UL , UR , dUL , d ! ! Compute the LEFT and RIGHT viscous fluxes ! ----------------------------------------- - FL = ViscousFlux( N , uL , duL ) - FR = ViscousFlux( N , uR , duR ) + FL = ( dimensionless % mu + edge % mu_a ) * ViscousFlux( N , uL , duL ) + FR = ( dimensionless % mu + edge % mu_a ) * ViscousFlux( N , uR , duR ) ! ! Perform the average and the projection along the edge normal ! ------------------------------------------------------------ @@ -233,7 +218,7 @@ module pure function BR1_RiemannSolver( self , N , invh_edge , UL , UR , dUL , d end function BR1_RiemannSolver - module pure function BR1_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function BR1_RiemannSolver_Dirichlet( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) ! ! ***************************************************************************************** ! For the Dirichlet boundary conditions, the BR1 scheme uses the interior values @@ -241,6 +226,7 @@ module pure function BR1_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , ! implicit none class(BR1Method_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer , intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(NCONS) @@ -257,7 +243,7 @@ module pure function BR1_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , ! ! Compute the two dimensional flux ! -------------------------------- - F = ViscousFlux( u , g ) + F = (dimensionless % mu + edge % mu_a ) * ViscousFlux( u , g ) ! ! Projection along the boundary normal ! ------------------------------------ @@ -265,7 +251,7 @@ module pure function BR1_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , end function BR1_RiemannSolver_Dirichlet - module pure function BR1_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function BR1_RiemannSolver_Adiabatic( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) ! ! ************************************************************************************************ ! For the Adiabatic dirichlet boundary conditions, the BR1 scheme uses the interior values @@ -273,6 +259,7 @@ module pure function BR1_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , ! implicit none class(BR1Method_t), intent (in) :: self + class(Edge_t) , intent(in) :: edge integer , intent (in) :: N real(kind=RP), intent (in) :: invh_edge real(kind=RP), intent (in) :: u ( 0:N , NCONS ) @@ -290,7 +277,7 @@ module pure function BR1_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , ! ! Compute the Adiabatic viscous flux based on the interior points ! --------------------------------------------------------------- - F = AdiabaticViscousFlux( N , u , u , g) + F = ( dimensionless % mu + edge % mu_a ) * AdiabaticViscousFlux( N , u , u , g) ! ! Perform the projection along the boundary normal ! ------------------------------------------------ diff --git a/Solver/src/Approximation/DGViscous_IP.f90 b/Solver/src/Approximation/DGViscous_IP.f90 index 122329b..095205a 100644 --- a/Solver/src/Approximation/DGViscous_IP.f90 +++ b/Solver/src/Approximation/DGViscous_IP.f90 @@ -37,7 +37,7 @@ module subroutine IP_ComputeGradient( self , mesh ) ! Just compute the interior gradient ! ---------------------------------- do eID = 1 , mesh % no_of_elements - call mesh % elements(eID) % ComputeInteriorGradient + mesh % elements(eID) % dQ = mesh % elements(eID) % ComputeInteriorGradient() end do end subroutine IP_ComputeGradient @@ -59,35 +59,13 @@ module pure function IP_ComputeInnerFluxes( self , e ) result (Fv) class(QuadElement_t), intent(in) :: e real(kind=RP) :: Fv(0 : e % spA % N , 0 : e % spA % N , 1:NCONS , 1:NDIM) ! -! --------------- -! Local variables -! --------------- -! - real(kind=RP) :: F_cartesian(0:e % spA % N,0:e % spA % N,1:NCONS,1:NDIM) - integer :: eq - integer :: N - - N = e % spA % N -! ! Compute the cartesian flux ! -------------------------- - F_cartesian = ViscousFlux( e % spA % N , e % Q , e % dQ) - - do eq = 1 , NCONS -! -! F flux (contravariant) -! ---------------------- - Fv(0:N,0:N,eq,IX) = F_cartesian(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,1) + F_cartesian(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,1) -! -! G flux (contravariant) -! ---------------------- - Fv(0:N,0:N,eq,IY) = F_cartesian(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,2) + F_cartesian(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,2) - end do - + Fv = (dimensionless % mu + e % mu_a) * ViscousFlux( e % spA % N , e % Q , e % dQ) end function IP_ComputeInnerFluxes - module pure function IP_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) + module pure function IP_RiemannSolver( self , edge , N , invh_edge , UL , UR , dUL , dUR , normal ) result ( FStar ) ! ! ***************************************************************************************** ! The IP Viscous Riemann Solver averages the fluxes obtained from both sides. @@ -98,6 +76,7 @@ module pure function IP_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dU ! implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) @@ -118,13 +97,13 @@ module pure function IP_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dU ! ! Compute the penalty parameter ! ----------------------------- - penalty = self % sigma0 * dimensionless % mu * N * N * invh_edge + penalty = self % sigma0 * (dimensionless % mu + edge % mu_a) * N * N * invh_edge ! ! Compute the LEFT and RIGHT viscous fluxes ! ----------------------------------------- - FL = ViscousFlux( N , uL , duL ) - FR = ViscousFlux( N , uR , duR ) + FL = (dimensionless % mu + edge % mu_a) * ViscousFlux( N , uL , duL ) + FR = (dimensionless % mu + edge % mu_a) * ViscousFlux( N , uR , duR ) ! ! Perform the average and the projection along the edge normal ! ------------------------------------------------------------ @@ -137,7 +116,7 @@ module pure function IP_RiemannSolver( self , N , invh_edge , UL , UR , dUL , dU end function IP_RiemannSolver - module pure function IP_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function IP_RiemannSolver_Dirichlet( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) ! ! ***************************************************************************************** ! For the Dirichlet boundary conditions, the IP scheme uses the interior values @@ -145,6 +124,7 @@ module pure function IP_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , ! implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: invh_edge real(kind=RP), intent(in) :: u(NCONS) @@ -158,22 +138,22 @@ module pure function IP_RiemannSolver_Dirichlet( self , N , invh_edge , u , g , ! --------------- ! real(kind=RP) :: F(1:NCONS,1:NDIM) - real(kind=RP) :: penalty + real(kind=RP) :: penalty ! ! Compute the two dimensional flux ! -------------------------------- - F = ViscousFlux( u , g ) + F = (dimensionless % mu + edge % mu_a) * ViscousFlux( u , g ) ! ! Projection along the boundary normal ! ------------------------------------ FStar = F(:,IX) * normal(IX) + F(:,IY) * normal(IY) - penalty = self % sigma0 * dimensionless % mu * N * N * invh_edge + penalty = self % sigma0 * (dimensionless % mu + edge % mu_a) * N * N * invh_edge Fstar = Fstar - penalty * ( u - uB ) end function IP_RiemannSolver_Dirichlet - module pure function IP_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , uB , normal ) result ( Fstar ) + module pure function IP_RiemannSolver_Adiabatic( self , edge , N , invh_edge , u , g , uB , normal ) result ( Fstar ) ! ! ************************************************************************************************ ! For the Adiabatic dirichlet boundary conditions, the IP scheme uses the interior values @@ -181,6 +161,7 @@ module pure function IP_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , ! implicit none class(IPMethod_t), intent (in) :: self + class(Edge_t) , intent(in) :: edge integer , intent (in) :: N real(kind=RP), intent (in) :: invh_edge real(kind=RP), intent (in) :: u ( 0:N , NCONS ) @@ -194,12 +175,12 @@ module pure function IP_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , ! --------------- ! real(kind=RP) :: F(0:N , 1:NCONS , 1:NDIM) - real(kind=RP) :: penalty + real(kind=RP) :: penalty integer :: eq ! ! Compute the Adiabatic viscous flux based on the interior points ! --------------------------------------------------------------- - F = AdiabaticViscousFlux( N , u , u , g) + F = (dimensionless % mu + edge % mu_a) * AdiabaticViscousFlux( N , u , u , g) ! ! Perform the projection along the boundary normal ! ------------------------------------------------ @@ -207,18 +188,19 @@ module pure function IP_RiemannSolver_Adiabatic( self , N , invh_edge , u , g , FStar(:,eq) = F(:,eq,IX) * normal(IX,:) + F(:,eq,IY) * normal(IY,:) end do - penalty = self % sigma0 * dimensionless % mu * N * N * invh_edge + penalty = self % sigma0 * (dimensionless % mu + edge % mu_a) * N * N * invh_edge Fstar = Fstar - penalty * ( u - uB ) end function IP_RiemannSolver_Adiabatic - module pure subroutine IP_GradientRiemannSolver( self , N , UL , UR , normal , GstarL , GstarR ) + module pure subroutine IP_GradientRiemannSolver( self , edge , N , UL , UR , normal , GstarL , GstarR ) ! ! ***************************************************************************************** ! ***************************************************************************************** ! implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: uL(0:N , 1:NCONS) real(kind=RP), intent(in) :: uR(0:N , 1:NCONS) @@ -242,18 +224,19 @@ module pure subroutine IP_GradientRiemannSolver( self , N , UL , UR , normal , G ! ! Compute the fluxes ! ------------------ - GStarL = 0.5_RP * self % epsilon * ViscousFlux( N , uL , falseGradient ) - GStarR = -0.5_RP * self % epsilon * ViscousFlux( N , uR , falseGradient ) + GStarL = 0.5_RP * self % epsilon * (dimensionless % mu + edge % mu_a) * ViscousFlux( N , uL , falseGradient ) + GStarR = -0.5_RP * self % epsilon * (dimensionless % mu + edge % mu_a) * ViscousFlux( N , uR , falseGradient ) end subroutine IP_GradientRiemannSolver - module pure function IP_GradientRiemannSolver_BoundaryCondition( self , u , uB , normal ) result ( Gstar ) + module pure function IP_GradientRiemannSolver_BoundaryCondition( self , edge , u , uB , normal ) result ( Gstar ) ! ! ***************************************************************************************** ! ***************************************************************************************** ! implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge real(kind=RP), intent(in) :: u(1:NCONS) real(kind=RP), intent(in) :: uB(1:NCONS) real(kind=RP), intent(in) :: normal(IX:IY) @@ -273,17 +256,18 @@ module pure function IP_GradientRiemannSolver_BoundaryCondition( self , u , uB , ! ! Compute the fluxes ! ------------------ - GStar = 1.0_RP * self % epsilon * ViscousFlux( u , falseGradient ) + GStar = 1.0_RP * self % epsilon * (dimensionless % mu + edge % mu_a) * ViscousFlux( u , falseGradient ) end function IP_GradientRiemannSolver_BoundaryCondition - module pure function IP_GradientRiemannSolver_Adiabatic( self , N , u , uB , normal ) result ( Gstar ) + module pure function IP_GradientRiemannSolver_Adiabatic( self , edge , N , u , uB , normal ) result ( Gstar ) ! ! ***************************************************************************************** ! ***************************************************************************************** ! implicit none class(IPMethod_t), intent(in) :: self + class(Edge_t) , intent(in) :: edge integer, intent(in) :: N real(kind=RP), intent(in) :: u(0:N , 1:NCONS) real(kind=RP), intent(in) :: uB(0:N , 1:NCONS) @@ -306,9 +290,7 @@ module pure function IP_GradientRiemannSolver_Adiabatic( self , N , u , uB , nor ! ! Compute the fluxes ! ------------------ -! GStar = self % epsilon * AdiabaticViscousFlux( N , u , u , falseGradient ) - !GStar = self % epsilon * ViscousFlux( N , u , falseGradient ) - Gstar = 0.0_RP + GStar = self % epsilon * (dimensionless % mu + edge % mu_a) * AdiabaticViscousFlux( N , u , u , falseGradient ) end function IP_GradientRiemannSolver_Adiabatic end submodule DGViscous_IP diff --git a/Solver/src/Approximation/DGWeakIntegrals.f90 b/Solver/src/Approximation/DGWeakIntegrals.f90 index cf5b457..db7d510 100644 --- a/Solver/src/Approximation/DGWeakIntegrals.f90 +++ b/Solver/src/Approximation/DGWeakIntegrals.f90 @@ -115,12 +115,33 @@ function Scalar_StdVolumeGreen( e , F) result ( volInt ) real(kind=RP), intent(in) :: F(0:e % spA % N , 0:e % spA % N, 1:NCONS , 1:NDIM) real(kind=RP) :: volInt(0:e % spA % N, 0:e % spA % N , 1:NCONS) ! ---------------------------------------------------------------------------------------------- + real(kind=RP) :: contravariant_F( 0 : e % spA % N , 0 : e % spA % N , 1:NCONS ) + real(kind=RP) :: contravariant_G( 0 : e % spA % N , 0 : e % spA % N , 1:NCONS ) integer, pointer :: N + integer :: eq N => e % spA % N + + do eq = 1 , NCONS +! +! F flux (contravariant) +! ---------------------- + contravariant_F(0:N,0:N,eq) = F(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,1) + F(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,1) + + end do - volInt = MatrixMultiplyInIndex_F(F(:,:,:,IX) , e % spA % hatD , N+1 , N+1 , NCONS , IX) & - + MatrixMultiplyInIndex_F(F(:,:,:,IY) , e % spA % hatD , N+1 , N+1 , NCONS , IY) + do eq = 1 , NCONS +! +! G flux (contravariant) +! ---------------------- + contravariant_G(0:N,0:N,eq) = F(0:N,0:N,eq,IX) * e % Ja(0:N,0:N,1,2) + F(0:N,0:N,eq,IY) * e % Ja(0:N,0:N,2,2) + + end do + + + + volInt = MatrixMultiplyInIndex_F(contravariant_F , e % spA % hatD , N+1 , N+1 , NCONS , IX) & + + MatrixMultiplyInIndex_F(contravariant_G , e % spA % hatD , N+1 , N+1 , NCONS , IY) end function Scalar_StdVolumeGreen diff --git a/Solver/src/BoundaryConditions/Dirichlet.incf b/Solver/src/BoundaryConditions/Dirichlet.incf index 513e30b..05f7113 100644 --- a/Solver/src/BoundaryConditions/Dirichlet.incf +++ b/Solver/src/BoundaryConditions/Dirichlet.incf @@ -59,11 +59,10 @@ #ifdef _DIMENSIONLESS_TAU rho = pressure / Temperature - V = sqrt(thermodynamics % gamma) * Mach #else rho = dimensionless % gammaMach2 * pressure / Temperature - V = Mach * sqrt( thermodynamics % Gamma * pressure / rho ) #endif + V = Mach * sqrt( thermodynamics % Gamma * pressure / rho ) ! ! Construct the state vector ! -------------------------- @@ -94,7 +93,7 @@ type is (StraightBdryEdge_t) allocate ( edge % uB ( 0:N , NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -118,7 +117,7 @@ type is (CurvedBdryEdge_t) allocate ( edge % uB ( 0:N , NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/EulerWall.incf b/Solver/src/BoundaryConditions/EulerWall.incf index 87c3478..ccbd387 100644 --- a/Solver/src/BoundaryConditions/EulerWall.incf +++ b/Solver/src/BoundaryConditions/EulerWall.incf @@ -21,7 +21,7 @@ allocate ( edge % FB ( 0:N , NCONS ) ) allocate ( edge % uB ( 0:N , NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -35,7 +35,7 @@ allocate ( edge % FB ( 0:N , NCONS ) ) allocate ( edge % uB ( 0:N , NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/Farfield.incf b/Solver/src/BoundaryConditions/Farfield.incf index 1ecaa16..c130756 100644 --- a/Solver/src/BoundaryConditions/Farfield.incf +++ b/Solver/src/BoundaryConditions/Farfield.incf @@ -90,7 +90,7 @@ type is (StraightBdryEdge_t) allocate( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -102,7 +102,7 @@ type is (CurvedBdryEdge_t) allocate( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/PressureInlet.incf b/Solver/src/BoundaryConditions/PressureInlet.incf index 3b63fcb..587684c 100644 --- a/Solver/src/BoundaryConditions/PressureInlet.incf +++ b/Solver/src/BoundaryConditions/PressureInlet.incf @@ -88,7 +88,7 @@ allocate ( edge % FB ( 0:N , 1:NCONS ) ) allocate ( edge % uB ( 0:N , 1:NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -100,7 +100,7 @@ allocate ( edge % FB ( 0:N , 1:NCONS ) ) allocate ( edge % uB ( 0:N , 1:NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/PressureOutlet.incf b/Solver/src/BoundaryConditions/PressureOutlet.incf index 06e9f10..1a13ea6 100644 --- a/Solver/src/BoundaryConditions/PressureOutlet.incf +++ b/Solver/src/BoundaryConditions/PressureOutlet.incf @@ -101,7 +101,7 @@ type is (StraightBdryEdge_t) allocate( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -113,7 +113,7 @@ type is (CurvedBdryEdge_t) allocate( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/Riemann.incf b/Solver/src/BoundaryConditions/Riemann.incf index 4114d9c..f2f3b26 100644 --- a/Solver/src/BoundaryConditions/Riemann.incf +++ b/Solver/src/BoundaryConditions/Riemann.incf @@ -104,7 +104,7 @@ allocate ( edge % FB(0:N,NCONS) ) allocate ( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -116,7 +116,7 @@ allocate ( edge % FB(0:N,NCONS) ) allocate ( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/ViscousWall.incf b/Solver/src/BoundaryConditions/ViscousWall.incf index b0d34f0..4364b33 100644 --- a/Solver/src/BoundaryConditions/ViscousWall.incf +++ b/Solver/src/BoundaryConditions/ViscousWall.incf @@ -95,7 +95,7 @@ type is (StraightBdryEdge_t) allocate( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -114,7 +114,7 @@ type is (CurvedBdryEdge_t) allocate( edge % uB(0:N,NCONS) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/BoundaryConditions/newDirichlet.incf b/Solver/src/BoundaryConditions/newDirichlet.incf index b548352..505a6a6 100644 --- a/Solver/src/BoundaryConditions/newDirichlet.incf +++ b/Solver/src/BoundaryConditions/newDirichlet.incf @@ -116,7 +116,7 @@ allocate ( edge % FB ( 0:N , 1:NCONS ) ) allocate ( edge % uB ( 0:N , 1:NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver @@ -125,7 +125,7 @@ allocate ( edge % FB ( 0:N , 1:NCONS ) ) allocate ( edge % uB ( 0:N , 1:NCONS ) ) - edge % BCWeakType = self % WeakType + edge % inviscidBCType = self % WeakType edge % RiemannSolver => self % RiemannSolver diff --git a/Solver/src/Generic/Setup.f90 b/Solver/src/Generic/Setup.f90 index 6293b89..480b8bf 100755 --- a/Solver/src/Generic/Setup.f90 +++ b/Solver/src/Generic/Setup.f90 @@ -63,6 +63,15 @@ module Setup_class real(kind=RP), allocatable :: sigma1IP ! ! ------------------------------------------------------------------------------ +! Artificial dissipation +! ------------------------------------------------------------------------------ +! + logical :: artificialDissipation + real(kind=RP), allocatable :: artificialDissipationIntensity + character(len=STR_LEN_SETUP) :: artificialDissipationIndicator + character(len=STR_LEN_SETUP) :: artificialDissipationType +! +! ------------------------------------------------------------------------------ ! Integration parameters ! ------------------------------------------------------------------------------ ! @@ -86,6 +95,7 @@ module Setup_class character(len=STR_LEN_SETUP) :: solution_file character(len=STR_LEN_SETUP) :: restart_file character(len=STR_LEN_SETUP) :: outputType + character(len=STR_LEN_SETUP) :: exportFormat contains procedure, nopass :: Initialization => Setup_Initialization @@ -113,6 +123,7 @@ subroutine Setup_Initialization character(len=STR_LEN_SETUP) :: case_name character(len=STR_LEN_SETUP) :: interp_nodes character(len=STR_LEN_SETUP) :: inviscid_form + integer, allocatable :: artificialDissipation ! ! Get case file from command line @@ -295,6 +306,27 @@ subroutine Setup_Initialization call readValue ( trim ( case_name ) , "Integration mode" , Setup % integrationMode ) call Setup_CheckWithDefault( Setup % integrationMode , "Steady" , "Integration mode" ) ! +! Request the artificial dissipation +! ---------------------------------- + call readValue ( trim ( case_name ) , "Artificial dissipation (0/1)" , artificialDissipation ) + call Setup_CheckWithDefault ( artificialDissipation , 0 , "Artificial dissipation (0/1)" ) + if ( artificialDissipation .eq. 1 ) then + Setup % artificialDissipation = .true. + + else + Setup % artificialDissipation = .false. + + end if + + call readValue ( trim ( case_name ) , "Artificial dissipation intensity" , Setup % artificialDissipationIntensity ) + call Setup_CheckWithDefault ( Setup % artificialDissipationIntensity , 1.0_RP , "Artificial dissipation intensity" ) + + call readValue ( trim ( case_name ) , "Artificial dissipation indicator" , Setup % artificialDissipationIndicator ) + call Setup_CheckWithDefault ( Setup % artificialDissipationIndicator , "Jumps-based" , "Artificial dissipation indicator" ) + + call readValue ( trim ( case_name ) , "Artificial dissipation type" , Setup % artificialDissipationType ) + call Setup_CheckWithDefault ( Setup % artificialDissipationType , "Physical" , "Artificial dissipation type" ) +! ! Request the integration scheme ! ------------------------------ call readValue ( trim ( case_name ) , "Integration scheme" , Setup % integrationMethod ) @@ -363,6 +395,11 @@ subroutine Setup_Initialization ! ------------------------------- call readValue ( trim ( case_name ) , "Number of representation points" , Setup % no_of_plotPoints ) call Setup_CheckWithDefault( Setup % no_of_plotPoints , 2 * Setup % N , "Number of representation points" ) +! +! Export format +! ------------- + call readValue ( trim ( case_name ) , "Export format" , Setup % exportFormat ) + call Setup_CheckWithDefault( Setup % exportFormat , "Tecplot" , "Export format" ) end subroutine Setup_Initialization diff --git a/Solver/src/Generic/Storage.f90 b/Solver/src/Generic/Storage.f90 index 038a43f..3effb0d 100755 --- a/Solver/src/Generic/Storage.f90 +++ b/Solver/src/Generic/Storage.f90 @@ -42,9 +42,12 @@ subroutine Storage_AllocateMemory( self , totalPolynomialOrder ) ! dQ -> NGRAD * NDIM * (N+1) * (N+1) * no_of_elements ! ------------------------------------------------------------- allocate ( self % Q ( NCONS * totalPolynomialOrder ) ) + self % Q = 0.0_RP allocate ( self % QDot ( NCONS * totalPolynomialOrder ) ) + self % QDot = 0.0_RP #ifdef NAVIER_STOKES allocate ( self % dQ ( NCONS * NDIM * totalPolynomialOrder ) ) + self % dQ = 0.0_RP #endif end subroutine Storage_AllocateMemory diff --git a/Solver/src/IO/ProblemFile.f90 b/Solver/src/IO/ProblemFile.f90 index b822458..d28a745 100755 --- a/Solver/src/IO/ProblemFile.f90 +++ b/Solver/src/IO/ProblemFile.f90 @@ -13,20 +13,28 @@ function UserDefinedInitialCondition(x , argin) result (val) associate ( gamma => thermodynamics % gamma , Mach => dimensionless % Mach ) - val(IRHO) = 1.0_RP - val(IRHOV) = 0.0_RP - val(IRHOE) = dimensionless % cv * 1.0_RP + 0.5_RP * gamma * Mach * Mach - - if ( x(IY) .lt. 0.0_RP ) then - val(IRHOU) = sqrt(gamma) * Mach - else - val(IRHOU) = -sqrt(gamma) * Mach - - end if - - val(IRHOU) = val(IRHOU) + 1.0e-3_RP * sqrt(gamma) * Mach * exp(-(x(IY) / (28.0_RP))**2.0_RP) * ( cos(2.0_RP * PI * x(IX)) + cos(20.0_RP * PI * x(IX) ) ) - - +! val(IRHO) = 1.0_RP +! val(IRHOV) = 0.0_RP +! val(IRHOE) = dimensionless % cv * 1.0_RP + 0.5_RP * gamma * Mach * Mach +! +! if ( x(IY) .lt. 0.0_RP ) then +! val(IRHOU) = sqrt(gamma) * Mach +! else +! val(IRHOU) = -sqrt(gamma) * Mach +! +! end if +! +! val(IRHOU) = val(IRHOU) + 1.0e-3_RP * sqrt(gamma) * Mach * exp(-(x(IY) / (28.0_RP))**2.0_RP) * ( cos(2.0_RP * PI * x(IX)) + cos(20.0_RP * PI * x(IX) ) ) +! + + if ( x(IX) .lt. 0.0_RP ) then + val = [0.537037037037037_RP, 0.829539386177859_RP, 0.0_RP, 1.657627118644068_RP] + + else + val = [1.000000000000000_RP, 0.829539386177859_RP, 0.0_RP, 2.844067796610170_RP] + + + end if end associate end function UserDefinedInitialCondition diff --git a/Solver/src/Mesh/QuadElement.f90 b/Solver/src/Mesh/QuadElement.f90 index 3500b6e..1207c46 100755 --- a/Solver/src/Mesh/QuadElement.f90 +++ b/Solver/src/Mesh/QuadElement.f90 @@ -10,6 +10,7 @@ module QuadElementClass private public QuadElement_t , QuadElement_p , Edge_t , StraightBdryEdge_t , CurvedBdryEdge_t , Edge_p + public Edge_ProjectSolutionType1 !//////////////////////////////////////////////////////////////////////////////////////////////////////// ! @@ -24,6 +25,7 @@ module QuadElementClass integer :: edgesDirection(EDGES_PER_QUAD) ! Direction (FORWARD/REVERSE) of the edges integer :: quadPosition(EDGES_PER_QUAD) ! Position of the quad for the edge (LEFT/RIGHT) real(kind=RP), allocatable :: Volume ! Volume of the element + real(kind=RP) :: mu_a real(kind=RP), allocatable :: x(:,:,:) ! Coordinates of the nodes ( xi , eta , X/Y ) real(kind=RP), allocatable :: Ja(:,:,:,:) ! Contravariant system metric matrix ( xi , eta , IROW , ICOL) real(kind=RP), allocatable :: Jac(:,:) ! Mapping jacobian ( xi , eta ) @@ -83,38 +85,43 @@ module QuadElementClass ! ********************************************************************************* ! type Edge_t - integer :: ID ! Edge ID - integer(kind=1) :: edgeType ! Edge Type: FACE_INTERIOR , or the marker value if boundary face - logical :: inverse ! Whether both edge projections are in the same or different direction - logical :: transform(QUADS_PER_EDGE) ! Whether the element contribution needs transformation to a higher degree - integer :: Nlow ! Lower polynomial degree - integer(kind=1), pointer :: edgeLocation(:) ! Edge location for the two (or one) sharing elements (ETOP,EBOTTOM,ELEFT,ERIGHT) - real(kind=RP) :: Area ! Area of the edge - real(kind=RP) :: invh ! Normal minimum distance approximation (inversed) - real(kind=RP), pointer :: n(:,:) ! Unitary normal: points from LEFT towards RIGHT, and outside the domain for bdryedges - real(kind=RP), pointer :: X(:,:) ! Coordinates: (X/Y, xi) - real(kind=RP), pointer :: dX(:,:) ! Tangent vector: (X/Y, xi) - real(kind=RP), pointer :: dS(:) ! Surface differential vector (X/Y, xi) - real(kind=RP), allocatable :: T_forward(:,:) ! Interpolation matrix from the low order to high order - real(kind=RP), allocatable :: T_backward(:,:) ! Interpolation matrix from the high order to low order - type(BoundaryData_t), allocatable :: storage(:) ! Solution interpolation to boundaries ( LEFT/RIGHT ) - type(Node_p) :: nodes(POINTS_PER_EDGE) ! Pointer to the two nodes - class(QuadElement_p), pointer :: quads(:) ! Pointers to the two (or one) shared quads - class(NodesAndWeights_t), pointer :: spA ! Pointer to the approximation nodal storage. In this case, is the largest from the two elements - class(NodesAndWeights_t), pointer :: spI ! Pointer to the integration nodal storage (if over-integration is active) + integer :: ID ! Edge ID + integer(kind=1) :: edgeType ! Edge Type: FACE_INTERIOR , or the marker value if boundary face + logical :: inverse ! Whether both edge projections are in the same or different direction + logical :: transform(QUADS_PER_EDGE) ! Whether the element contribution needs transformation to a higher degree + integer :: Nlow ! Lower polynomial degree + integer(kind=1), pointer :: edgeLocation(:) ! Edge location for the two (or one) sharing elements (ETOP,EBOTTOM,ELEFT,ERIGHT) + real(kind=RP) :: Area ! Area of the edge + real(kind=RP) :: mu_a + real(kind=RP) :: invh ! Normal minimum distance approximation (inversed) + real(kind=RP), pointer :: n(:,:) ! Unitary normal: points from LEFT towards RIGHT, and outside the domain for bdryedges + real(kind=RP), pointer :: X(:,:) ! Coordinates: (X/Y, xi) + real(kind=RP), pointer :: dX(:,:) ! Tangent vector: (X/Y, xi) + real(kind=RP), pointer :: dS(:) ! Surface differential vector (X/Y, xi) + real(kind=RP), allocatable :: T_forward(:,:) ! Interpolation matrix from the low order to high order + real(kind=RP), allocatable :: T_backward(:,:) ! Interpolation matrix from the high order to low order + type(BoundaryData_t), allocatable :: storage(:) ! Solution interpolation to boundaries ( LEFT/RIGHT ) + type(Node_p) :: nodes(POINTS_PER_EDGE) ! Pointer to the two nodes + class(QuadElement_p), pointer :: quads(:) ! Pointers to the two (or one) shared quads + class(NodesAndWeights_t), pointer :: spA ! Pointer to the approximation nodal storage. In this case, is the largest from the two elements + class(NodesAndWeights_t), pointer :: spI ! Pointer to the integration nodal storage (if over-integration is active) + procedure(ProjectSolutionFCN), nopass, pointer :: ProjectSolution => NULL() + procedure(ProjectFluxesFCN), nopass, pointer :: ProjectFluxes => NULL() + procedure(ProjectGradientFluxesFCN), nopass, pointer :: ProjectGradientFluxes => NULL() contains - procedure :: SetCurve => Edge_SetCurve ! Procedure that computes the coordinates, the tangent, and the normal. - procedure :: Invert => Edge_Invert ! Function to invert the edge orientation - procedure :: evaluateX => Edge_AnalyticalX ! Function to compute an edge point in a local coordinate "xi" - procedure :: evaluatedX => Edge_AnalyticaldX ! Function to compute an edge point in a local coordinate "xi" - procedure :: evaluatedS => Edge_AnalyticaldS - procedure :: getX => Edge_getX - procedure :: getdX => Edge_getdX - procedure :: getdS => Edge_getdS + procedure :: SetCurve => Edge_SetCurve ! Procedure that computes the coordinates, the tangent, and the normal. + procedure :: Invert => Edge_Invert ! Function to invert the edge orientation + procedure :: evaluateX => Edge_AnalyticalX ! Function to compute an edge point in a local coordinate "xi" + procedure :: evaluatedX => Edge_AnalyticaldX ! Function to compute an edge point in a local coordinate "xi" + procedure :: evaluatedS => Edge_AnalyticaldS + procedure :: getX => Edge_getX + procedure :: getdX => Edge_getdX + procedure :: getdS => Edge_getdS + procedure :: ComputeJumps => Edge_ComputeJumps end type Edge_t type, extends(Edge_t) :: StraightBdryEdge_t - integer(kind=1) :: BCWeakType + integer(kind=1) :: inviscidBCType logical :: associated = .false. procedure(RiemannSolverFunction), nopass, pointer :: RiemannSolver => NULL() real(kind=RP), pointer :: uB(:,:) => NULL() ! Solution at the boundary (used by the inviscid Riemann solver) @@ -124,10 +131,12 @@ module QuadElementClass real(kind=RP), pointer :: uSB(:,:) => NULL() ! Solution at the boundary (used by the solution Riemann solver) real(kind=RP), pointer :: gB(:,:,:) => NULL() ! Solution gradient at the boundary #endif + contains + procedure :: ComputeJumps => StraightBdryEdge_ComputeJumps end type StraightBdryEdge_t type, extends(Edge_t) :: CurvedBdryEdge_t - integer(kind=1) :: BCWeakType + integer(kind=1) :: inviscidBCType logical :: associated = .false. procedure(RiemannSolverFunction), nopass, pointer :: RiemannSolver => NULL() real(kind=RP), pointer :: uB(:,:) => NULL() ! Solution at the boundary (used by the Riemann solver) @@ -138,13 +147,14 @@ module QuadElementClass real(kind=RP), pointer :: gB(:,:,:) => NULL() ! Solution gradient at the boundary #endif contains - procedure :: SetCurve => CurvilinearEdge_SetCurve ! Procedure that computes the coordinates, the tangent, and the normal - procedure :: getX => Curvilinear_getX - procedure :: getdX => Curvilinear_getdX - procedure :: getdS => Curvilinear_getdS - procedure :: evaluateX => Curvilinear_InterpolantX - procedure :: evaluatedX => Curvilinear_InterpolantdX - procedure :: evaluatedS => Curvilinear_InterpolantdS + procedure :: SetCurve => CurvilinearEdge_SetCurve ! Procedure that computes the coordinates, the tangent, and the normal + procedure :: getX => Curvilinear_getX + procedure :: getdX => Curvilinear_getdX + procedure :: getdS => Curvilinear_getdS + procedure :: evaluateX => Curvilinear_InterpolantX + procedure :: evaluatedX => Curvilinear_InterpolantdX + procedure :: evaluatedS => Curvilinear_InterpolantdS + procedure :: ComputeJumps => CurvedBdryEdge_ComputeJumps end type CurvedBdryEdge_t ! @@ -159,6 +169,46 @@ module QuadElementClass procedure :: LinkWithElements => Edge_linkWithElements end type Edge_p ! +! + abstract interface + pure subroutine ProjectSolutionFCN( ed , QL , QR , dQL , dQR ) + use SMConstants + import Edge_t + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + end subroutine ProjectSolutionFCN + end interface + + abstract interface + pure subroutine ProjectFluxesFCN( ed , F , FL , FR ) + use SMConstants + import Edge_t + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + end subroutine ProjectFluxesFCN + end interface + + abstract interface + pure subroutine ProjectGradientFluxesFCN( ed , GL_edge , GR_edge , GL , GR) + use SMConstants + import Edge_t + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + end subroutine ProjectGradientFluxesFCN + end interface +! +! ! ======== contains ! ======== diff --git a/Solver/src/Mesh/QuadElement_EdgeProcedures.incf b/Solver/src/Mesh/QuadElement_EdgeProcedures.incf index 38db68a..bbcbeba 100644 --- a/Solver/src/Mesh/QuadElement_EdgeProcedures.incf +++ b/Solver/src/Mesh/QuadElement_EdgeProcedures.incf @@ -126,6 +126,41 @@ call TripleMatrixProduct( A = el2 % spA % Minv , B = self % f % T_forward , C = self % f % spA % M , val = self % f % T_backward , trB = .true. ) end if +! +! Point to the correct procedure +! ------------------------------ + if ( (.not. self % f % transform(LEFT)) .and. (.not. self % f % transform(RIGHT)) .and. (.not. self % f % inverse) ) then + self % f % ProjectSolution => Edge_ProjectSolutionType1 + self % f % ProjectFluxes => Edge_ProjectFluxesType1 + self % f % ProjectGradientFluxes => Edge_ProjectGradientFluxesType1 + + elseif ( (.not. self % f % transform(LEFT)) .and. (.not. self % f % transform(RIGHT)) .and. (self % f % inverse) ) then + self % f % ProjectSolution => Edge_ProjectSolutionType2 + self % f % ProjectFluxes => Edge_ProjectFluxesType2 + self % f % ProjectGradientFluxes => Edge_ProjectGradientFluxesType2 + + elseif ( ( self % f % transform(LEFT) ) .and. ( .not. self % f % inverse ) ) then + self % f % ProjectSolution => Edge_ProjectSolutionType3 + self % f % ProjectFluxes => Edge_ProjectFluxesType3 + self % f % ProjectGradientFluxes => Edge_ProjectGradientFluxesType3 + + elseif ( ( self % f % transform(LEFT) ) .and. ( self % f % inverse ) ) then + self % f % ProjectSolution => Edge_ProjectSolutionType4 + self % f % ProjectFluxes => Edge_ProjectFluxesType4 + self % f % ProjectGradientFluxes => Edge_ProjectGradientFluxesType4 + + elseif ( ( self % f % transform(RIGHT) ) .and. ( .not. self % f % inverse ) ) then + self % f % ProjectSolution => Edge_ProjectSolutionType5 + self % f % ProjectFluxes => Edge_ProjectFluxesType5 + self % f % ProjectGradientFluxes => Edge_ProjectGradientFluxesType5 + + elseif ( ( self % f % transform(RIGHT) ) .and. ( self % f % inverse ) ) then + self % f % ProjectSolution => Edge_ProjectSolutionType6 + self % f % ProjectFluxes => Edge_ProjectFluxesType6 + self % f % ProjectGradientFluxes => Edge_ProjectGradientFluxesType6 + + end if + elseif (present(elb) .and. (.not. present(el1)) .and. (.not. present(el2))) then ! Boundary edge @@ -235,6 +270,552 @@ end select end subroutine Edge_Invert +! +!////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +! +! PROJECTION PROCEDURES +! --------------------- +!////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +! + pure subroutine Edge_ProjectSolutionType1( ed , QL , QR , dQL , dQR) + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + + QL = ed % storage(LEFT) % Q + QR = ed % storage(RIGHT) % Q + + dQL = ed % storage(LEFT) % dQ + dQR = ed % storage(RIGHT) % dQ + + end subroutine Edge_ProjectSolutionType1 + + pure subroutine Edge_ProjectSolutionType2( ed , QL , QR , dQL , dQR ) + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + + QL = ed % storage(LEFT) % Q + QR = ed % storage(RIGHT) % Q( ed % spA % N : 0 : -1 , 1 : NCONS ) + +#ifdef NAVIER_STOKES + dQL = ed % storage(LEFT) % dQ + dQR = ed % storage(RIGHT) % dQ ( ed % spA % N : 0 : -1 , 1 : NDIM , 1 : NCONS ) +#else + dQL = 0.0_RP + dQR = 0.0_RP +#endif + + end subroutine Edge_ProjectSolutionType2 + + pure subroutine Edge_ProjectSolutionType3( ed , QL , QR , dQL , dQR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) +! +! --------------- +! Local variables +! --------------- +! + integer :: eq + + call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) + QR = ed % storage(RIGHT) % Q + +#ifdef NAVIER_STOKES + do eq = 1 , NCONS + call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % dQ(:,:,eq) , dQL(:,:,eq) ) + end do + dQR = ed % storage(RIGHT) % dQ +#else + dQL = 0.0_RP + dQR = 0.0_RP +#endif + + end subroutine Edge_ProjectSolutionType3 + + pure subroutine Edge_ProjectSolutionType4( ed , QL , QR , dQL , dQR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) +! +! --------------- +! Local variables +! --------------- +! + integer :: eq + + call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % Q , QL) + QR = ed % storage(RIGHT) % Q( ed % spA % N : 0 : -1 , 1 : NCONS ) + +#ifdef NAVIER_STOKES + do eq = 1 , NCONS + call Mat_x_Mat( ed % T_forward , ed % storage(LEFT) % dQ(:,:,eq) , dQL(:,:,eq) ) + end do + dQR = ed % storage(RIGHT) % dQ( ed % spA % N : 0 : -1 , 1 : NDIM , 1 : NCONS ) +#else + dQL = 0.0_RP + dQR = 0.0_RP +#endif + + end subroutine Edge_ProjectSolutionType4 + + pure subroutine Edge_ProjectSolutionType5( ed , QL , QR , dQL , dQR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) +! +! --------------- +! Local variables +! --------------- +! + integer :: eq + + QL = ed % storage(LEFT) % Q + call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) + +#ifdef NAVIER_STOKES + dQL = ed % storage(LEFT) % dQ + do eq = 1 , NCONS + call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % dQ(:,:,eq) , dQR(:,:,eq) ) + end do +#else + dQL = 0.0_RP + dQR = 0.0_RP +#endif + + end subroutine Edge_ProjectSolutionType5 + + pure subroutine Edge_ProjectSolutionType6( ed , QL , QR , dQL , dQR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(out) :: QL( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: QR( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: dQL( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) + real(kind=RP), intent(out) :: dQR( 0 : ed % spA % N , 1 : NDIM , 1 : NCONS ) +! +! --------------- +! Local variables +! --------------- +! + integer :: eq + + QL = ed % storage(LEFT) % Q + call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % Q , QR) + QR = QR( ed % spA % N : 0 : -1 , 1 : NCONS ) + +#ifdef NAVIER_STOKES + dQL = ed % storage(LEFT) % dQ + do eq = 1 , NCONS + call Mat_x_Mat( ed % T_forward , ed % storage(RIGHT) % dQ(:,:,eq) , dQR(:,:,eq) ) + end do + dQR = ed % storage(RIGHT) % dQ( ed % spA % N : 0 : -1 , 1 : NDIM , 1 : NCONS ) +#else + dQL = 0.0_RP + dQR = 0.0_RP +#endif + + end subroutine Edge_ProjectSolutionType6 + + pure subroutine Edge_ProjectFluxesType1( ed , F , FL , FR) + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + + FL = F + FR = -F + + end subroutine Edge_ProjectFluxesType1 + + pure subroutine Edge_ProjectFluxesType2( ed , F , FL , FR ) + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + + FL = F + FR = -F( ed % spA % N : 0 : -1 , 1 : NCONS ) + + end subroutine Edge_ProjectFluxesType2 + + pure subroutine Edge_ProjectFluxesType3( ed , F , FL , FR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + + call Mat_x_Mat( ed % T_backward , F , FL) + FR = -F + + end subroutine Edge_ProjectFluxesType3 + + pure subroutine Edge_ProjectFluxesType4( ed , F , FL , FR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + + call Mat_x_Mat( ed % T_backward , F , FL ) + FR = -F( ed % spA % N : 0 : -1 , 1 : NCONS ) + + end subroutine Edge_ProjectFluxesType4 + + pure subroutine Edge_ProjectFluxesType5( ed , F , FL , FR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + + FL = F + call Mat_x_Mat( ed % T_backward , -F , FR ) + + end subroutine Edge_ProjectFluxesType5 + + pure subroutine Edge_ProjectFluxesType6( ed , F , FL , FR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: F ( 0 : ed % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS ) + real(kind=RP), intent(out) :: FR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS ) + + FL = F + call Mat_x_Mat( ed % T_backward , F , FR ) + FR = -FR( ed % Nlow : 0 : -1 , 1 : NCONS ) + + end subroutine Edge_ProjectFluxesType6 + + pure subroutine Edge_ProjectGradientFluxesType1( ed , GL_edge , GR_edge , GL , GR) + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + + GL = GL_edge + GR = GR_edge + + end subroutine Edge_ProjectGradientFluxesType1 + + pure subroutine Edge_ProjectGradientFluxesType2( ed , GL_edge , GR_edge , GL , GR) + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + + GL = GL_edge + GR = GR_edge( ed % spA % N : 0 : -1 , 1 : NCONS , 1 : NDIM ) + + end subroutine Edge_ProjectGradientFluxesType2 + + pure subroutine Edge_ProjectGradientFluxesType3( ed , GL_edge , GR_edge , GL , GR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + + call Mat_x_Mat( ed % T_backward , GL_edge(:,:,IX) , GL(:,:,IX) ) + call Mat_x_Mat( ed % T_backward , GL_edge(:,:,IY) , GL(:,:,IY) ) + GR = GR_edge + + end subroutine Edge_ProjectGradientFluxesType3 + + pure subroutine Edge_ProjectGradientFluxesType4( ed , GL_edge , GR_edge , GL , GR) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + + call Mat_x_Mat( ed % T_backward , GL_edge(:,:,IX) , GL(:,:,IX) ) + call Mat_x_Mat( ed % T_backward , GL_edge(:,:,IY) , GL(:,:,IY) ) + GR = GR_edge( ed % spA % N : 0 : -1 , 1 : NCONS , 1 : NDIM) + + end subroutine Edge_ProjectGradientFluxesType4 + + pure subroutine Edge_ProjectGradientFluxesType5( ed , GL_edge , GR_edge , GL , GR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + + GL = GL_edge + call Mat_x_Mat( ed % T_backward , GR_edge(:,:,IX) , GR(:,:,IX) ) + call Mat_x_Mat( ed % T_backward , GR_edge(:,:,IY) , GR(:,:,IY) ) + + end subroutine Edge_ProjectGradientFluxesType5 + + pure subroutine Edge_ProjectGradientFluxesType6( ed , GL_edge , GR_edge , GL , GR ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + real(kind=RP), intent(in) :: GL_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(in) :: GR_edge ( 0 : ed % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GL( 0 : ed % storage(LEFT) % spA % N , 1 : NCONS , 1 : NDIM ) + real(kind=RP), intent(out) :: GR( 0 : ed % storage(RIGHT) % spA % N , 1 : NCONS , 1 : NDIM ) + + GL = GL_edge + call Mat_x_Mat( ed % T_backward , GR_edge(:,:,IX) , GR(:,:,IX) ) + call Mat_x_Mat( ed % T_backward , GR_edge(:,:,IY) , GR(:,:,IY) ) + + GR = GR( ed % Nlow : 0 : -1 , 1 : NCONS , 1 : NDIM ) + + end subroutine Edge_ProjectGradientFluxesType6 + + pure function Edge_ComputeJumps ( ed , var ) result ( jumps ) + use MatrixOperations + implicit none + class(Edge_t), intent(in) :: ed + integer , intent(in) :: var + real(kind=RP) :: jumps +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: QL( 0 : ed % spA % N ) + real(kind=RP) :: QR( 0 : ed % spA % N ) + real(kind=RP) :: nodalJumps ( 0 : ed % spA % N ) + + if ( ed % transform(LEFT) .and. ed % inverse ) then +! +! --------------------------------------------------------- +!> First case: LEFT needs p-Adaption, and RIGHT is reversed. +! --------------------------------------------------------- +! +! Transform the LEFT edge +! ----------------------- + QL = MatrixTimesVector_F( ed % T_forward , ed % storage(LEFT) % Q(:,var) , ed % spA % N + 1 ) +! +! Invert the RIGHT edge +! --------------------- + QR = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , var ) + + elseif ( ed % transform(LEFT) ) then +! +! ---------------------------------- +!> Second case: LEFT needs p-Adaption. +! ---------------------------------- +! +! Transform the LEFT +! ----------------------- + QL = MatrixTimesVector_F( ed % T_forward , ed % storage(LEFT) % Q(:,var) , ed % spA % N + 1 ) +! +! Get the RIGHT edge state +! ------------------------ + QR = ed % storage(RIGHT) % Q(:,var) + + elseif ( ed % transform(RIGHT) .and. ed % inverse ) then +! +! --------------------------------------------------------- +!> Third case: RIGHT needs p-Adaption, and RIGHT is reversed. +! --------------------------------------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(LEFT) % Q(:,var) +! +! Transform the RIGHT edge +! ------------------------ + QR = MatrixTimesVector_F( ed % T_forward , ed % storage(RIGHT) % Q(:,var) , ed % spA % N + 1 ) +! +! Invert the RIGHT edge +! --------------------- + QR(0:ed % spA % N) = QR(ed % spA % N : 0 : -1) + + elseif ( ed % transform(RIGHT) ) then +! +! ----------------------------------- +!> Fourth case: RIGHT needs p-Adaption. +! ----------------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(LEFT) % Q(:,var) +! +! Transform the RIGHT edge +! ------------------------ + QR = MatrixTimesVector_F( ed % T_forward , ed % storage(RIGHT) % Q(:,var) , ed % spA % N + 1 ) + + elseif ( ed % inverse ) then +! +! ----------------------------- +!> Fifth case: RIGHT is reversed. +! ----------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(LEFT) % Q(:,var) +! +! Invert the RIGHT edge +! --------------------- + QR = ed % storage(RIGHT) % Q(ed % spA % N : 0 : -1 , var) + + else +! +! ----------------------------------------------------------------------- +!> Sixth case: Default case: neither p-Adaption nor inversion are required. +! ----------------------------------------------------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(LEFT) % Q(:,var) +! +! Get the RIGHT edge +! ------------------ + QR = ed % storage(RIGHT) % Q(:,var) + + end if +! +! Compute the Jumps +! ----------------- + nodalJumps = 2.0_RP * ( QL - QR ) / ( QL + QR ) + jumps = sum( ed % spA % w * nodalJumps * nodalJumps * ed % dS(0)) + + end function Edge_ComputeJumps + + pure function StraightBdryEdge_ComputeJumps ( ed , var ) result ( jumps ) + use MatrixOperations + implicit none + class(StraightBdryEdge_t), intent(in) :: ed + integer , intent(in) :: var + real(kind=RP) :: jumps +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: QL( 0 : ed % spA % N ) + real(kind=RP) :: QR( 0 : ed % spA % N ) + real(kind=RP) :: nodalJumps ( 0 : ed % spA % N ) + + if ( ed % inverse ) then +! +! ----------------------------- +!> First case: RIGHT is reversed. +! ----------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(1) % Q(:,var) +! +! Invert the RIGHT edge +! --------------------- + QR(0:ed % spA % N) = ed % uSB(ed % spA % N : 0 : -1 , var ) + + else +! +! ----------------------------------------------------------------------- +!> Second case: Default case: neither p-Adaption nor inversion are required. +! ----------------------------------------------------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(1) % Q(:,var) +! +! Get the RIGHT edge +! ------------------ + QR = ed % uSB(:,var) + + end if +! +! Compute the Jumps +! ----------------- + nodalJumps = 2.0_RP * ( QL - QR ) / ( QL + QR ) + jumps = sum( ed % spA % w * nodalJumps * nodalJumps * ed % dS(0)) + + end function StraightBdryEdge_ComputeJumps + + pure function CurvedBdryEdge_ComputeJumps ( ed , var ) result ( jumps ) + use MatrixOperations + implicit none + class(CurvedBdryEdge_t), intent(in) :: ed + integer , intent(in) :: var + real(kind=RP) :: jumps +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: QL( 0 : ed % spA % N ) + real(kind=RP) :: QR( 0 : ed % spA % N ) + real(kind=RP) :: nodalJumps( 0 : ed % spA % N ) + + if ( ed % inverse ) then +! +! ----------------------------- +!> First case: RIGHT is reversed. +! ----------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(1) % Q(:,var) +! +! Invert the RIGHT edge +! --------------------- + QR(0:ed % spA % N) = ed % uSB(ed % spA % N : 0 : -1 , var ) + + else +! +! ----------------------------------------------------------------------- +!> Second case: Default case: neither p-Adaption nor inversion are required. +! ----------------------------------------------------------------------- +! +! Get the LEFT edge +! ----------------- + QL = ed % storage(1) % Q(:,var) +! +! Get the RIGHT edge +! ------------------ + QR = ed % uSB(:,var) + + end if +! +! Compute the Jumps +! ----------------- + nodalJumps = 2.0_RP * ( QL - QR ) / ( QL + QR ) + jumps = sum( ed % spA % w * nodalJumps * nodalJumps * ed % dS(0)) + + end function CurvedBdryEdge_ComputeJumps subroutine BoundaryData_Initialize( self , spA ) implicit none diff --git a/Solver/src/Mesh/QuadElement_QuadProcedures.incf b/Solver/src/Mesh/QuadElement_QuadProcedures.incf index dc5dab7..d335384 100644 --- a/Solver/src/Mesh/QuadElement_QuadProcedures.incf +++ b/Solver/src/Mesh/QuadElement_QuadProcedures.incf @@ -18,7 +18,7 @@ end subroutine QuadElement_SetStorage #ifdef NAVIER_STOKES - subroutine QuadElement_ComputeInteriorGradient( self ) + function QuadElement_ComputeInteriorGradient( self ) result (dQ) use MatrixOperations ! ! ********************************************************************** @@ -30,6 +30,7 @@ ! implicit none class(QuadElement_t) :: self + real(kind=RP) :: dQ(0:self % spA % N , 0:self % spA % N , 1:NDIM , 1:NCONS) ! ------------------------------------------------------------- real(kind=RP) :: dxiQ (0:self % spA % N,0:self % spA % N, 1:NCONS) real(kind=RP) :: detaQ(0:self % spA % N,0:self % spA % N, 1:NCONS) @@ -44,14 +45,14 @@ ! ! x-direction gradient ! -------------------- - self % dQ(0:N,0:N,IX,var) = (dxiQ(0:N,0:N,var) * self % Ja(0:N,0:N,1,1) + detaQ(0:N,0:N,var) * self % Ja(0:N,0:N,1,2)) / self % jac(0:N,0:N) + dQ(0:N,0:N,IX,var) = (dxiQ(0:N,0:N,var) * self % Ja(0:N,0:N,1,1) + detaQ(0:N,0:N,var) * self % Ja(0:N,0:N,1,2)) / self % jac(0:N,0:N) ! ! y-direction gradient ! -------------------- - self % dQ(0:N,0:N,IY,var) = (dxiQ(0:N,0:N,var) * self % Ja(0:N,0:N,2,1) + detaQ(0:N,0:N,var) * self % Ja(0:N,0:N,2,2)) / self % jac(0:N,0:N) + dQ(0:N,0:N,IY,var) = (dxiQ(0:N,0:N,var) * self % Ja(0:N,0:N,2,1) + detaQ(0:N,0:N,var) * self % Ja(0:N,0:N,2,2)) / self % jac(0:N,0:N) end do end associate - end subroutine QuadElement_ComputeInteriorGradient + end function QuadElement_ComputeInteriorGradient #endif diff --git a/Solver/src/Mesh/QuadMesh_Zones.incf b/Solver/src/Mesh/QuadMesh_Zones.incf index 483cca0..8e07a35 100644 --- a/Solver/src/Mesh/QuadMesh_Zones.incf +++ b/Solver/src/Mesh/QuadMesh_Zones.incf @@ -84,8 +84,8 @@ edge1 % inverse = areInverted edge2 % associated = .true. edge1 % associated = .true. - edge1 % BCWeakType = WEAK_RIEMANN - edge2 % BCWeakType = WEAK_RIEMANN + edge1 % inviscidBCType = WEAK_RIEMANN + edge2 % inviscidBCType = WEAK_RIEMANN edge1 % uB(0:,1:) => edge2 % storage(1) % Q(0:,1:) edge2 % uB(0:,1:) => edge1 % storage(1) % Q(0:,1:) diff --git a/Solver/src/Physics/PhysicsNS.f90 b/Solver/src/Physics/PhysicsNS.f90 index 5bc2367..a98cf75 100755 --- a/Solver/src/Physics/PhysicsNS.f90 +++ b/Solver/src/Physics/PhysicsNS.f90 @@ -59,12 +59,13 @@ module PhysicsNS end type Dimensionless_t abstract interface - pure function RiemannSolverFunction( QL , QR , n ) result ( val ) + pure function RiemannSolverFunction( N , QL , QR , normal ) result ( Fstar ) use SMConstants - real(kind=RP), dimension(NCONS), intent(in) :: QL - real(kind=RP), dimension(NCONS), intent(in) :: QR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: val + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS), intent(in) :: QL + real(kind=RP), dimension(0:N,NCONS), intent(in) :: QR + real(kind=RP), dimension(NDIM,0:N) , intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar end function RiemannSolverFunction end interface @@ -227,31 +228,34 @@ end function getTemperatureGradient2D end interface interface - module pure function ExactRiemannSolver(qL , qR , n) result (Fstar) + module pure function ExactRiemannSolver( N , qL , qR , normal) result (Fstar) use MatrixOperations implicit none - real(kind=RP), dimension(NCONS), intent(in) :: qL - real(kind=RP), dimension(NCONS), intent(in) :: qR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: Fstar + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qL + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qR + real(kind=RP), dimension(NDIM,0:N ), intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar end function ExactRiemannSolver - module pure function RoeFlux(qL, qR , n) result(Fstar) + module pure function RoeFlux(N , qL , qR , normal) result(Fstar) use MatrixOperations implicit none - real(kind=RP), dimension(NCONS), intent(in) :: qL - real(kind=RP), dimension(NCONS), intent(in) :: qR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: Fstar + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qL + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qR + real(kind=RP), dimension(NDIM,0:N ), intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar end function RoeFlux - module pure function HLLFlux(qL , qR , n) result(Fstar) + module pure function HLLFlux(N , qL , qR , normal) result(Fstar) use MatrixOperations implicit none - real(kind=RP), dimension(NCONS), intent(in) :: qL - real(kind=RP), dimension(NCONS), intent(in) :: qR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: Fstar + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qL + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qR + real(kind=RP), dimension(NDIM,0:N ), intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar end function HLLFlux end interface @@ -423,8 +427,17 @@ subroutine InitializePhysics() dimensionless % sqrtGammaMach = sqrt(thermodynamics % gamma) * dimensionless % Mach dimensionless % gammaMach2 = thermodynamics % gamma * dimensionless % Mach * dimensionless % Mach dimensionless % invSqrtGammaMach = 1.0_RP / (sqrt(thermodynamics % gamma) * dimensionless % Mach) - dimensionless % mu = 1.0_RP / Setup % reynolds_number - dimensionless % kappa = thermodynamics % gogm1 / (Setup % prandtl_number * Setup % reynolds_number) + + if ( setup % reynolds_number .eq. 0.0_RP ) then + dimensionless % mu = 0.0_RP + dimensionless % kappa = 0.0_RP + + else + dimensionless % mu = 1.0_RP / Setup % reynolds_number + dimensionless % kappa = thermodynamics % gogm1 / (Setup % prandtl_number * Setup % reynolds_number) + + end if + dimensionless % cp = thermodynamics % gogm1 dimensionless % cv = thermodynamics % invgm1 dimensionless % Re = Setup % reynolds_number diff --git a/Solver/src/Physics/RiemannSolvers.f90 b/Solver/src/Physics/RiemannSolvers.f90 index b57914f..c88c36b 100644 --- a/Solver/src/Physics/RiemannSolvers.f90 +++ b/Solver/src/Physics/RiemannSolvers.f90 @@ -9,31 +9,40 @@ ! Riemann solvers ! **************************************************** ! - module pure function ExactRiemannSolver(qL , qR , n) result (Fstar) + module pure function ExactRiemannSolver(N , qL , qR , normal) result (Fstar) use MatrixOperations implicit none - real(kind=RP), dimension(NCONS), intent(in) :: qL - real(kind=RP), dimension(NCONS), intent(in) :: qR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: Fstar -! --------------------------------------------------------------- + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS) , intent(in) :: qL + real(kind=RP), dimension(0:N,NCONS) , intent(in) :: qR + real(kind=RP), dimension(NDIM,0:N ) , intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar +! +! --------------- +! Local variables +! --------------- +! real(kind=RP), dimension(NCONS) :: qnL , qnR real(kind=RP), dimension(NPRIM) :: wL , wR real(kind=RP) :: pstar , ustar real(kind=RP) :: rhostar , uFan , pFan real(kind=RP) :: Ft + integer :: i + + do i = 0 , N +! ! 0/ Gather variables ! ---------------- - qnL(IRHO) = qL(IRHO) - qnL(IRHOU) = qL(IRHOU) * n(IX) + qL(IRHOV) * n(IY) - qnL(IRHOV) = -qL(IRHOU) * n(IY) + qL(IRHOV) * n(IX) - qnL(IRHOE) = qL(IRHOE) + qnL(IRHO) = qL(i,IRHO) + qnL(IRHOU) = qL(i,IRHOU) * normal(IX,i) + qL(i,IRHOV) * normal(IY,i) + qnL(IRHOV) = -qL(i,IRHOU) * normal(IY,i) + qL(i,IRHOV) * normal(IX,i) + qnL(IRHOE) = qL(i,IRHOE) - qnR(IRHO) = qR(IRHO) - qnR(IRHOU) = qR(IRHOU) * n(IX) + qR(IRHOV) * n(IY) - qnR(IRHOV) = -qR(IRHOU) * n(IY) + qR(IRHOV) * n(IX) - qnR(IRHOE) = qR(IRHOE) + qnR(IRHO) = qR(i,IRHO) + qnR(IRHOU) = qR(i,IRHOU) * normal(IX,i) + qR(i,IRHOV) * normal(IY,i) + qnR(IRHOV) = -qR(i,IRHOU) * normal(IY,i) + qR(i,IRHOV) * normal(IX,i) + qnR(IRHOE) = qR(i,IRHOE) associate( gamma => Thermodynamics % gamma , gm1 => Thermodynamics % gm1 ) @@ -49,29 +58,29 @@ module pure function ExactRiemannSolver(qL , qR , n) result (Fstar) wR(IP) = gm1 * (qnR(IRHOE) - 0.5_RP * (qnR(IRHOU) * wR(IU) + qnR(IRHOV) * wR(IV) ) ) wR(IA) = sqrt(gamma * wR(IP) / wR(IRHO)) end associate - +! ! 1/ Compute the star region ! ----------------------- call ExactRiemann_ComputePStar (wL , wR , pstar , ustar ) - +! ! 2/ Check to which region belongs the solution ! ------------------------------------------ associate ( cv => Dimensionless % cv , cp => Dimensionless % cp , gamma => Thermodynamics % gamma ) if ( ustar .ge. 0.0_RP ) then if ( ( pstar .le. wL(IP) ) .and. ( wL(IU) .ge. wL(IA) ) ) then - Fstar(IRHO) = qL(IRHOU) - Fstar(IRHOU) = qL(IRHOU) * wL(IU) + wL(IP) - Fstar(IRHOV) = qL(IRHOV) * wL(IU) - Fstar(IRHOE) = wL(IU)*( cp * wL(IP) + 0.5_RP * (qL(IRHOU)*wL(IU) + qL(IRHOV)*wL(IV)) ) + Fstar(i,IRHO) = qnL(IRHOU) + Fstar(i,IRHOU) = qnL(IRHOU) * wL(IU) + wL(IP) + Fstar(i,IRHOV) = qnL(IRHOV) * wL(IU) + Fstar(i,IRHOE) = wL(IU)*( cp * wL(IP) + 0.5_RP * (qnL(IRHOU)*wL(IU) + qnL(IRHOV)*wL(IV)) ) elseif ( ( pstar .le. wL(IP) ) .and. ( wL(IU) .lt. wL(IA)) .and. (ustar .lt. wL(IA) * (pstar/wL(IP))**(0.5_RP / cp ) ) ) then rhostar = wL(IRHO) * ( pstar / wL(IP) ) ** ( 1.0_RP / gamma ) - Fstar(IRHO) = rhostar * ustar - Fstar(IRHOU) = rhostar * ustar * ustar + pstar - Fstar(IRHOV) = Fstar(IRHO) * wL(IV) - Fstar(IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wL(IV) * wL(IV)) ) + Fstar(i,IRHO) = rhostar * ustar + Fstar(i,IRHOU) = rhostar * ustar * ustar + pstar + Fstar(i,IRHOV) = Fstar(i,IRHO) * wL(IV) + Fstar(i,IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wL(IV) * wL(IV)) ) elseif ( ( pstar .le. wL(IP) ) .and. ( wL(IU) .lt. wL(IA) ) .and. ( ustar .ge. wL(IA) * (pstar / wL(IP))**(0.5_RP / cp))) then @@ -79,64 +88,64 @@ module pure function ExactRiemannSolver(qL , qR , n) result (Fstar) rhostar = wL(IRHO) * (uFan / wL(IA)) ** (2.0_RP * cv) pFan = wL(IP) * (rhostar / wL(IRHO)) ** (gamma) - Fstar(IRHO) = rhostar * uFan - Fstar(IRHOU) = rhostar * uFan * uFan + pFan - Fstar(IRHOV) = rhostar * uFan * wL(IV) - Fstar(IRHOE) = uFan * ( cp * pFan + 0.5_RP * rhostar * ( uFan*uFan + wL(IV)*wL(IV) ) ) + Fstar(i,IRHO) = rhostar * uFan + Fstar(i,IRHOU) = rhostar * uFan * uFan + pFan + Fstar(i,IRHOV) = rhostar * uFan * wL(IV) + Fstar(i,IRHOE) = uFan * ( cp * pFan + 0.5_RP * rhostar * ( uFan*uFan + wL(IV)*wL(IV) ) ) elseif ( ( pstar .gt. wL(IP) ) .and. ( wL(IU) .ge. wL(IA) * sqrt( (gamma+1.0_RP)/(2.0_RP * gamma)*pstar/wL(IP) + 0.5_RP/cp) ) ) then - Fstar(IRHO) = qL(IRHOU) - Fstar(IRHOU) = qL(IRHOU) * wL(IU) + wL(IP) - Fstar(IRHOV) = qL(IRHOV) * wL(IU) - Fstar(IRHOE) = wL(IU)*( cp * wL(IP) + 0.5_RP * (qL(IRHOU)*wL(IU) + qL(IRHOV)*wL(IV)) ) + Fstar(i,IRHO) = qnL(IRHOU) + Fstar(i,IRHOU) = qnL(IRHOU) * wL(IU) + wL(IP) + Fstar(i,IRHOV) = qnL(IRHOV) * wL(IU) + Fstar(i,IRHOE) = wL(IU)*( cp * wL(IP) + 0.5_RP * (qnL(IRHOU)*wL(IU) + qnL(IRHOV)*wL(IV)) ) elseif ( ( pstar .gt. wL(IP) ) .and. ( wL(IU) .lt. wL(IA) * sqrt( (gamma+1.0_RP)/(2.0_RP * gamma)*pstar/wL(IP) + 0.5_RP/cp) ) ) then rhostar = wL(IRHO) * ( (pstar/wL(IP) + (gamma-1.0_RP)/(gamma+1.0_RP))/(pstar/wL(IP)*(gamma-1.0_RP)/(gamma+1.0_RP) + 1.0_RP)) - Fstar(IRHO) = rhostar * ustar - Fstar(IRHOU) = rhostar * ustar * ustar + pstar - Fstar(IRHOV) = rhostar * ustar * wL(IV) - Fstar(IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wL(IV)*wL(IV) ) ) + Fstar(i,IRHO) = rhostar * ustar + Fstar(i,IRHOU) = rhostar * ustar * ustar + pstar + Fstar(i,IRHOV) = rhostar * ustar * wL(IV) + Fstar(i,IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wL(IV)*wL(IV) ) ) end if else ! ( ustar .lt. 0.0_RP ) if ( ( pstar .le. wR(IP) ) .and. ( wR(IU) + wR(IA) .le. 0.0_RP ) ) then - Fstar(IRHO) = qR(IRHOU) - Fstar(IRHOU) = qR(IRHOU) * wR(IU) + wR(IP) - Fstar(IRHOV) = qR(IRHOV) * wR(IU) - Fstar(IRHOE) = wR(IU)*( cp * wR(IP) + 0.5_RP * (qR(IRHOU)*wR(IU) + qR(IRHOV)*wR(IV)) ) + Fstar(i,IRHO) = qnR(IRHOU) + Fstar(i,IRHOU) = qnR(IRHOU) * wR(IU) + wR(IP) + Fstar(i,IRHOV) = qnR(IRHOV) * wR(IU) + Fstar(i,IRHOE) = wR(IU)*( cp * wR(IP) + 0.5_RP * (qnR(IRHOU)*wR(IU) + qnR(IRHOV)*wR(IV)) ) elseif ( ( pstar .le. wR(IP) ) .and. (wR(IU) + wR(IA) .ge. 0.0_RP) .and. (ustar + wR(IA)*(pstar/wR(IP))**(0.5_RP / cp) .ge. 0.0_RP) ) then rhostar = wR(IRHO) * ( pstar / wR(IP) ) ** ( 1.0_RP / gamma ) - Fstar(IRHO) = rhostar * ustar - Fstar(IRHOU) = rhostar * ustar * ustar + pstar - Fstar(IRHOV) = Fstar(IRHO) * wR(IV) - Fstar(IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wR(IV) * wR(IV)) ) + Fstar(i,IRHO) = rhostar * ustar + Fstar(i,IRHOU) = rhostar * ustar * ustar + pstar + Fstar(i,IRHOV) = Fstar(i,IRHO) * wR(IV) + Fstar(i,IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wR(IV) * wR(IV)) ) elseif ( ( pstar .le. wR(IP) ) .and. (wR(IU) + wR(IA) .gt. 0.0_RP) .and. (ustar + wR(IA) * (pstar/wR(IP))**(0.5_RP / cp) .lt. 0.0_RP ) ) then uFan = -2.0_RP * wR(IA) / ( gamma + 1.0_RP) + (gamma-1.0_RP)/(gamma+1.0_RP) * wR(IU) rhostar = wR(IRHO) * (-uFan / wR(IA)) ** (2.0_RP * cv) pFan = wR(IP) * (rhostar / wR(IRHO)) ** (gamma) - Fstar(IRHO) = rhostar * uFan - Fstar(IRHOU) = rhostar * uFan * uFan + pFan - Fstar(IRHOV) = rhostar * uFan * wR(IV) - Fstar(IRHOE) = uFan * ( cp * pFan + 0.5_RP * rhostar * ( uFan*uFan + wR(IV)*wR(IV) ) ) + Fstar(i,IRHO) = rhostar * uFan + Fstar(i,IRHOU) = rhostar * uFan * uFan + pFan + Fstar(i,IRHOV) = rhostar * uFan * wR(IV) + Fstar(i,IRHOE) = uFan * ( cp * pFan + 0.5_RP * rhostar * ( uFan*uFan + wR(IV)*wR(IV) ) ) elseif ( (pstar .gt. wR(IP)) .and. (wR(IU) + wR(IA)*( 0.5_RP*(gamma+1.0_RP)/gamma*pstar/wR(IP) + 0.5_RP/cp ) .le. 0.0_RP) ) then - Fstar(IRHO) = qR(IRHOU) - Fstar(IRHOU) = qR(IRHOU) * wR(IU) + wR(IP) - Fstar(IRHOV) = qR(IRHOV) * wR(IU) - Fstar(IRHOE) = wR(IU)*( cp * wR(IP) + 0.5_RP * (qR(IRHOU)*wR(IU) + qR(IRHOV)*wR(IV)) ) + Fstar(i,IRHO) = qnR(IRHOU) + Fstar(i,IRHOU) = qnR(IRHOU) * wR(IU) + wR(IP) + Fstar(i,IRHOV) = qnR(IRHOV) * wR(IU) + Fstar(i,IRHOE) = wR(IU)*( cp * wR(IP) + 0.5_RP * (qnR(IRHOU)*wR(IU) + qnR(IRHOV)*wR(IV)) ) elseif ( (pstar .gt. wR(IP)) .and. (wR(IU) + wR(IA)*( (gamma+1.0_RP)/(2.0_RP*gamma)*pstar/wR(IP) + 0.5_RP/cp) .gt. 0.0_RP) ) then rhostar = wR(IRHO) * ( (pstar/wR(IP) + (gamma-1.0_RP)/(gamma+1.0_RP))/(pstar/wR(IP)*(gamma-1.0_RP)/(gamma+1.0_RP) + 1.0_RP)) - Fstar(IRHO) = rhostar * ustar - Fstar(IRHOU) = rhostar * ustar * ustar + pstar - Fstar(IRHOV) = rhostar * ustar * wR(IV) - Fstar(IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wR(IV)*wR(IV) ) ) + Fstar(i,IRHO) = rhostar * ustar + Fstar(i,IRHOU) = rhostar * ustar * ustar + pstar + Fstar(i,IRHOV) = rhostar * ustar * wR(IV) + Fstar(i,IRHOE) = ustar * ( cp * pstar + 0.5_RP * rhostar * (ustar * ustar + wR(IV)*wR(IV) ) ) end if end if @@ -144,26 +153,33 @@ module pure function ExactRiemannSolver(qL , qR , n) result (Fstar) ! 3/ Return to the 3D Space ! ---------------------- - Ft = Fstar(IRHOU) * n(IX) - Fstar(IRHOV) * n(IY) - Fstar(IRHOV) = Fstar(IRHOU) * n(IY) + Fstar(IRHOV) * n(IX) - Fstar(IRHOU) = Ft + Ft = Fstar(i,IRHOU) * normal(IX,i) - Fstar(i,IRHOV) * normal(IY,i) + Fstar(i,IRHOV) = Fstar(i,IRHOU) * normal(IY,i) + Fstar(i,IRHOV) * normal(IX,i) + Fstar(i,IRHOU) = Ft + + end do #ifdef _DIMENSIONLESS_TAU - Fstar = Fstar * dimensionless % invSqrtGammaMach + Fstar = Fstar * dimensionless % invSqrtGammaMach #endif end function ExactRiemannSolver - module pure function RoeFlux(qL, qR , n) result(Fstar) + module pure function RoeFlux(N , qL, qR , normal) result(Fstar) use MatrixOperations implicit none - real(kind=RP), dimension(NCONS), intent(in) :: qL - real(kind=RP), dimension(NCONS), intent(in) :: qR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: Fstar -! --------------------------------------------------------------- - real(kind=RP) :: sqrtRhoL , invRhoL , rhoL , uL , vL , HL , TL , pL - real(kind=RP) :: sqrtRhoR , invRhoR , rhoR , uR , vR , HR , TR , pR + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qL + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qR + real(kind=RP), dimension(NDIM,0:N ), intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: sqrtRhoL , invRhoL , rhoL , uL , vL , HL , pL + real(kind=RP) :: sqrtRhoR , invRhoR , rhoR , uR , vR , HR , pR real(kind=RP) :: invrho , u , v , H , a real(kind=RP) :: dq(NCONS) real(kind=RP) :: lambda(NCONS) @@ -172,26 +188,27 @@ module pure function RoeFlux(qL, qR , n) result(Fstar) real(kind=RP) :: Ft integer :: negativeWaves integer :: wave + integer :: i + do i = 0 , N +! ! 0/ Gather variables ! ---------------- - rhoL = qL(IRHO) + rhoL = qL(i,IRHO) invRhoL = 1.0_RP / rhoL sqrtRhoL = sqrt(rhoL) - uL = ( qL(IRHOU) * n(IX) + qL(IRHOV) * n(IY) ) * invRhoL - vL = ( -qL(IRHOU) * n(IY) + qL(IRHOV) * n(IX) ) * invRhoL - pL = thermodynamics % gm1 * (qL(IRHOE) - 0.5_RP * rhoL * (uL * uL + vL * vL) ) - TL = dimensionless % gammaMach2 * pL * invRhoL + uL = ( qL(i,IRHOU) * normal(IX,i) + qL(i,IRHOV) * normal(IY,i) ) * invRhoL + vL = ( -qL(i,IRHOU) * normal(IY,i) + qL(i,IRHOV) * normal(IX,i) ) * invRhoL + pL = thermodynamics % gm1 * (qL(i,IRHOE) - 0.5_RP * rhoL * (uL * uL + vL * vL) ) HL = dimensionless % cp * pL * invRhoL + 0.5_RP * ( uL*uL + vL*vL ) - rhoR = qR(IRHO) + rhoR = qR(i,IRHO) invRhoR = 1.0_RP / rhoR sqrtRhoR = sqrt(rhoR) - uR = ( qR(IRHOU) * n(IX) + qR(IRHOV) * n(IY) ) * invRhoR - vR = ( -qR(IRHOU) * n(IY) + qR(IRHOV) * n(IX) ) * invRhoR - pR = thermodynamics % gm1 * (qR(IRHOE) - 0.5_RP * rhoR * (uR * uR + vR * vR) ) - TR = dimensionless % gammaMach2 * pR * invRhoR + uR = ( qR(i,IRHOU) * normal(IX,i) + qR(i,IRHOV) * normal(IY,i) ) * invRhoR + vR = ( -qR(i,IRHOU) * normal(IY,i) + qR(i,IRHOV) * normal(IX,i) ) * invRhoR + pR = thermodynamics % gm1 * (qR(i,IRHOE) - 0.5_RP * rhoR * (uR * uR + vR * vR) ) HR = dimensionless % cp * pR * invRhoR + 0.5_RP * ( uR*uR + vR*vR ) ! ! 1/ Compute Roe averages @@ -236,7 +253,7 @@ module pure function RoeFlux(qL, qR , n) result(Fstar) dq(IRHO) = rhoR - rhoL dq(IRHOU) = rhoR*uR - rhoL*uL dq(IRHOV) = rhoR*vR - rhoL*vL - dq(IRHOE) = qR(IRHOE) - qL(IRHOE) + dq(IRHOE) = qR(i,IRHOE) - qL(i,IRHOE) alpha(3) = dq(IRHOV) - v*dq(IRHO) alpha(2) = gm1 * (dq(IRHO) * (H - u*u) + u*dq(IRHOU) - dq(IRHOE) + (dq(IRHOV) - v*dq(IRHO))*v) / ( a*a ) @@ -246,60 +263,69 @@ module pure function RoeFlux(qL, qR , n) result(Fstar) ! ! 5/ Compute the flux ! ---------------- - Fstar = F_inviscidFlux(rhoL,uL,vL,pL,HL) + Fstar(i,:) = F_inviscidFlux(rhoL,uL,vL,pL,HL) do wave = 1 , negativeWaves - Fstar = Fstar + alpha(wave) * lambda(wave) * K(1:NCONS , wave) + Fstar(i,:) = Fstar(i,:) + alpha(wave) * lambda(wave) * K(1:NCONS , wave) end do ! ! 6/ Return to the 3D space ! ---------------------- - Ft = Fstar(IRHOU) * n(IX) - Fstar(IRHOV) * n(IY) - Fstar(IRHOV) = Fstar(IRHOU) * n(IY) + Fstar(IRHOV) * n(IX) - Fstar(IRHOU) = Ft + Ft = Fstar(i,IRHOU) * normal(IX,i) - Fstar(i,IRHOV) * normal(IY,i) + Fstar(i,IRHOV) = Fstar(i,IRHOU) * normal(IY,i) + Fstar(i,IRHOV) * normal(IX,i) + Fstar(i,IRHOU) = Ft + + end do #ifdef _DIMENSIONLESS_TAU - Fstar = Fstar * dimensionless % invSqrtGammaMach + Fstar = Fstar * dimensionless % invSqrtGammaMach #endif end function RoeFlux - module pure function HLLFlux(qL , qR , n) result(Fstar) + module pure function HLLFlux(N , qL , qR , normal) result(Fstar) use MatrixOperations implicit none - real(kind=RP), dimension(NCONS), intent(in) :: qL - real(kind=RP), dimension(NCONS), intent(in) :: qR - real(kind=RP), dimension(NDIM) , intent(in) :: n - real(kind=RP), dimension(NCONS) :: Fstar -! --------------------------------------------------------------- + integer , intent(in) :: N + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qL + real(kind=RP), dimension(0:N,NCONS), intent(in) :: qR + real(kind=RP), dimension(NDIM,0:N) , intent(in) :: normal + real(kind=RP), dimension(0:N,NCONS) :: Fstar +! +! --------------- +! Local variables +! --------------- +! real(kind=RP) :: rhoL , invRhoL , uL , vL , HL , aL , pL , rhoeL , sqrtRhoL real(kind=RP) :: rhoR , invRhoR , uR , vR , HR , aR , pR , rhoeR , sqrtRhoR real(kind=RP) :: invrho , u , v , H , a real(kind=RP) :: SL , SR , invdS real(kind=RP) :: Ft + integer :: i - + do i = 0 , N +! ! 0/ Gather variables ! ---------------- - rhoL = qL(IRHO) + rhoL = qL(i,IRHO) invRhoL = 1.0_RP / rhoL sqrtRhoL = sqrt(rhoL) - uL = ( qL(IRHOU) * n(IX) + qL(IRHOV) * n(IY) ) * invRhoL - vL = ( -qL(IRHOU) * n(IY) + qL(IRHOV) * n(IX) ) * invRhoL - pL = thermodynamics % gm1 * (qL(IRHOE) - 0.5_RP * rhoL * (uL * uL + vL * vL) ) + uL = ( qL(i,IRHOU) * normal(IX,i) + qL(i,IRHOV) * normal(IY,i) ) * invRhoL + vL = ( -qL(i,IRHOU) * normal(IY,i) + qL(i,IRHOV) * normal(IX,i) ) * invRhoL + pL = thermodynamics % gm1 * (qL(i,IRHOE) - 0.5_RP * rhoL * (uL * uL + vL * vL) ) HL = dimensionless % cp * pL * invRhoL + 0.5_RP * ( uL*uL + vL*vL ) aL = sqrt( thermodynamics % gamma * pL * invRhoL ) - rhoeL = qL(IRHOE) + rhoeL = qL(i,IRHOE) - rhoR = qR(IRHO) + rhoR = qR(i,IRHO) invRhoR = 1.0_RP / rhoR sqrtRhoR = sqrt(rhoR) - uR = ( qR(IRHOU) * n(IX) + qR(IRHOV) * n(IY) ) * invRhoR - vR = ( -qR(IRHOU) * n(IY) + qR(IRHOV) * n(IX) ) * invRhoR - pR = thermodynamics % gm1 * (qR(IRHOE) - 0.5_RP * rhoR * (uR * uR + vR * vR) ) + uR = ( qR(i,IRHOU) * normal(IX,i) + qR(i,IRHOV) * normal(IY,i) ) * invRhoR + vR = ( -qR(i,IRHOU) * normal(IY,i) + qR(i,IRHOV) * normal(IX,i) ) * invRhoR + pR = thermodynamics % gm1 * (qR(i,IRHOE) - 0.5_RP * rhoR * (uR * uR + vR * vR) ) HR = dimensionless % cp * pR * invRhoR + 0.5_RP * ( uR*uR + vR*vR ) aR = sqrt( thermodynamics % gamma * pR * invRhoR ) - rhoeR = qR(IRHOE) + rhoeR = qR(i,IRHOE) ! 1/ Compute Roe averages ! -------------------- @@ -323,28 +349,30 @@ module pure function HLLFlux(qL , qR , n) result(Fstar) ! 3/ Compute the fluxes depending on the speeds ! ------------------------------------------ if ( SL .ge. 0.0_RP ) then - Fstar = F_inviscidFlux(rhoL,uL,vL,pL,HL) + Fstar(i,:) = F_inviscidFlux(rhoL,uL,vL,pL,HL) elseif ( SR .le. 0.0_RP ) then - Fstar = F_inviscidFlux(rhoR,uR,vR,pR,HR) + Fstar(i,:) = F_inviscidFlux(rhoR,uR,vR,pR,HR) elseif ( (SL .lt. 0.0_RP) .or. (SR .gt. 0.0_RP) ) then invdS = 1.0_RP / (SR - SL) - Fstar(IRHO) = (SR*rhoL*uL - SL*rhoR*uR + SL*SR*(rhoR-rhoL)) * invdS - Fstar(IRHOU) = (SR*(rhoL*uL*uL+pL) - SL*(rhoR*uR*uR+pR) + SL*SR*(rhoR*uR-rhoL*uL)) * invdS - Fstar(IRHOV) = (SR*rhoL*uL*vL - SL*rhoR*uR*vR + SR*SL*(rhoR*vR-rhoL*vL)) * invdS - Fstar(IRHOE) = (SR*uL*rhoL*HL - SL*uR*rhoR*HR + SR*SL*(rhoeR-rhoeL)) * invdS + Fstar(i,IRHO) = (SR*rhoL*uL - SL*rhoR*uR + SL*SR*(rhoR-rhoL)) * invdS + Fstar(i,IRHOU) = (SR*(rhoL*uL*uL+pL) - SL*(rhoR*uR*uR+pR) + SL*SR*(rhoR*uR-rhoL*uL)) * invdS + Fstar(i,IRHOV) = (SR*rhoL*uL*vL - SL*rhoR*uR*vR + SR*SL*(rhoR*vR-rhoL*vL)) * invdS + Fstar(i,IRHOE) = (SR*uL*rhoL*HL - SL*uR*rhoR*HR + SR*SL*(rhoeR-rhoeL)) * invdS end if ! ! 6/ Return to the 3D space ! ---------------------- - Ft = Fstar(IRHOU) * n(IX) - Fstar(IRHOV) * n(IY) - Fstar(IRHOV) = Fstar(IRHOU) * n(IY) + Fstar(IRHOV) * n(IX) - Fstar(IRHOU) = Ft + Ft = Fstar(i,IRHOU) * normal(IX,i) - Fstar(i,IRHOV) * normal(IY,i) + Fstar(i,IRHOV) = Fstar(i,IRHOU) * normal(IY,i) + Fstar(i,IRHOV) * normal(IX,i) + Fstar(i,IRHOU) = Ft + + end do #ifdef _DIMENSIONLESS_TAU - Fstar = Fstar * dimensionless % invSqrtGammaMach + Fstar = Fstar * dimensionless % invSqrtGammaMach #endif end function HLLFlux diff --git a/Solver/src/Physics/ViscousFluxes.f90 b/Solver/src/Physics/ViscousFluxes.f90 index 7c4e6bd..fbdb11d 100644 --- a/Solver/src/Physics/ViscousFluxes.f90 +++ b/Solver/src/Physics/ViscousFluxes.f90 @@ -33,7 +33,7 @@ module pure function viscousFlux0D( q , dq) result(F) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU) - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV) @@ -42,19 +42,19 @@ module pure function viscousFlux0D( q , dq) result(F) dyT = gm1 * invRho * ( -q(IRHOE) * invRho * dq(IY,IRHO) + dq(IY,IRHOE) + u * ( u * dq(IY,IRHO) - dq(IY,IRHOU)) + v * (v * dq(IY,IRHO) - dq(IY,IRHOV)) ) F(IRHO,:) = 0.0_RP - F(IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU)) + lambda * mu * divV - F(IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV)) + lambda * mu * divV + F(IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU)) + lambda * divV + F(IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV)) + lambda * divV - F(IRHOU,IY) = - mu * vDivRho * dq(IX,IRHO) & - - mu * uDivRho * dq(IY,IRHO) & - + mu * invRho * dq(IY,IRHOU) & - + mu * invRho * dq(IX,IRHOV) + F(IRHOU,IY) = - vDivRho * dq(IX,IRHO) & + - uDivRho * dq(IY,IRHO) & + + invRho * dq(IY,IRHOU) & + + invRho * dq(IX,IRHOV) F(IRHOV,IX) = F(IRHOU,IY) - F(IRHOE,IX) = u * F(IRHOU,IX) + v * F(IRHOU,IY) + kappa * dxT - F(IRHOE,IY) = u * F(IRHOV,IX) + v * F(IRHOV,IY) + kappa * dyT + F(IRHOE,IX) = u * F(IRHOU,IX) + v * F(IRHOU,IY) + ( dimensionless % cp / dimensionless % Pr ) * dxT + F(IRHOE,IY) = u * F(IRHOV,IX) + v * F(IRHOV,IY) + ( dimensionless % cp / dimensionless % Pr ) * dyT end associate @@ -79,7 +79,7 @@ module pure function viscousFlux1D( N , q , dq ) result ( F ) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU) - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV) @@ -89,19 +89,19 @@ module pure function viscousFlux1D( N , q , dq ) result ( F ) dyT = gm1 * invRho * ( -q(:,IRHOE) * invRho * dq(:,IY,IRHO) + dq(:,IY,IRHOE) + u * ( u * dq(:,IY,IRHO) - dq(:,IY,IRHOU)) + v * (v * dq(:,IY,IRHO) - dq(:,IY,IRHOV)) ) F(:,IRHO,:) = 0.0_RP - F(:,IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU)) + lambda * mu * divV - F(:,IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV)) + lambda * mu * divV + F(:,IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU)) + lambda * divV + F(:,IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV)) + lambda * divV - F(:,IRHOU,IY) = - mu * vDivRho * dq(:,IX,IRHO) & - - mu * uDivRho * dq(:,IY,IRHO) & - + mu * invRho * dq(:,IY,IRHOU) & - + mu * invRho * dq(:,IX,IRHOV) + F(:,IRHOU,IY) = - vDivRho * dq(:,IX,IRHO) & + - uDivRho * dq(:,IY,IRHO) & + + invRho * dq(:,IY,IRHOU) & + + invRho * dq(:,IX,IRHOV) F(:,IRHOV,IX) = F(:,IRHOU,IY) - F(:,IRHOE,IX) = u * F(:,IRHOU,IX) + v * F(:,IRHOU,IY) + kappa * dxT - F(:,IRHOE,IY) = u * F(:,IRHOV,IX) + v * F(:,IRHOV,IY) + kappa * dyT + F(:,IRHOE,IX) = u * F(:,IRHOU,IX) + v * F(:,IRHOU,IY) + ( dimensionless % cp / dimensionless % Pr ) * dxT + F(:,IRHOE,IY) = u * F(:,IRHOV,IX) + v * F(:,IRHOV,IY) + ( dimensionless % cp / dimensionless % Pr ) * dyT end associate @@ -126,7 +126,7 @@ module pure function viscousFlux2D( N , q , dq ) result ( F ) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(:,:,IX,IRHO) + invRho * dq(:,:,IX,IRHOU) - vDivRho * dq(:,:,IY,IRHO) + invRho * dq(:,:,IY,IRHOV) @@ -135,19 +135,19 @@ module pure function viscousFlux2D( N , q , dq ) result ( F ) dyT = gm1 * invRho * ( -q(:,:,IRHOE) * invRho * dq(:,:,IY,IRHO) + dq(:,:,IY,IRHOE) + u * ( u * dq(:,:,IY,IRHO) - dq(:,:,IY,IRHOU)) + v * (v * dq(:,:,IY,IRHO) - dq(:,:,IY,IRHOV)) ) F(:,:,IRHO,:) = 0.0_RP - F(:,:,IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(:,:,IX,IRHO) + invRho * dq(:,:,IX,IRHOU)) + lambda * mu * divV - F(:,:,IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(:,:,IY,IRHO) + invRho * dq(:,:,IY,IRHOV)) + lambda * mu * divV + F(:,:,IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(:,:,IX,IRHO) + invRho * dq(:,:,IX,IRHOU)) + lambda * divV + F(:,:,IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(:,:,IY,IRHO) + invRho * dq(:,:,IY,IRHOV)) + lambda * divV - F(:,:,IRHOU,IY) = - mu * vDivRho * dq(:,:,IX,IRHO) & - - mu * uDivRho * dq(:,:,IY,IRHO) & - + mu * invRho * dq(:,:,IY,IRHOU) & - + mu * invRho * dq(:,:,IX,IRHOV) + F(:,:,IRHOU,IY) = - vDivRho * dq(:,:,IX,IRHO) & + - uDivRho * dq(:,:,IY,IRHO) & + + invRho * dq(:,:,IY,IRHOU) & + + invRho * dq(:,:,IX,IRHOV) F(:,:,IRHOV,IX) = F(:,:,IRHOU,IY) - F(:,:,IRHOE,IX) = u * F(:,:,IRHOU,IX) + v * F(:,:,IRHOU,IY) + kappa * dxT - F(:,:,IRHOE,IY) = u * F(:,:,IRHOV,IX) + v * F(:,:,IRHOV,IY) + kappa * dyT + F(:,:,IRHOE,IX) = u * F(:,:,IRHOU,IX) + v * F(:,:,IRHOU,IY) + ( dimensionless % cp / dimensionless % Pr ) * dxT + F(:,:,IRHOE,IY) = u * F(:,:,IRHOV,IX) + v * F(:,:,IRHOV,IY) + ( dimensionless % cp / dimensionless % Pr ) * dyT end associate @@ -172,7 +172,7 @@ module pure function viscousFluxBC0D( q , qB , dq) result(F) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU) - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV) @@ -182,19 +182,19 @@ module pure function viscousFluxBC0D( q , qB , dq) result(F) F(IRHO,:) = 0.0_RP - F(IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU)) + lambda * mu * divV - F(IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV)) + lambda * mu * divV + F(IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU)) + lambda * divV + F(IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV)) + lambda * divV - F(IRHOU,IY) = - mu * vDivRho * dq(IX,IRHO) & - - mu * uDivRho * dq(IY,IRHO) & - + mu * invRho * dq(IY,IRHOU) & - + mu * invRho * dq(IX,IRHOV) + F(IRHOU,IY) = - vDivRho * dq(IX,IRHO) & + - uDivRho * dq(IY,IRHO) & + + invRho * dq(IY,IRHOU) & + + invRho * dq(IX,IRHOV) F(IRHOV,IX) = F(IRHOU,IY) - F(IRHOE,IX) = ( qB(IRHOU) * F(IRHOU,IX) + qB(IRHOV) * F(IRHOU,IY) ) / qB(IRHO) + kappa * dxT - F(IRHOE,IY) = ( qB(IRHOU) * F(IRHOV,IX) + qB(IRHOV) * F(IRHOV,IY) ) / qB(IRHO) + kappa * dyT + F(IRHOE,IX) = ( qB(IRHOU) * F(IRHOU,IX) + qB(IRHOV) * F(IRHOU,IY) ) / qB(IRHO) + ( dimensionless % cp / dimensionless % Pr ) * dxT + F(IRHOE,IY) = ( qB(IRHOU) * F(IRHOV,IX) + qB(IRHOV) * F(IRHOV,IY) ) / qB(IRHO) + ( dimensionless % cp / dimensionless % Pr ) * dyT end associate @@ -220,7 +220,7 @@ module pure function viscousFluxBC1D( N , q , qb , dq ) result ( F ) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU) - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV) @@ -230,19 +230,19 @@ module pure function viscousFluxBC1D( N , q , qb , dq ) result ( F ) dyT = gm1 * invRho * ( -q(:,IRHOE) * invRho * dq(:,IY,IRHO) + dq(:,IY,IRHOE) + u * ( u * dq(:,IY,IRHO) - dq(:,IY,IRHOU)) + v * (v * dq(:,IY,IRHO) - dq(:,IY,IRHOV)) ) F(:,IRHO,:) = 0.0_RP - F(:,IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU)) + lambda * mu * divV - F(:,IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV)) + lambda * mu * divV + F(:,IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU)) + lambda * divV + F(:,IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV)) + lambda * divV - F(:,IRHOU,IY) = - mu * vDivRho * dq(:,IX,IRHO) & - - mu * uDivRho * dq(:,IY,IRHO) & - + mu * invRho * dq(:,IY,IRHOU) & - + mu * invRho * dq(:,IX,IRHOV) + F(:,IRHOU,IY) = - vDivRho * dq(:,IX,IRHO) & + - uDivRho * dq(:,IY,IRHO) & + + invRho * dq(:,IY,IRHOU) & + + invRho * dq(:,IX,IRHOV) F(:,IRHOV,IX) = F(:,IRHOU,IY) - F(:,IRHOE,IX) = (qB(:,IRHOU) * F(:,IRHOU,IX) + qB(:,IRHOV) * F(:,IRHOU,IY)) / qB(:,IRHO) + kappa * dxT - F(:,IRHOE,IY) = (qB(:,IRHOU) * F(:,IRHOV,IX) + qB(:,IRHOV) * F(:,IRHOV,IY)) / qB(:,IRHO) + kappa * dyT + F(:,IRHOE,IX) = (qB(:,IRHOU) * F(:,IRHOU,IX) + qB(:,IRHOV) * F(:,IRHOU,IY)) / qB(:,IRHO) + ( dimensionless % cp / dimensionless % Pr ) * dxT + F(:,IRHOE,IY) = (qB(:,IRHOU) * F(:,IRHOV,IX) + qB(:,IRHOV) * F(:,IRHOV,IY)) / qB(:,IRHO) + ( dimensionless % cp / dimensionless % Pr ) * dyT end associate @@ -267,20 +267,20 @@ module pure function AdiabaticViscousFlux0D( q , qB , dq) result(F) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU) - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV) F(IRHO,:) = 0.0_RP - F(IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU)) + lambda * mu * divV - F(IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV)) + lambda * mu * divV + F(IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(IX,IRHO) + invRho * dq(IX,IRHOU)) + lambda * divV + F(IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(IY,IRHO) + invRho * dq(IY,IRHOV)) + lambda * divV ! - - - F(IRHOU,IY) = - mu * vDivRho * dq(IX,IRHO) & ! du dv | u drho 1 drhou v drho 1 drhov | - - mu * uDivRho * dq(IY,IRHO) & ! mu -- + -- = mu | - --- ---- + --- ----- - --- ---- + --- ----- | - + mu * invRho * dq(IY,IRHOU) & ! dy dx | rho dy rho dy rho dx rho dx | - + mu * invRho * dq(IX,IRHOV) ! - - + F(IRHOU,IY) = - vDivRho * dq(IX,IRHO) & ! du dv | u drho 1 drhou v drho 1 drhov | + - uDivRho * dq(IY,IRHO) & ! mu -- + -- = mu | - --- ---- + --- ----- - --- ---- + --- ----- | + + invRho * dq(IY,IRHOU) & ! dy dx | rho dy rho dy rho dx rho dx | + + invRho * dq(IX,IRHOV) ! - - F(IRHOV,IX) = F(IRHOU,IY) @@ -312,19 +312,19 @@ module pure function AdiabaticViscousFlux1D( N , q , qB , dq ) result ( F ) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU) - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV) F(:,IRHO ,: ) = 0.0_RP - F(:,IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU)) + lambda * mu * divV - F(:,IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV)) + lambda * mu * divV + F(:,IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU)) + lambda * divV + F(:,IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV)) + lambda * divV - F(:,IRHOU,IY) = - mu * vDivRho * dq(:,IX,IRHO) & - - mu * uDivRho * dq(:,IY,IRHO) & - + mu * invRho * dq(:,IY,IRHOU) & - + mu * invRho * dq(:,IX,IRHOV) + F(:,IRHOU,IY) = - vDivRho * dq(:,IX,IRHO) & + - uDivRho * dq(:,IY,IRHO) & + + invRho * dq(:,IY,IRHOU) & + + invRho * dq(:,IX,IRHOV) F(:,IRHOV,IX) = F(:,IRHOU,IY) @@ -356,19 +356,19 @@ module pure function AdiabaticViscousFlux2D( N , q , qB , dq ) result ( F ) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(:,:,IX,IRHO) + invRho * dq(:,:,IX,IRHOU) - vDivRho * dq(:,:,IY,IRHO) + invRho * dq(:,:,IY,IRHOV) F(:,:,IRHO , :) = 0.0_RP - F(:,:,IRHOU,IX) = 2.0_RP * mu * ( - uDivRho * dq(:,:,IX,IRHO) + invRho * dq(:,:,IX,IRHOU)) + lambda * mu * divV - F(:,:,IRHOV,IY) = 2.0_RP * mu * ( - vDivRho * dq(:,:,IY,IRHO) + invRho * dq(:,:,IY,IRHOV)) + lambda * mu * divV + F(:,:,IRHOU,IX) = 2.0_RP * ( - uDivRho * dq(:,:,IX,IRHO) + invRho * dq(:,:,IX,IRHOU)) + lambda * divV + F(:,:,IRHOV,IY) = 2.0_RP * ( - vDivRho * dq(:,:,IY,IRHO) + invRho * dq(:,:,IY,IRHOV)) + lambda * divV - F(:,:,IRHOU,IY) = - mu * vDivRho * dq(:,:,IX,IRHO) & - - mu * uDivRho * dq(:,:,IY,IRHO) & - + mu * invRho * dq(:,:,IY,IRHOU) & - + mu * invRho * dq(:,:,IX,IRHOV) + F(:,:,IRHOU,IY) = - vDivRho * dq(:,:,IX,IRHO) & + - uDivRho * dq(:,:,IY,IRHO) & + + invRho * dq(:,:,IY,IRHOU) & + + invRho * dq(:,:,IX,IRHOV) F(:,:,IRHOV,IX) = F(:,:,IRHOU,IY) @@ -399,7 +399,7 @@ module pure function ComputeViscousTensor ( N , Q , dQ ) result ( tau ) uDivRho = u * invRho vDivRho = v * invRho - associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , kappa => dimensionless % kappa , gm1 => thermodynamics % gm1 ) + associate ( mu => dimensionless % mu , lambda => thermodynamics % lambda , gm1 => thermodynamics % gm1 ) divV = - uDivRho * dq(:,IX,IRHO) + invRho * dq(:,IX,IRHOU) - vDivRho * dq(:,IY,IRHO) + invRho * dq(:,IY,IRHOV) diff --git a/Solver/src/Plotting/Paraview.f90 b/Solver/src/Plotting/Paraview.f90 index 79878bd..c603736 100644 --- a/Solver/src/Plotting/Paraview.f90 +++ b/Solver/src/Plotting/Paraview.f90 @@ -14,6 +14,8 @@ character(len=STR_LEN_PLOTTER) :: str class(Charlist), pointer :: next => NULL() end type Charlist + + integer :: point_position = 0 ! ! ! ======== @@ -42,11 +44,21 @@ module subroutine Paraview_ExportMesh( self , mesh , Name ) self % Name = trim(auxname) call Paraview_OpenFile ( self , isMesh = .true. ) - +! +! Compute number of cells and points +! ---------------------------------- + self % npoints = 0 + self % ncells = 0 do eID = 1 , mesh % no_of_elements - call Paraview_NewMeshZone( self , mesh , eID ) + self % npoints = self % npoints + (mesh % elements(eID) % spA % N + 1)**2 + self % ncells = self % ncells + (mesh % elements(eID) % spA % N)**2 end do - +! + + write ( self % fID , '(A)') "DATASET UNSTRUCTURED_GRID" + + call Paraview_WriteMesh( self , mesh ) + close ( self % fID ) end subroutine Paraview_ExportMesh @@ -59,15 +71,34 @@ module subroutine Paraview_Export( self , mesh , Name) character(len=*) :: Name integer :: eID -! self % Name = trim(Name) + self % Name = trim(Name) // ".vtk" + + call Paraview_OpenFile ( self , isMesh = .false. ) ! -! call Paraview_OpenFile ( self , isMesh = .false. ) +! Compute number of cells and points +! ---------------------------------- + self % npoints = 0 + self % ncells = 0 + do eID = 1 , mesh % no_of_elements + self % npoints = self % npoints + (mesh % elements(eID) % spA % N + 1)**2 + self % ncells = self % ncells + (mesh % elements(eID) % spA % N)**2 + end do +! +! Create the zone +! --------------- + write ( self % fID , '(A)') "DATASET UNSTRUCTURED_GRID" +! +! Write the mesh +! -------------- + call Paraview_WriteMesh( self , mesh ) ! -! do eID = 1 , mesh % no_of_elements -! call Paraview_NewZone( self , mesh , eID ) -! end do +! Write the variables +! ------------------- + call Paraview_WriteVariables(self , mesh) ! -! close( self % fID ) +! Close file +! ---------- + close ( self % fID ) ! end subroutine Paraview_Export @@ -160,415 +191,320 @@ subroutine Paraview_OpenFile( self , IsMesh ) end subroutine Paraview_OpenFile - subroutine Paraview_NewMeshZone(self , mesh , eID ) + subroutine Paraview_WriteMesh(self , mesh ) use QuadMeshClass use Physics implicit none class(Paraview_t) :: self class(QuadMesh_t) :: mesh - integer :: eID ! ! --------------- ! Local variables ! --------------- ! - integer :: iXi , iEta - integer :: var - real(kind=RP) :: Q(1:NCONS) + integer :: eID + integer :: iXi , iEta + integer :: var + real(kind=RP) :: Q(1:NCONS) - associate ( N => mesh % elements(eID) % spA % N ) -! -! New header -! ---------- - write( self % fID , '(A)' ) "DATASET UNSTRUCTURED_GRID" - write( self % fID , '(A,I0,A)') "POINTS ",(N+1)*(N+1)," float" - - do iEta = 0 , N - do iXi = 0 , N - write( self % fID , '(ES17.10,1X,ES17.10,1X,ES17.10)') mesh % elements(eID) % x(iXi,iEta,IX) * RefValues % L & - , mesh % elements(eID) % x(iXi,iEta,IY) * RefValues % L & - , 0.0_RP + write ( self % fID , *) + write ( self % fID , '(A,I0,A)') "POINTS " , self % npoints , " float" + + + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N + do iXi = 0 , N + write( self % fID , '(ES17.10,1X,ES17.10,1X,ES17.10)') mesh % elements(eID) % x(iXi,iEta,IX) * RefValues % L & + , mesh % elements(eID) % x(iXi,iEta,IY) * RefValues % L & + , 0.0_RP + end do end do + end associate end do write( self % fID , * ) ! One blank line - write( self % fID , '(A,I0,1X,I0)' ) "CELLS ", N*N , 5*N*N - do iEta = 1 , N - do iXi = 1 , N - write(self % fID , '(I0,1X,I0,1X,I0,1X,I0,1X,I0)') 4,pointPosition(iXi,iEta,N) + + write( self % fID , '(A,I0,1X,I0)' ) "CELLS ", self % ncells,5*self % ncells + + point_position = -1 + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 1 , N + do iXi = 1 , N + write(self % fID , '(I0,1X,I0,1X,I0,1X,I0,1X,I0)') 4,pointPosition(iXi,iEta,N) + point_position + end do end do + point_position = point_position + (N+1)*(N+1) + end associate end do write( self % fID , * ) ! One blank line - write( self % fID , '(A,I0)' ) "CELL_TYPES ", N*N - do iEta = 1 , N - do iXi = 1 , N - write(self % fID , '(I0)') 9 + write( self % fID , '(A,I0)' ) "CELL_TYPES ", self % ncells + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 1 , N + do iXi = 1 , N + write(self % fID , '(I0)') 9 + end do end do + end associate end do - end associate - end subroutine Paraview_NewMeshZone + end subroutine Paraview_WriteMesh - subroutine Paraview_NewZone( self , mesh , eID , zoneType) + subroutine Paraview_WriteVariables( self , mesh ) use QuadMeshClass use Physics - use Setup_Class implicit none - class(Paraview_t) :: self - class(QuadMesh_t) :: mesh - integer :: eID - character(len=*), optional :: zoneType - character(len=STR_LEN_PLOTTER) :: zType -! -! associate ( N => mesh % elements(eID) % spA % N ) + class(Paraview_t), intent(in) :: self + class(QuadMesh_t), intent(in) :: mesh ! -! if ( present ( zoneType ) ) then -! zType = zoneType -! -! else -! zType = trim( Setup % outputType ) -! -! end if -! -! if ( trim(zType) .eq. "DGSEM" ) then -! call Paraview_StandardZone( self , mesh , eID ) -! elseif ( trim(zType) .eq. "Interpolated") then -! call Paraview_InterpolatedZone( self , mesh , eID ) -! else -! print*, "Unknown output type " , trim(Setup % outputType) -! print*, "Options available are: " -! print*, " * DGSEM" -! print*, " * Interpolated" -! stop "Stopped." -! end if -! -! end associate +! --------------- +! Local variables +! --------------- ! - end subroutine Paraview_NewZone - - subroutine Paraview_StandardZone(self , mesh , eID ) - use QuadMeshClass - use Physics - implicit none - class(Paraview_t) :: self - class(QuadMesh_t) :: mesh - integer :: eID -! -------------------------------------------------------------------------- - real(kind=RP), pointer :: rho (:,:) , rhou (:,:) , rhov (:,:) , rhoe (:,:) - real(kind=RP), pointer :: rhot(:,:) , rhout(:,:) , rhovt(:,:) , rhoet(:,:) + integer :: iVar , eID , iXi , iEta + real(kind=RP) :: Q(NCONS) + real(kind=RP) :: dQ(NDIM,NCONS) + + write(self % fID,*) + write(self % fID , '(A,I0)') "POINT_DATA " , self % npoints + do iVar = 1 , self % no_of_variables + + select case ( trim( self % variables(iVar) ) ) + + case ("rho") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rho float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % Q(iXi,iEta,IRHO) * refValues % rho + end do ; end do + end associate + end do + + case ("rhou") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhou float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % Q(iXi,iEta,IRHOU) * refValues % rho * refValues % a + end do ; end do + end associate + end do + + case ("rhov") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhov float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % Q(iXi,iEta,IRHOV) * refValues % rho * refValues % a + end do ; end do + end associate + end do + + case ("rhoe") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhoe float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % Q(iXi,iEta,IRHOE) * refValues % p + end do ; end do + end associate + end do + + case ("rhot") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhot float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % QDot(iXi,iEta,IRHO) * refValues % rho / refValues % tc + end do ; end do + end associate + end do + + case ("rhout") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhout float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % QDot(iXi,iEta,IRHOU) * refValues % rho * refValues % a / refValues % tc + end do ; end do + end associate + end do + + case ("rhovt") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhovt float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % QDot(iXi,iEta,IRHOV) * refValues % rho * refValues % a / refValues % tc + end do ; end do + end associate + end do + + case ("rhoet") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS rhoet float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % QDot(iXi,iEta,IRHOE) * refValues % p / refValues % tc + end do ; end do + end associate + end do + + case ("u") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS u float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % Q(iXi,iEta,IRHOU) / mesh % elements(eID) % Q(iXi,iEta,IRHO) * refValues % a + end do ; end do + end associate + end do + + case ("v") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS v float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + write(self % fID , '(ES17.10)') mesh % elements(eID) % Q(iXi,iEta,IRHOV) / mesh % elements(eID) % Q(iXi,iEta,IRHO) * refValues % a + end do ; end do + end associate + end do + + case ("p") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS p float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + write(self % fID , '(ES17.10)') getPressure(Q) * refValues % p + end do ; end do + end associate + end do + + case ("Mach") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS Mach float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + write(self % fID , '(ES17.10)') ( sqrt(Q(IRHOU)*Q(IRHOU)+Q(IRHOV)*Q(IRHOV))/ Q(IRHO)) / getSoundSpeed(Q) + end do ; end do + end associate + end do + + case ("s") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS s float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + write(self % fID , '(ES17.10)') getPressure(Q) * refValues % p / (Q(IRHO) * refValues % rho) ** ( Thermodynamics % gamma) + end do ; end do + end associate + end do + #ifdef NAVIER_STOKES - real(kind=RP), target :: du(0:mesh % elements(eID) % spA % N , 0:mesh % elements(eID) % spA % N , 1:NDIM , 1:NDIM) - real(kind=RP), pointer :: ux(:,:) , uy(:,:) , vx(:,:) , vy(:,:) -#endif - integer :: iXi , iEta - integer :: var - real(kind=RP) :: Q(1:NCONS) - -! associate ( N => mesh % elements(eID) % spA % N ) -!! -!! New header -!! ---------- -! write( self % fID , '(A,I0,A)' , advance="no" ) "ZONE N=",(N+1)*(N+1),", " -! write( self % fID , '(A,I0,A)' , advance="no" ) "E=",(N)*(N),", " -! write( self % fID , '(A)' ) "DATAPACKING=POINT, ZONETYPE=FEQUADRILATERAL" -!! -!! Point to the quantities -!! ----------------------- -! rho(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHO) -! rhou(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHOU) -! rhov(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHOV) -! rhoe(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHOE) -! rhot(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHO) -! rhout(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHOU) -! rhovt(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHOV) -! rhoet(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHOE) -!#ifdef NAVIER_STOKES -! if ( self % no_of_variables .ne. 0 ) then -! du = getStrainTensor( N , mesh % elements(eID) % Q , mesh % elements(eID) % dQ ) -! ux(0:,0:) => du(0:,0:,IX,IX) -! uy(0:,0:) => du(0:,0:,IY,IX) -! vx(0:,0:) => du(0:,0:,IX,IY) -! vy(0:,0:) => du(0:,0:,IY,IY) -! end if -!#endif -! -! -! do iEta = 0 , N -! do iXi = 0 , N -! write( self % fID , '(ES17.10,1X,ES17.10,1X,ES17.10)',advance="no") mesh % elements(eID) % x(iXi,iEta,IX) * RefValues % L & -! , mesh % elements(eID) % x(iXi,iEta,IY) * RefValues % L & -! , 0.0_RP -!! -!! Save quantities -!! --------------- -! do var = 1 , self % no_of_variables -! -! select case ( trim( self % variables(var) ) ) -! case ("rho") -! write(self % fID,'(1X,ES17.10)',advance="no") rho(iXi,iEta) * refValues % rho -! -! case ("rhou") -! write(self % fID,'(1X,ES17.10)',advance="no") rhou(iXi,iEta) * refValues % rho * refValues % a -! -! case ("rhov") -! write(self % fID,'(1X,ES17.10)',advance="no") rhov(iXi,iEta) * refValues % rho * refValues % a -! -! case ("rhoe") -! write(self % fID,'(1X,ES17.10)',advance="no") rhoe(iXi,iEta) * refValues % rho * refValues % p -! -! case ("rhot") -! write(self % fID,'(1X,ES17.10)',advance="no") rhot(iXi,iEta) * refValues % rho / refValues % tc -! -! case ("rhout") -! write(self % fID,'(1X,ES17.10)',advance="no") rhout(iXi,iEta) * refValues % rho * refValues % a / refValues % tc -! -! case ("rhovt") -! write(self % fID,'(1X,ES17.10)',advance="no") rhovt(iXi,iEta) * refValues % rho * refValues % a / refValues % tc -! -! case ("rhoet") -! write(self % fID,'(1X,ES17.10)',advance="no") rhoet(iXi,iEta) * refValues % rho * refValues % p / refValues % tc -! -! case ("u") -! write(self % fID,'(1X,ES17.10)',advance="no") rhou(iXi,iEta)/rho(iXi,iEta) * refValues % a -! -! case ("v") -! write(self % fID,'(1X,ES17.10)',advance="no") rhov(iXi,iEta)/rho(iXi,iEta) * refValues % a -! -! case ("p") -! Q = [rho(iXi,iEta) , rhou(iXi,iEta) , rhov(iXi,iEta) , rhoe(iXi,iEta) ] -! write(self % fID,'(1X,ES17.10)',advance="no") getPressure( Q ) * refValues % p -! -! case ("Mach") -! Q = [rho(iXi,iEta) , rhou(iXi,iEta) , rhov(iXi,iEta) , rhoe(iXi,iEta) ] -! write(self % fID,'(1X,ES17.10)',advance="no") sqrt(rhou(iXi,iEta)*rhou(iXi,iEta)+rhov(iXi,iEta)*rhov(iXi,iEta))/rho(iXi,iEta)/ getSoundSpeed( Q ) -! -! case ("s") -! write(self % fID,'(1X,ES17.10)',advance="no") Thermodynamics % gm1 * ( rhoe(iXi,iEta) - & -! 0.5*rhou(iXi,iEta)*rhou(iXi,iEta)/rho(iXi,iEta) - 0.5*rhov(iXi,iEta)*rhov(iXi,iEta)/rho(iXi,iEta) ) & -! / (rho(iXi,iEta) * refValues % rho)**(Thermodynamics % gamma) * refValues % p -!#ifdef NAVIER_STOKES -! case ( "ux" ) -! write(self % fID,'(1X,ES17.10)',advance="no") ux(iXi,iEta) * refValues % a -! -! case ( "uy" ) -! write(self % fID,'(1X,ES17.10)',advance="no") uy(iXi,iEta) * refValues % a -! -! case ( "vx" ) -! write(self % fID,'(1X,ES17.10)',advance="no") ux(iXi,iEta) * refValues % a -! -! case ( "vy" ) -! write(self % fID,'(1X,ES17.10)',advance="no") vy(iXi,iEta) * refValues % a -! -! case ("vort") -! write(self % fID,'(1X,ES17.10)',advance="no") ( vx(iXi,iEta) - uy(iXi,iEta) ) * refValues % a / refValues % L -!#endif -! -! end select -! -! end do -! -!! Jump to next line -!! ----------------- -! write( self % fID , *) -! -! end do -! end do -! -! write( self % fID , * ) ! One blank line -! -! do iEta = 1 , N -! do iXi = 1 , N -! write(self % fID , '(I0,1X,I0,1X,I0,1X,I0)') pointPosition(iXi,iEta,N) -! end do -! end do -! -! end associate -! - end subroutine Paraview_StandardZone - subroutine Paraview_InterpolatedZone(self , mesh , eID ) - use QuadMeshClass - use Physics - use Setup_Class - use MatrixOperations - use InterpolationAndDerivatives - implicit none - class(Paraview_t) :: self - class(QuadMesh_t) :: mesh - integer :: eID -! -------------------------------------------------------------------------- - real(kind=RP), pointer :: rhoDG (:,:) , rhouDG (:,:) , rhovDG (:,:) , rhoeDG (:,:) - real(kind=RP), pointer :: rhotDG (:,:) , rhoutDG (:,:) , rhovtDG (:,:) , rhoetDG (:,:) - real(kind=RP), pointer :: rho (:,:) , rhou (:,:) , rhov (:,:) , rhoe (:,:) - real(kind=RP), pointer :: rhot (:,:) , rhout (:,:) , rhovt (:,:) , rhoet (:,:) -#ifdef NAVIER_STOKES - real(kind=RP), target :: du(0:mesh % elements(eID) % spA % N , 0:mesh % elements(eID) % spA % N , 1:NDIM , 1:NDIM) - real(kind=RP), pointer :: uxDG (:,:) , uyDG (:,:) , vxDG (:,:) , vyDG (:,:) - real(kind=RP), pointer :: ux (:,:) , uy (:,:) , vx (:,:) , vy (:,:) + case ("ux") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS ux float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + dQ = mesh % elements(eID) % dQ(iXi,iEta,:,:) + write(self % fID , '(ES17.10)') (- Q(IRHOU) / Q(IRHO) * dQ(IX,IRHO) + dQ(IX,IRHOU))/Q(IRHO) + end do ; end do + end associate + end do + + case ("uy") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS uy float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + dQ = mesh % elements(eID) % dQ(iXi,iEta,:,:) + write(self % fID , '(ES17.10)') (- Q(IRHOU) / Q(IRHO) * dQ(IY,IRHO) + dQ(IY,IRHOU))/Q(IRHO) + end do ; end do + end associate + end do + + case ("vx") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS vx float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + dQ = mesh % elements(eID) % dQ(iXi,iEta,:,:) + write(self % fID , '(ES17.10)') (- Q(IRHOV) / Q(IRHO) * dQ(IX,IRHO) + dQ(IX,IRHOV))/Q(IRHO) + end do ; end do + end associate + end do + + case ("vy") + write(self % fID,*) + write(self % fID , '(A)') "SCALARS vy float" + write(self % fID , '(A)') "LOOKUP_TABLE default" + do eID = 1 , mesh % no_of_elements + associate ( N => mesh % elements(eID) % spA % N ) + do iEta = 0 , N ; do iXi = 0 , N + Q = mesh % elements(eID) % Q(iXi,iEta,:) + dQ = mesh % elements(eID) % dQ(iXi,iEta,:,:) + write(self % fID , '(ES17.10)') (- Q(IRHOV) / Q(IRHO) * dQ(IY,IRHO) + dQ(IY,IRHOV))/Q(IRHO) + end do ; end do + end associate + end do + #endif - real(kind=RP), allocatable :: xi(:) , T(:,:) , x(:) - integer :: iXi , iEta - integer :: Nout - integer :: var - real(kind=RP) :: Q(1:NCONS) - -! associate ( N => mesh % elements(eID) % spA % N , spA => mesh % elements(eID) % spA) -!! -!! Construct the interpolation framework -!! ------------------------------------- -! Nout = Setup % no_of_plotPoints - 1 -! -!! -!! New header -!! ---------- -! write( self % fID , '(A,I0,A)' , advance="no" ) "ZONE N=",(Nout+1)*(Nout+1),", " -! write( self % fID , '(A,I0,A)' , advance="no" ) "E=",(Nout)*(Nout),", " -! write( self % fID , '(A)' ) "DATAPACKING=POINT, ZONETYPE=FEQUADRILATERAL" -! -! allocate( xi ( 0 : Nout ) , T (0 : Nout , 0 : N ) , x(NDIM)) -! -! xi = reshape( (/( (1.0_RP * iXi) / Nout , iXi = 0 , Nout )/) , (/Nout + 1/) ) -! call PolynomialInterpolationMatrix( N , Nout , spA % xi , spA % wb , xi , T ) -!! -!! Point to the quantities -!! ----------------------- -! rhoDG(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHO) -! rhouDG(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHOU) -! rhovDG(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHOV) -! rhoeDG(0:,0:) => mesh % elements(eID) % Q(0:,0:,IRHOE) -! rhotDG(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHO) -! rhoutDG(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHOU) -! rhovtDG(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHOV) -! rhoetDG(0:,0:) => mesh % elements(eID) % QDot(0:,0:,IRHOE) -! -!#ifdef NAVIER_STOKES -! if ( self % no_of_variables .ne. 0 ) then -! du = getStrainTensor( N , mesh % elements(eID) % Q , mesh % elements(eID) % dQ ) -! uxDG(0:,0:) => du(0:,0:,IX,IX) -! uyDG(0:,0:) => du(0:,0:,IY,IX) -! vxDG(0:,0:) => du(0:,0:,IX,IY) -! vyDG(0:,0:) => du(0:,0:,IY,IY) -! end if -!#endif -!! -!! Obtain the interpolation to a new set of equispaced points -!! ---------------------------------------------------------- -! allocate ( rho(0:Nout , 0:Nout) , rhou(0:Nout , 0:Nout) , rhov(0:Nout , 0:Nout) , rhoe(0:Nout , 0:Nout) ) -! allocate ( rhot(0:Nout , 0:Nout) , rhout(0:Nout , 0:Nout) , rhovt(0:Nout , 0:Nout) , rhoet(0:Nout , 0:Nout) ) -!#ifdef NAVIER_STOKES -! allocate ( ux(0:Nout , 0:Nout) , uy(0:Nout , 0:Nout) , vx(0:Nout , 0:Nout) , vy(0:Nout,0:Nout) ) -!#endif -! call TripleMatrixProduct ( T , rhoDG , T , rho , trC = .true. ) -! call TripleMatrixProduct ( T , rhouDG , T , rhou , trC = .true. ) -! call TripleMatrixProduct ( T , rhovDG , T , rhov , trC = .true. ) -! call TripleMatrixProduct ( T , rhoeDG , T , rhoe , trC = .true. ) -! call TripleMatrixProduct ( T , rhotDG , T , rhot , trC = .true. ) -! call TripleMatrixProduct ( T , rhoutDG , T , rhout , trC = .true. ) -! call TripleMatrixProduct ( T , rhovtDG , T , rhovt , trC = .true. ) -! call TripleMatrixProduct ( T , rhoetDG , T , rhoet , trC = .true. ) -!#ifdef NAVIER_STOKES -! call TripleMatrixProduct ( T , uxDG , T , ux , trC = .true. ) -! call TripleMatrixProduct ( T , uyDG , T , uy , trC = .true. ) -! call TripleMatrixProduct ( T , vxDG , T , vx , trC = .true. ) -! call TripleMatrixProduct ( T , vyDG , T , vy , trC = .true. ) -!#endif -! -! do iEta = 0 , Nout -! do iXi = 0 , Nout -! -! x = mesh % elements(eID) % compute_X ( xi(iXi) , xi(iEta) ) -! write( self % fID , '(ES17.10,1X,ES17.10,1X,ES17.10)',advance="no") x(IX) * RefValues % L & -! , x(IY) * RefValues % L & -! , 0.0_RP -!! -!! Save quantities -!! --------------- -! do var = 1 , self % no_of_variables -! -! select case ( trim( self % variables(var) ) ) -! case ("rho") -! write(self % fID,'(1X,ES17.10)',advance="no") rho(iXi,iEta) * refValues % rho -! -! case ("rhou") -! write(self % fID,'(1X,ES17.10)',advance="no") rhou(iXi,iEta) * refValues % rho * refValues % a -! -! case ("rhov") -! write(self % fID,'(1X,ES17.10)',advance="no") rhov(iXi,iEta) * refValues % rho * refValues % a -! -! case ("rhoe") -! write(self % fID,'(1X,ES17.10)',advance="no") rhoe(iXi,iEta) * refValues % rho * refValues % p -! -! case ("rhot") -! write(self % fID,'(1X,ES17.10)',advance="no") rhot(iXi,iEta) * refValues % rho / refValues % tc -! -! case ("rhout") -! write(self % fID,'(1X,ES17.10)',advance="no") rhout(iXi,iEta) * refValues % rho * refValues % a / refValues % tc -! -! case ("rhovt") -! write(self % fID,'(1X,ES17.10)',advance="no") rhovt(iXi,iEta) * refValues % rho * refValues % a / refValues % tc -! -! case ("rhoet") -! write(self % fID,'(1X,ES17.10)',advance="no") rhoet(iXi,iEta) * refValues % rho * refValues % p / refValues % tc -! -! case ("u") -! write(self % fID,'(1X,ES17.10)',advance="no") rhou(iXi,iEta)/rho(iXi,iEta) * refValues % a -! -! case ("v") -! write(self % fID,'(1X,ES17.10)',advance="no") rhov(iXi,iEta)/rho(iXi,iEta) * refValues % a -! -! case ("p") -! Q = [rho(iXi,iEta) , rhou(iXi,iEta) , rhov(iXi,iEta) , rhoe(iXi,iEta) ] -! write(self % fID,'(1X,ES17.10)',advance="no") getPressure( Q ) * refValues % p -! -! case ("Mach") -! Q = [rho(iXi,iEta) , rhou(iXi,iEta) , rhov(iXi,iEta) , rhoe(iXi,iEta) ] -! write(self % fID,'(1X,ES17.10)',advance="no") sqrt(rhou(iXi,iEta)*rhou(iXi,iEta)+rhov(iXi,iEta)*rhov(iXi,iEta))/rho(iXi,iEta)/ getSoundSpeed( Q ) -! -! case ("s") -! write(self % fID,'(1X,ES17.10)',advance="no") Thermodynamics % gm1 * ( rhoe(iXi,iEta) - & -! 0.5*rhou(iXi,iEta)*rhou(iXi,iEta)/rho(iXi,iEta) - 0.5*rhov(iXi,iEta)*rhov(iXi,iEta)/rho(iXi,iEta) ) & -! / (rho(iXi,iEta) * refValues % rho)**(Thermodynamics % gamma) * refValues % p -!#ifdef NAVIER_STOKES -! case ( "ux" ) -! write(self % fID,'(1X,ES17.10)',advance="no") ux(iXi,iEta) * refValues % a -! -! case ( "uy" ) -! write(self % fID,'(1X,ES17.10)',advance="no") uy(iXi,iEta) * refValues % a -! -! case ( "vx" ) -! write(self % fID,'(1X,ES17.10)',advance="no") ux(iXi,iEta) * refValues % a -! -! case ( "vy" ) -! write(self % fID,'(1X,ES17.10)',advance="no") vy(iXi,iEta) * refValues % a -! -! case ("vort") -! write(self % fID,'(1X,ES17.10)',advance="no") ( vx(iXi,iEta) - uy(iXi,iEta) ) * refValues % a / refValues % L -!#endif -! -! end select -! -! end do -! -!! Jump to next line -!! ----------------- -! write( self % fID , *) -! -! end do -! end do -! -! write( self % fID , * ) ! One blank line + + + end select + end do + + end subroutine Paraview_WriteVariables ! -! do iEta = 1 , Nout -! do iXi = 1 , Nout -! write(self % fID , '(I0,1X,I0,1X,I0,1X,I0)') pointPosition(iXi,iEta,Nout) -! end do -! end do +!/////////////////////////////////////////////////////////////////////////////////////////// ! -! end associate +! AUXILIAR SUBROUTINES +! -------------------- +!/////////////////////////////////////////////////////////////////////////////////////////// ! - end subroutine Paraview_InterpolatedZone - function pointPosition(iXi , iEta , N) result( val ) implicit none integer :: iXi @@ -576,10 +512,10 @@ function pointPosition(iXi , iEta , N) result( val ) integer :: N integer :: val(POINTS_PER_QUAD) - val(1) = (N+1)*(iEta-1) + iXi - val(2) = (N+1)*(iEta-1) + iXi -1 - val(3) = (N+1)*iEta + iXi -1 - val(4) = (N+1)*iEta + iXi + val(1) = (N+1)*(iEta-1) + iXi + 1 + val(2) = (N+1)*(iEta-1) + iXi + val(3) = (N+1)*iEta + iXi + val(4) = (N+1)*iEta + iXi + 1 end function pointPosition subroutine LinkedList_Destruct( self ) diff --git a/Solver/src/Plotting/Plotter.f90 b/Solver/src/Plotting/Plotter.f90 index 295b500..081e984 100644 --- a/Solver/src/Plotting/Plotter.f90 +++ b/Solver/src/Plotting/Plotter.f90 @@ -29,6 +29,8 @@ module Plotter type, extends(Plotter_t) :: Paraview_t integer :: no_of_variables = 0 + integer :: npoints + integer :: ncells integer :: fID character(len=STR_LEN_PLOTTER), allocatable :: variables(:) character(len=STR_LEN_PLOTTER) :: Name @@ -94,14 +96,25 @@ end subroutine Paraview_Export ! ======== ! subroutine ConstructPlotter( self ) + use Setup_class implicit none class(Plotter_t), allocatable :: self ! ! TODO: Select format in case file ! -------------------------------- - allocate ( Tecplot_t :: self ) - !allocate ( Paraview_t :: self ) + if ( (trim(Setup % exportFormat) .eq. "Tecplot") .or. (trim(Setup % exportFormat) .eq. "plt") ) then + allocate ( Tecplot_t :: self ) + + elseif (( trim(Setup % exportFormat) .eq. "Paraview") .or. (trim(Setup % exportFormat) .eq. "vtk" ) ) then + allocate ( Paraview_t :: self ) + + else + print*, "Plotter was not loaded" + errorMessage(STD_OUT) + stop "Stopped" + + end if call self % Construct diff --git a/Solver/src/Plotting/Tecplot.f90 b/Solver/src/Plotting/Tecplot.f90 index 3f446f6..2b0ef08 100644 --- a/Solver/src/Plotting/Tecplot.f90 +++ b/Solver/src/Plotting/Tecplot.f90 @@ -319,7 +319,7 @@ subroutine Tecplot_StandardZone(self , mesh , eID ) write(self % fID,'(1X,ES17.10)',advance="no") rhov(iXi,iEta) * refValues % rho * refValues % a case ("rhoe") - write(self % fID,'(1X,ES17.10)',advance="no") rhoe(iXi,iEta) * refValues % rho * refValues % p + write(self % fID,'(1X,ES17.10)',advance="no") rhoe(iXi,iEta) * refValues % p case ("rhot") write(self % fID,'(1X,ES17.10)',advance="no") rhot(iXi,iEta) * refValues % rho / refValues % tc @@ -331,7 +331,7 @@ subroutine Tecplot_StandardZone(self , mesh , eID ) write(self % fID,'(1X,ES17.10)',advance="no") rhovt(iXi,iEta) * refValues % rho * refValues % a / refValues % tc case ("rhoet") - write(self % fID,'(1X,ES17.10)',advance="no") rhoet(iXi,iEta) * refValues % rho * refValues % p / refValues % tc + write(self % fID,'(1X,ES17.10)',advance="no") rhoet(iXi,iEta) * refValues % p / refValues % tc case ("u") write(self % fID,'(1X,ES17.10)',advance="no") rhou(iXi,iEta)/rho(iXi,iEta) * refValues % a @@ -366,6 +366,9 @@ subroutine Tecplot_StandardZone(self , mesh , eID ) case ("vort") write(self % fID,'(1X,ES17.10)',advance="no") ( vx(iXi,iEta) - uy(iXi,iEta) ) * refValues % a / refValues % L + + case ("muart") + write(self % fID,'(1X,ES17.10)',advance="no") mesh % elements(eID) % mu_a #endif end select @@ -501,7 +504,7 @@ subroutine Tecplot_InterpolatedZone(self , mesh , eID ) write(self % fID,'(1X,ES17.10)',advance="no") rhov(iXi,iEta) * refValues % rho * refValues % a case ("rhoe") - write(self % fID,'(1X,ES17.10)',advance="no") rhoe(iXi,iEta) * refValues % rho * refValues % p + write(self % fID,'(1X,ES17.10)',advance="no") rhoe(iXi,iEta) * refValues % p case ("rhot") write(self % fID,'(1X,ES17.10)',advance="no") rhot(iXi,iEta) * refValues % rho / refValues % tc @@ -513,7 +516,7 @@ subroutine Tecplot_InterpolatedZone(self , mesh , eID ) write(self % fID,'(1X,ES17.10)',advance="no") rhovt(iXi,iEta) * refValues % rho * refValues % a / refValues % tc case ("rhoet") - write(self % fID,'(1X,ES17.10)',advance="no") rhoet(iXi,iEta) * refValues % rho * refValues % p / refValues % tc + write(self % fID,'(1X,ES17.10)',advance="no") rhoet(iXi,iEta) * refValues % p / refValues % tc case ("u") write(self % fID,'(1X,ES17.10)',advance="no") rhou(iXi,iEta)/rho(iXi,iEta) * refValues % a @@ -548,6 +551,9 @@ subroutine Tecplot_InterpolatedZone(self , mesh , eID ) case ("vort") write(self % fID,'(1X,ES17.10)',advance="no") ( vx(iXi,iEta) - uy(iXi,iEta) ) * refValues % a / refValues % L + + case ("muart") + write(self % fID,'(1X,ES17.10)',advance="no") mesh % elements(eID) % mu_a #endif end select diff --git a/TestCases/Benchmarktests/FlowOverCircle/SETUP/CylinderRe1000.HiOCase b/TestCases/Benchmarktests/FlowOverCircle/SETUP/CylinderRe1000.HiOCase index ae9b315..9be03b0 100644 --- a/TestCases/Benchmarktests/FlowOverCircle/SETUP/CylinderRe1000.HiOCase +++ b/TestCases/Benchmarktests/FlowOverCircle/SETUP/CylinderRe1000.HiOCase @@ -34,7 +34,7 @@ !------------------------------------: Viscous discretization Viscous discretization: IP Interior penalty method: SIPG - Jumps penalty parameter: 5.0 + Jumps penalty parameter: 10.0 Gradient jumps penalty parameter: 0.0 !-----------------------------------: Time integration @@ -61,10 +61,10 @@ ! # define p-Refinement - 1: 2 - 2: 8 - 3: 6 - 4: 4 + 1: 6 + 2: 10 + 3: 8 + 4: 8 # end diff --git a/TestCases/Benchmarktests/Nozzle/SETUP/Nozzle.HiOCase b/TestCases/Benchmarktests/Nozzle/SETUP/Nozzle.HiOCase index c2cd2c7..69166b0 100644 --- a/TestCases/Benchmarktests/Nozzle/SETUP/Nozzle.HiOCase +++ b/TestCases/Benchmarktests/Nozzle/SETUP/Nozzle.HiOCase @@ -1,36 +1,36 @@ ! ! *************************** -!  * Cylinder parameter file * -! *************************** +!  * Cylinder parameter file * +! *************************** ! !-------------------------------------: Geometry - Mesh file : ./MESH/Nozzle.HiOMesh - Boundary file : _this_ + Mesh file : ./MESH/Nozzle.HiOMesh + Boundary file : _this_ !-----------------------------------------: Physics - Gas : Air - Reference pressure : 101325.0 - Reference Temperature : 273.15 - Reynolds length : 1.0 - Reynolds number : 1600.0 - Prandtl number : 0.72 - Mach number : 0.01 - + Gas : Air + Reference pressure : 101325.0 + Reference Temperature : 273.15 + Reynolds length : 1.0 + Reynolds number : 0.0 + Prandtl number : 0.72 + Mach number : 0.01 + !------------------------------------: DG Discretization - Interpolation nodes : Legendre-Gauss - Default polynomial order : 4 + Interpolation nodes : Legendre-Gauss + Default polynomial order : 4 !------------------------------------: Initialization - Initial condition : Steady - Restart file : ./RESULTS/NozzleSteady_46000000.HiORst + Initial condition : Steady + Restart file : ./RESULTS/NozzleSteady_46000000.HiORst !---------------------------------------------: Inviscid discretization - Inviscid strategy : Standard - Inviscid formulation : Form I - Inviscid Riemann Flux : Roe - ! Number of integration points : 6 + Inviscid discretization : Standard + Inviscid formulation : Green form + Inviscid Riemann Flux : Roe + ! Number of integration points : 6 !---------------------------------------------: Viscous discretization Viscous strategy : Interior-penalty @@ -38,6 +38,11 @@ Jumps penalty parameter : 10.0 Gradient jumps penalty parameter : 0.0 +!---------------------------------------------: Artificial dissipation + Artificial dissipation (0/1) : 1 + Artificial dissipation type : Green volume + Artificial dissipation matrix : Residuals-based + !----------------------------------------: Time integration Integration mode : Steady Integration scheme : Williamson RK5 @@ -48,40 +53,40 @@ Initial time : 0.0 !-----------------------------------: Output variables - Solution file : ./RESULTS/NozzleSteady.HiORst - Autosave interval : 10000 - Output interval : 10 + Solution file : ./RESULTS/NozzleSteady.HiORst + Autosave interval : 10000 + Output interval : 10 Output file type : Interpolated Number of representation points : 10 - Save variables : rho_rhou_rhov_rhoe_rhot_rhout_rhovt_rhoet_u_v_p_Mach + Save variables : rho_rhou_rhov_rhoe_rhot_rhout_rhovt_rhoet_u_v_p_Mach ! ! ********************************** -!  * Cylinder boundary mapping file * -! ********************************** +!  * Cylinder boundary mapping file * +! ********************************** ! # define zone 1 - Name = Nozzle-bottom - Type = Euler wall + Name = Nozzle-bottom + Type = Euler wall Riemann solver = Exact # end # define zone 2 - Name = Nozzle-top - Type = Euler wall + Name = Nozzle-top + Type = Euler wall Riemann solver = Exact # end # define zone 3 - Name = Outflow - Type = Pressure outlet + Name = Outflow + Type = Pressure outlet Riemann solver = Exact # end # define zone 4 - Name = Inflow - Type = Pressure inlet + Name = Inflow + Type = Pressure inlet # end diff --git a/TestCases/SupersonicTests/Cylinder/MESH/Cylinder.HiOMesh b/TestCases/SupersonicTests/Cylinder/MESH/Cylinder.HiOMesh new file mode 100644 index 0000000..0d9bc9f Binary files /dev/null and b/TestCases/SupersonicTests/Cylinder/MESH/Cylinder.HiOMesh differ diff --git a/TestCases/SupersonicTests/Cylinder/SETUP/CylinderSupersonic.HiOCase b/TestCases/SupersonicTests/Cylinder/SETUP/CylinderSupersonic.HiOCase new file mode 100644 index 0000000..dc6a189 --- /dev/null +++ b/TestCases/SupersonicTests/Cylinder/SETUP/CylinderSupersonic.HiOCase @@ -0,0 +1,129 @@ + +! *************************** +!  * Cylinder parameter file * +! *************************** +! + +!-------------------------------------: Geometry + Mesh file: ./MESH/Cylinder.HiOMesh + Boundary file: _this_ + + +!-------------------------------------: Physics + Gas : Air + Reference pressure : 101325.0 + Reference Temperature : 273.15 + Reynolds length : 1.0 + Reynolds number : 0.0 + Prandtl number : 0.72 + Mach number : 2.0 + +!------------------------------------: DG Discretization + Interpolation nodes : Legendre-Gauss + Default polynomial order : 8 + +!------------------------------------: Initialization + Initial condition : Restart + Restart file : ./RESULTS/Cylinder.HiORst + +!------------------------------------: Inviscid discretization + Inviscid discretization : Standard + Inviscid formulation : Green form + Inviscid Riemann solver : Roe + +!------------------------------------: Viscous discretization + Viscous discretization: IP + Interior penalty method: IIPG + Jumps penalty parameter: 1.0 + Gradient jumps penalty parameter: 0.0 + +!------------------------------------: Artificial dissipation + Artificial dissipation (0/1): 1 + Artificial dissipation type: Laplacian + Artificial dissipation indicator: Jumps-based + Artificial dissipation intensity: 0.2 + +!-----------------------------------: Time integration + Integration mode : Steady + Integration scheme : Williamson RK5 + Time step : 1.0e-2 + CFL Number : 0.05 + Simulation time : 0.01 + Number of iterations : 1000000 + Initial time : 0.0 + +!----------------------------------: Output variables + Solution file: ./RESULTS/Cylinder.HiORst + Autosave interval: 10000 + Output interval: 10 + Output file type: Interpolated + Number of representation points: 10 + Save variables: rho_rhou_rhov_rhoe_u_v_p_Mach_vort_rhot_rhout_rhovt_rhoet_muart + +! +! ********************************** +!  * Cylinder boundary mapping file * +! ********************************** +! + +# define p-Refinement + 1: 5 + 2: 5 + 3: 5 + 4: 5 +# end + + +# define zone 1 + Name = Symmetry +! Type = newDirichlet +! Type = Dirichlet + Type = Euler wall +# end + +# define zone 2 + Name = Symmetry + Type = Euler wall +! Type = Dirichlet +# end + +# define zone 3 + Name = Inflow + Type = Dirichlet +# end + +# define zone 4 + Name = Outflow +! Type = Dirichlet + Type = Pressure outlet + Type = Reflective +# end + +# define zone 5 + Name = Wall +! Type = Dirichlet +! Type = Viscous wall + Type = Euler wall + Wall type = Adiabatic +# end + +# define probe 1 + Name = Wake_v + x position = 0.75 + y position = 0.0 + Variable = v +# end + +# define surface monitor 1 + Name = lift + Marker = 5 + Variable = lift + Reference surface = 1.0 +# end + +# define surface monitor 2 + Name = drag + Marker = 5 + Variable = drag + Reference surface = 1.0 +# end diff --git a/TestCases/SupersonicTests/Ramp/MESH/Ramp.HiOMesh b/TestCases/SupersonicTests/Ramp/MESH/Ramp.HiOMesh new file mode 100644 index 0000000..3133ff1 Binary files /dev/null and b/TestCases/SupersonicTests/Ramp/MESH/Ramp.HiOMesh differ diff --git a/TestCases/SupersonicTests/Ramp/SETUP/Ramp.HiOCase b/TestCases/SupersonicTests/Ramp/SETUP/Ramp.HiOCase new file mode 100644 index 0000000..9d7cd45 --- /dev/null +++ b/TestCases/SupersonicTests/Ramp/SETUP/Ramp.HiOCase @@ -0,0 +1,99 @@ +! +! *************************** +!  * Cylinder parameter file * +! *************************** +! + +!-------------------------------------: Geometry + Mesh file : ./MESH/Ramp.HiOMesh + Boundary file : _this_ + + +!-----------------------------------------: Physics + Gas : Air + Reference pressure : 101325.0 + Reference Temperature : 273.15 + Reynolds length : 1.0 + Reynolds number : 0.0 + Prandtl number : 0.72 + Mach number : 2.0 + +!------------------------------------: DG Discretization + Interpolation nodes : Legendre-Gauss + Default polynomial order : 3 + +!------------------------------------: Initialization + Initial condition : Uniform + Restart file : ./RESULTS/Ramp.HiORst + +!------------------------------------: Inviscid discretization + Inviscid discretization : Standard + Inviscid formulation : Green form + Inviscid Riemann solver : Roe + +!------------------------------------: Viscous discretization + Viscous discretization: IP + Interior penalty method: IIPG + Jumps penalty parameter: 1.0 + Gradient jumps penalty parameter: 0.0 + +!------------------------------------: Artificial dissipation + Artificial dissipation (0/1): 1 + Artificial dissipation type: Laplacian + Artificial dissipation intensity: 0.5 + Artificial dissipation indicator: Jumps-based + + +!----------------------------------------: Time integration + Integration mode : Steady + Integration scheme : Williamson RK5 + Time step : 2.0e-0 + CFL Number : 0.05 + Simulation time : 1.0 + Number of iterations : 10000 + Initial time : 0.0 + +!-----------------------------------: Output variables + Solution file : ./RESULTS/Ramp.HiORst + Autosave interval : 10000 + Output file type: Interpolated + Number of representation points: 10 + Output interval : 10 + Save variables : rho_rhou_rhov_rhoe_rhot_rhout_rhovt_rhoet_u_v_p_Mach + +! +! ********************************** +!  * Cylinder boundary mapping file * +! ********************************** +! + + + +# define zone 1 + Name = Ramp + Type = Euler wall +# end + +# define zone 2 + Name = Symmetry + Type = Euler wall +# end + +# define zone 3 + Name = Outflow + Type = newDirichlet +! Type = Pressure outlet +! Outflow type = Partially non reflective +# end + +# define zone 4 + Name = Inflow + Type = newDirichlet +! Type = newDirichlet +! Mode = Specify total pressure +# end + +# define volume monitor 1 + Name = dSnorm + Variable = dSnorm +# end diff --git a/TestCases/SupersonicTests/ShockWave/MESH/ShockWave.HiOMesh b/TestCases/SupersonicTests/ShockWave/MESH/ShockWave.HiOMesh new file mode 100644 index 0000000..bcdbc03 Binary files /dev/null and b/TestCases/SupersonicTests/ShockWave/MESH/ShockWave.HiOMesh differ diff --git a/TestCases/SupersonicTests/ShockWave/RESULTS/InitialConditionP5.HiORst b/TestCases/SupersonicTests/ShockWave/RESULTS/InitialConditionP5.HiORst new file mode 100644 index 0000000..b47246f Binary files /dev/null and b/TestCases/SupersonicTests/ShockWave/RESULTS/InitialConditionP5.HiORst differ diff --git a/TestCases/SupersonicTests/ShockWave/SETUP/ShockWave.HiOCase b/TestCases/SupersonicTests/ShockWave/SETUP/ShockWave.HiOCase new file mode 100644 index 0000000..37cf2fa --- /dev/null +++ b/TestCases/SupersonicTests/ShockWave/SETUP/ShockWave.HiOCase @@ -0,0 +1,95 @@ +! +! ****************************** +!  * Free stream parameter file * +! ****************************** +! + +!-------------------------------------: Geometry + Mesh file : ./MESH/ShockWave.HiOMesh + Boundary file : _this_ + + +!-----------------------------------------: Physics + Gas : Air + Reference pressure : 101325.0 + Reference Temperature : 273.15 + Reynolds length : 1.0 + Reynolds number : 00.0 + Prandtl number : 0.72 + Mach number : 1.5 + +!------------------------------------: DG Discretization + Interpolation nodes : Legendre-Gauss + Default polynomial order : 5 + +!------------------------------------: Initialization + Initial condition : Restart + Restart file : ./RESULTS/InitialConditionP5.HiORst + +!---------------------------------------------: Inviscid discretization + Inviscid discretization : Standard + Inviscid formulation : Green form + Inviscid Riemann solver : Exact + ! Number of integration points : 6 + +!-------------------------------------------: Viscous discretization + Viscous discretization : IP + Interior penalty method : SIPG + Jumps penalty parameter : 10.0 + Gradient jumps penalty parameter : 0.0 + +!------------------------------------: Artificial dissipation + Artificial dissipation (0/1): 1 + Artificial dissipation type: Laplacian + Artificial dissipation indicator: Jumps-based + + +!----------------------------------------: Time integration + Integration mode : Steady + Integration scheme : Williamson RK5 + Time step : 1.0e-4 + CFL Number : 0.01 + Simulation time : 1.0 + Number of iterations : 1000 + Initial time : 0.0 + +!----------------------------------: Output variables + Solution file : ./RESULTS/ShockWave.HiORst + Autosave interval : 1000000 + Output interval : 10 + Output file type : Interpolated + Number of representation points : 12 + Save variables : rho_rhou_rhov_rhoe_rhot_rhout_rhovt_rhoet_u_v_Mach_p + +! +! ************************************* +!  * Free stream boundary mapping file * +! ************************************* +! + + + +# define zone 1 + Name = Symmetry_bottom + Type = Euler wall +# end + +# define zone 2 + Name = Symmetry_top +! Type = Dirichlet + Type = Euler wall +# end + +# define zone 3 + Name = Inflow + Type = Dirichlet + Mach number = 1.5 + pressure = 4.121694915254237e+04 + Temperature = 2.068979544126242e+02 +# end + +# define zone 4 + Name = Outflow + Type = Pressure outlet + Outflow type = Reflective +# end