Merge pull request #291 from loganoz/omarinoDev

Add speed of sound squared to Qbase in CAA solver

Author: oscarmarino

Solver/src/libs/discretization/DGSEMClass.f90

@@ -144,7 +144,9 @@ SUBROUTINE ConstructDGSem( self, meshFileName_, controlVariables, &
       logical                     :: useWeightsPartition                ! Partitioning mesh using DOF of elements as weights
       logical                     :: genMonitor
 	  logical                     :: isReconstruct=.false.
-	  real(kind=RP)               :: QbaseUniform(1:NCONS)
+#if defined(ACOUSTIC)
+	  real(kind=RP)               :: QbaseUniform(1:NCONSB)
+#endif
       character(len=*), parameter :: TWOD_OFFSET_DIR_KEY = "2d mesh offset direction"
       procedure(UserDefinedInitialCondition_f) :: UserDefinedInitialCondition
 #if (!defined(NAVIERSTOKES))
@@ -400,16 +402,16 @@ SUBROUTINE ConstructDGSem( self, meshFileName_, controlVariables, &
 !
 #if defined(ACOUSTIC)
       ! start by default with no flow conditions
-      QbaseUniform = [1.0_RP,0.0_RP,0.0_RP,0.0_RP,1.0_RP/(dimensionless % gammaM2)]
+      QbaseUniform = [1.0_RP,0.0_RP,0.0_RP,0.0_RP,1.0_RP/(dimensionless % gammaM2),1.0_RP/POW2(dimensionless%Mach)]
       call self % mesh % SetUniformBaseFlow(QbaseUniform)
-      call self % mesh % ProlongBaseSolutionToFaces(NCONS)
+      call self % mesh % ProlongBaseSolutionToFaces(NCONSB)
 #ifdef _HAS_MPI_
 !$omp single
-      call self % mesh % UpdateMPIFacesBaseSolution(NCONS)
+      call self % mesh % UpdateMPIFacesBaseSolution(NCONSB)
       ! not efficient, but only done once
       ! we can wait for the communication with more computation in between, but will need to be in a different subroutine
       call mpi_barrier(MPI_COMM_WORLD, ierr)
-      call self % mesh % GatherMPIFacesBaseSolution(NCONS)
+      call self % mesh % GatherMPIFacesBaseSolution(NCONSB)
 !$omp end single
 #endif
 !

Solver/src/libs/discretization/HyperbolicDiscretizationClass.f90

@@ -36,7 +36,11 @@ pure subroutine HyperbolicFlux0D_f(Q, F, rho_, Qbase)
          real(kind=RP), intent(in)           :: Q(1:NCONS     )
          real(kind=RP), intent(out)          :: F(1:NCONS, 1:NDIM)
          real(kind=RP), intent(in), optional :: rho_
+#if defined(ACOUSTIC)
+         real(kind=RP), intent(in), optional :: Qbase(1:NCONSB)
+#else
          real(kind=RP), intent(in), optional :: Qbase(1:NCONS)
+#endif
       end subroutine HyperbolicFlux0D_f
    end interface
 !

Solver/src/libs/mesh/HexElementClass.f90

@@ -952,6 +952,7 @@ subroutine HexElement_ProlongBaseSolutionToFaces(self, nEqn, fFR, fBK, fBOT, fR,
          spAxi   => NodalStorage(N(1))
          spAeta  => NodalStorage(N(2))
          spAzeta => NodalStorage(N(3))
+
 !
 !        *************************
 !        Prolong solution to faces

Solver/src/libs/mesh/HexMesh.f90

@@ -1514,7 +1514,7 @@ subroutine HexMesh_UpdateMPIFacesBaseSolution(self, nEqn)
          do k = 1, nShared
 			domain = MPIfaces % listDomain(k)
 			if (MPIfaces % faces(domain) % no_of_faces > 0) then
-				call MPIfaces % faces(domain) % RecvQ(domain, nEqn, all_reqs(k))
+				call MPIfaces % faces(domain) % RecvQBase(domain, nEqn, all_reqs(k))
 			end if 
          end do
 !
@@ -1542,7 +1542,7 @@ subroutine HexMesh_UpdateMPIFacesBaseSolution(self, nEqn)
                thisSide = self % MPIfaces % faces(domain) % elementSide(mpifID)
                associate(f => self % faces(fID))
                do j = 0, f % Nf(2)  ; do i = 0, f % Nf(1)
-                  self % MPIfaces % faces(domain) % Qsend(counter:counter+nEqn-1) = f % storage(thisSide) % Qbase(:,i,j)
+                  self % MPIfaces % faces(domain) % QBaseSend(counter:counter+nEqn-1) = f % storage(thisSide) % Qbase(:,i,j)
                   counter = counter + nEqn
                end do               ; end do
                end associate
@@ -1552,7 +1552,7 @@ subroutine HexMesh_UpdateMPIFacesBaseSolution(self, nEqn)
 !           Send solution
 !           -------------
 !
-            call MPIfaces % faces(domain) % SendQ(domain, nEqn, all_reqs(idx_send))
+            call MPIfaces % faces(domain) % SendQBase(domain, nEqn, all_reqs(idx_send))
 			idx_send = idx_send + 1
          end do
 !
@@ -1743,7 +1743,7 @@ subroutine HexMesh_GatherMPIFacesBaseSolution(self, nEqn)
                thisSide = self % MPIfaces % faces(domain) % elementSide(mpifID)
                associate(f => self % faces(fID))
                do j = 0, f % Nf(2)  ; do i = 0, f % Nf(1)
-                  f % storage(otherSide(thisSide)) % Qbase(:,i,j) = self % MPIfaces % faces(domain) % Qrecv(counter:counter+nEqn-1)
+                  f % storage(otherSide(thisSide)) % Qbase(:,i,j) = self % MPIfaces % faces(domain) % QBaseRecv(counter:counter+nEqn-1)
                   counter = counter + nEqn
                end do               ; end do
                end associate
@@ -2544,7 +2544,7 @@ subroutine HexMesh_SetConnectivitiesAndLinkFaces(self,nodes,facesList)
 #elif defined(MULTIPHASE)
             call ConstructMPIFacesStorage(self % MPIfaces, NCONS, NCONS, MPI_NDOFS)
 #elif defined(ACOUSTIC)
-            call ConstructMPIFacesStorage(self % MPIfaces, NCONS, NCONS, MPI_NDOFS)
+            call ConstructMPIFacesStorage(self % MPIfaces, NCONS, NCONS, MPI_NDOFS, NCONSB_in=NCONSB)
 #endif
 
 #endif
@@ -4041,15 +4041,15 @@ END SUBROUTINE HexMesh_LoadSolution
     Subroutine HexMesh_SetUniformBaseFlow(self,Q_in)
         Implicit None
          CLASS(HexMesh)                  :: self
-         real(kind=RP), dimension(1:NCONS), intent(in)  :: Q_in
+         real(kind=RP), dimension(1:NCONSB), intent(in)  :: Q_in
 !
 !        ---------------
 !        Local variables
 !        ---------------
          INTEGER                        :: eID, eq
 
          do eID = 1, size(self % elements)
-            do eq = 1,NCONS
+            do eq = 1,NCONSB
                 self % elements(eID) % storage % Qbase(eq,:,:,:) = Q_in(eq)
             end do
          end do

Solver/src/libs/mesh/StorageClass.f90

@@ -812,7 +812,7 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
          ALLOCATE( self % S_SA  (NCONS,0:Nx,0:Ny,0:Nz) )
 #endif
 #if defined (ACOUSTIC)
-         ALLOCATE( self % Qbase  (NCONS,0:Nx,0:Ny,0:Nz) )
+         ALLOCATE( self % Qbase  (NCONSB,0:Nx,0:Ny,0:Nz) )
 #endif
          if (computeGradients) then
             ALLOCATE( self % U_xNS (NGRAD,0:Nx,0:Ny,0:Nz) )
@@ -1406,7 +1406,7 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
              ALLOCATE( self % QdotNS   (NCONS,0:Nf(1),0:Nf(2)) )
          end if
 #if defined (ACOUSTIC)
-         ALLOCATE( self % Qbase (NCONS,0:Nf(1),0:Nf(2)) )
+         ALLOCATE( self % Qbase (NCONSB,0:Nf(1),0:Nf(2)) )
 #endif
 !        Biggest Interface flux memory size is u\vec{n}
 !        ----------------------------------------------