Merge pull request #295 from loganoz/293-acoustics-save-lamb-vector-from-flow-simulations-and-read-it-from-acoustics-solver

293 acoustics save lamb vector from flow simulations and read it from acoustics solver

Author: albertjimenezramos

Solver/src/NavierStokesSolver/main.f90

@@ -34,7 +34,7 @@ PROGRAM HORSES3DMainNS
       TYPE( FTValueDictionary)            :: controlVariables
       TYPE( DGSem )                       :: sem
       TYPE( TimeIntegrator_t )            :: timeIntegrator
-      LOGICAL                             :: success, saveGradients, saveSensor, saveLES, saveSource
+      LOGICAL                             :: success, saveGradients, saveSensor, saveLES, saveSource, saveLambVector
 	  logical                             :: generateMonitor = .TRUE.
 	  logical                             :: optimizePartitionLevel=.FALSE.
 	  logical                             :: isMLRK=.FALSE.
@@ -204,7 +204,8 @@ PROGRAM HORSES3DMainNS
          saveSensor       = controlVariables % logicalValueForKey(saveSensorToSolutionKey)
          saveLES          = controlVariables % logicalValueForKey(saveLESToSolutionKey)
          saveSource       = controlVariables % logicalValueForKey(saveSourceToSolutionKey)
-         CALL sem % mesh % SaveSolution(sem % numberOfTimeSteps, timeIntegrator % time, solutionFileName, saveGradients, saveSensor, saveLES, saveSource)
+         saveLambVector   = controlVariables % logicalValueForKey(saveLambVectorToSolutionKey)
+         CALL sem % mesh % SaveSolution(sem % numberOfTimeSteps, timeIntegrator % time, solutionFileName, saveGradients, saveSensor, saveLES, saveSource, saveLambVector)
          if ( sem % particles % active ) then
             call sem % particles % ExportToVTK ( sem % numberOfTimeSteps, sem % monitors % solution_file )
          end if

Solver/src/libs/foundation/Setup.f90

@@ -10,6 +10,7 @@ Module mainKeywordsModule
          CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: restartFileNameKey         = "restart file name"
          CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: solutionFileNameKey        = "solution file name"
          CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveGradientsToSolutionKey = "save gradients with solution"
+         CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveLambVectorToSolutionKey = "save lamb vector" 
          CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveSensorToSolutionKey    = "save sensor with solution"
          CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveLESToSolutionKey       = "save les with solution"
          CHARACTER(LEN=KEYWORD_LENGTH), PARAMETER :: saveSourceToSolutionKey    = "save source with solution"

Solver/src/libs/mesh/HexMesh.f90

@@ -9,7 +9,7 @@ Module MappedGeometryClass
       private
       public MappedGeometry, MappedGeometryFace
       public SetMappingsToCrossProduct
-      public vDot, vCross
+      public vDot, vCross, ComputeVorticity
 !
 !     ---------
 !     Constants
@@ -860,5 +860,19 @@ end subroutine SetMappingsToCrossProduct
 !
 !///////////////////////////////////////////////////////////////////////////////
 !
+      pure subroutine ComputeVorticity(U_x, U_y, U_z, vorticity)
 
+         real(kind=RP), intent(in)  :: U_x(3)
+         real(kind=RP), intent(in)  :: U_y(3)
+         real(kind=RP), intent(in)  :: U_z(3)
+         real(kind=RP), intent(out) :: vorticity(3)
+
+         vorticity(1) = U_y(3) - U_z(2)
+         vorticity(2) = U_z(1) - U_x(3)
+         vorticity(3) = U_x(2) - U_y(1)
+
+      end subroutine ComputeVorticity
+!
+!///////////////////////////////////////////////////////////////////////////////
+!
 END Module MappedGeometryClass

Solver/src/libs/mesh/MappedGeometry.f90

@@ -94,6 +94,8 @@ module StorageClass
 #endif
 #ifdef ACOUSTIC
       real(kind=RP),           allocatable :: Qbase(:,:,:,:)         ! Base flow State vector
+      real(kind=RP),           allocatable :: Lambbase(:,:,:,:)      ! Lamb vector base
+      real(kind=RP),           allocatable :: Lamb_NS(:,:,:,:)          ! Lamb vector from NS
 #endif
 #ifdef CAHNHILLIARD
       real(kind=RP), dimension(:,:,:,:),   allocatable :: c     ! CHE concentration
@@ -195,6 +197,8 @@ module StorageClass
       real(kind=RP), dimension(:,:),       allocatable :: wallNodeDistance ! for BC walls, distance to the first fluid node
 #ifdef ACOUSTIC
       real(kind=RP), dimension(:,:,:),     allocatable :: Qbase ! Base flow State vector
+      real(kind=RP), dimension(:,:,:),     allocatable :: Lambbase ! Lamb vector base
+      real(kind=RP), dimension(:,:,:),     allocatable :: Lamb_NS ! Lamb vector
 #endif
 #ifdef MULTIPHASE
       real(kind=RP), dimension(:,:),       allocatable :: invMa2
@@ -813,6 +817,8 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
 #endif
 #if defined (ACOUSTIC)
          ALLOCATE( self % Qbase  (NCONSB,0:Nx,0:Ny,0:Nz) )
+         ALLOCATE( self % Lambbase  (NDIM,0:Nx,0:Ny,0:Nz) )
+         ALLOCATE( self % Lamb_NS  (NDIM,0:Nx,0:Ny,0:Nz) )
 #endif
          if (computeGradients) then
             ALLOCATE( self % U_xNS (NGRAD,0:Nx,0:Ny,0:Nz) )
@@ -893,6 +899,8 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
 #endif
 #if defined (ACOUSTIC)
          self % Qbase  = 0.0_RP
+         self % Lambbase = 0.0_RP
+         self % Lamb_NS = 0.0_RP
 #endif
          if (computeGradients) then
             self % U_xNS = 0.0_RP
@@ -1002,6 +1010,8 @@ elemental subroutine ElementStorage_Assign(to, from)
 #endif
 #if defined (ACOUSTIC)
          to % Qbase   = from % Qbase
+         to % Lambbase   = from % Lambbase
+         to % Lamb_NS   = from % Lamb_NS
 #endif
 
 #ifndef ACOUSTIC
@@ -1098,6 +1108,8 @@ elemental subroutine ElementStorage_Destruct(self)
 #endif
 #if defined (ACOUSTIC)
          safedeallocate(self % Qbase)
+         safedeallocate(self % Lambbase)
+         safedeallocate(self % Lamb_NS)
 #endif
 
          if (self % computeGradients) then
@@ -1407,6 +1419,8 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
          end if
 #if defined (ACOUSTIC)
          ALLOCATE( self % Qbase (NCONSB,0:Nf(1),0:Nf(2)) )
+         ALLOCATE( self % Lambbase (NDIM,0:Nf(1),0:Nf(2)) )
+         ALLOCATE( self % Lamb_NS (NDIM,0:Nf(1),0:Nf(2)) )
 #endif
 !        Biggest Interface flux memory size is u\vec{n}
 !        ----------------------------------------------
@@ -1470,6 +1484,8 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
          end if
 #if defined (ACOUSTIC)
          self % Qbase    = 0.0_RP
+         self % Lambbase    = 0.0_RP
+         self % Lamb_NS    = 0.0_RP
 #endif
 
          self % rho    = 0.0_RP
@@ -1568,6 +1584,8 @@ elemental subroutine FaceStorage_Destruct(self)
          safedeallocate(self % rho )
 #if defined (ACOUSTIC)
          safedeallocate(self % Qbase )
+         safedeallocate(self % Lambbase )
+         safedeallocate(self % Lamb_NS )
 #endif
 
          self % anJacobian      = .FALSE.
@@ -1732,6 +1750,8 @@ elemental subroutine FaceStorage_Assign(to, from)
          to % QNS = from % QNS
 #if defined (ACOUSTIC)
          to % Qbase = from % Qbase
+         to % Lambbase = from % Lambbase
+         to % Lamb_NS = from % Lamb_NS
 #endif
          if (to % computeGradients) then
             to % U_xNS = from % U_xNS

Solver/src/libs/mesh/StorageClass.f90

@@ -84,7 +84,7 @@ subroutine Monitors_Construct( Monitors, mesh, controlVariables )
          character(len=STR_LEN_MONITORS) :: line
          character(len=STR_LEN_MONITORS) :: solution_file                                            
          logical, save                   :: FirstCall = .TRUE.
-         logical                         :: saveGradients
+         logical                         :: saveGradients, saveLambVector
 		 logical                         :: monitorMem =.FALSE.
 !
 !        Setup the buffer
@@ -147,7 +147,8 @@ subroutine Monitors_Construct( Monitors, mesh, controlVariables )
 
 #if defined(NAVIERSTOKES) || defined(INCNS)
          saveGradients    = controlVariables % logicalValueForKey(saveGradientsToSolutionKey)
-         call Monitors % stats     % Construct(mesh, saveGradients)
+         saveLambVector   = controlVariables % logicalValueForKey(saveLambVectorToSolutionKey)
+         call Monitors % stats     % Construct(mesh, saveGradients, saveLambVector)
 
 
          allocate ( Monitors % surfaceMonitors ( Monitors % no_of_surfaceMonitors )  )

Solver/src/libs/monitors/Monitors.f90

@@ -51,6 +51,7 @@ module StatisticsMonitor
       real(kind=RP)  :: starting_time
       integer        :: no_of_samples
       logical        :: saveGradients
+      logical        :: saveLambVector
       contains
          procedure   :: Construct    => StatisticsMonitor_Construct
          procedure   :: Update       => StatisticsMonitor_Update
@@ -67,23 +68,23 @@ module StatisticsMonitor
 !
 !//////////////////////////////////////////////////////////////////////////////////
 !
-      subroutine StatisticsMonitor_Construct(self, mesh, saveGradients)
+      subroutine StatisticsMonitor_Construct(self, mesh, saveGradients, saveLambVector)
          use ParamfileRegions
          use PhysicsStorage, only: NCONS, NGRAD
-         use HexMeshClass,   only: no_of_stats_variables
          implicit none
          class(StatisticsMonitor_t)    :: self
          class(HexMesh)                :: mesh
-         logical, intent(in)           :: saveGradients
+         logical, intent(in)           :: saveGradients, saveLambVector
          integer, allocatable          :: Nsample, i0, Ndump, Nreset
          real(kind=RP), allocatable    :: t0
          integer                       :: eID
          character(len=LINE_LENGTH)    :: paramFile
 
          NO_OF_VARIABLES = NO_OF_VARIABLES_Sij + NCONS
          self % saveGradients = saveGradients
+         self % saveLambVector = saveLambVector
          if (saveGradients) NO_OF_VARIABLES = NO_OF_VARIABLES + NGRAD * NDIM
-         no_of_stats_variables = NO_OF_VARIABLES
+         if (saveLambVector) NO_OF_VARIABLES = NO_OF_VARIABLES + NDIM
 !
 !        Search for the parameters in the case file
 !        ------------------------------------------
@@ -168,7 +169,11 @@ subroutine StatisticsMonitor_WriteFile(self, mesh, iter, t, solution_file)
          character(len=LINE_LENGTH)    :: fileName
 
          write(fileName,'(A,A)') trim(solution_file),'.stats.hsol'
-         if ( self % state .ne. OFF) call mesh % SaveStatistics(iter, t, trim(fileName), self % saveGradients)
+         if ( self % state .ne. OFF) then
+            call mesh % SaveStatistics(iter, t, trim(fileName), self % saveGradients)
+            write(fileName,'(A,A)') trim(solution_file),'.Lamb.stats.hsol'
+            if (self % saveLambVector) call mesh % SaveLambVectorStatistics(iter, t, trim(fileName), self % saveGradients)
+         end if
 
       end subroutine StatisticsMonitor_WriteFile
 
@@ -225,6 +230,7 @@ end subroutine StatisticsMonitor_Update
 
       subroutine StatisticsMonitor_UpdateValues(self, mesh)
          use PhysicsStorage
+         use MappedGeometryClass, only: vCross, ComputeVorticity
          implicit none
          class(StatisticsMonitor_t)    :: self
          class(HexMesh)              :: mesh
@@ -237,16 +243,23 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
          integer  :: i, j, k
          real(RP) :: ratio, inv_nsamples_plus_1
          real(RP) :: rfactor1, rfactor2
-         integer, dimension(5) :: limits
+         integer, dimension(6) :: limits
+         real(kind=RP), dimension(NDIM) :: vorticity, velocity, LambVector
+         real(kind=RP), pointer, contiguous :: U_x(:), U_y(:), U_z(:)
 
 
 #ifdef NAVIERSTOKES
          !  if gradients are not saved, limits(2) is equal to limits(5), the latter wont be used
             limits(1) = NO_OF_VARIABLES_Sij + IRHO
             limits(2) = NO_OF_VARIABLES_Sij + NCONS
-            limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD
-            limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD
-            limits(5) = NO_OF_VARIABLES
+            limits(3:6) = limits(2)
+            if (self % saveGradients) then
+               limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD ! U_x
+               limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD ! U_y
+               limits(5) = NO_OF_VARIABLES_Sij + NCONS + NDIM * NGRAD ! U_z
+               limits(6) = limits(5)
+            end if
+            if (self % saveLambVector) limits(6) = limits(5) + NDIM
 
             inv_nsamples_plus_1 = 1.0_RP / (self % no_of_samples + 1)
             ratio = self % no_of_samples * inv_nsamples_plus_1
@@ -272,7 +285,18 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
                       data(limits(2)+1:limits(3),i,j,k) = data(limits(2)+1:limits(3),i,j,k) * ratio + e % storage % U_x(:,i,j,k) * inv_nsamples_plus_1
                       data(limits(3)+1:limits(4),i,j,k) = data(limits(3)+1:limits(4),i,j,k) * ratio + e % storage % U_y(:,i,j,k) * inv_nsamples_plus_1
                       data(limits(4)+1:limits(5),i,j,k) = data(limits(4)+1:limits(5),i,j,k) * ratio + e % storage % U_z(:,i,j,k) * inv_nsamples_plus_1
-                  end if 
+                  end if
+                  ! Save Lamb vector: u x w
+                  if (self % saveLambVector) then
+                     ! Attention: qm in APE is -(w x u)'
+                     U_x => e % storage % U_x(:,i,j,k)
+                     U_y => e % storage % U_y(:,i,j,k)
+                     U_z => e % storage % U_z(:,i,j,k)
+                     call ComputeVorticity(U_x, U_y, U_z, vorticity)
+                     velocity = e % storage % Q(IRHOU:IRHOW,i,j,k) / e % storage % Q(IRHO,i,j,k)
+                     call vCross(velocity, vorticity, LambVector)
+                     data(limits(5)+1:limits(6),i,j,k) = data(limits(5)+1:limits(6),i,j,k) * ratio + LambVector * inv_nsamples_plus_1
+                  end if
                end do                  ; end do                   ; end do
 
                end associate
@@ -283,9 +307,14 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
          !  if gradients are not saved, limits(2) is equal to limits(5), the latter wont be used
             limits(1) = NO_OF_VARIABLES_Sij + INSRHO
             limits(2) = NO_OF_VARIABLES_Sij + NCONS
-            limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD
-            limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD
-            limits(5) = NO_OF_VARIABLES
+            limits(3:6) = limits(2)
+            if (self % saveGradients) then
+               limits(3) = NO_OF_VARIABLES_Sij + NCONS + NGRAD ! U_x
+               limits(4) = NO_OF_VARIABLES_Sij + NCONS + 2*NGRAD ! U_y
+               limits(5) = NO_OF_VARIABLES_Sij + NCONS + NDIM * NGRAD ! U_z
+               limits(6) = limits(5)
+            end if
+            if (self % saveLambVector) limits(6) = limits(5) + NDIM
 
             inv_nsamples_plus_1 = 1.0_RP / (self % no_of_samples + 1)
             ratio = self % no_of_samples * inv_nsamples_plus_1
@@ -312,6 +341,17 @@ subroutine StatisticsMonitor_UpdateValues(self, mesh)
                       data(limits(3)+1:limits(4),i,j,k) = data(limits(3)+1:limits(4),i,j,k) * ratio + e % storage % U_y(:,i,j,k) * inv_nsamples_plus_1
                       data(limits(4)+1:limits(5),i,j,k) = data(limits(4)+1:limits(5),i,j,k) * ratio + e % storage % U_z(:,i,j,k) * inv_nsamples_plus_1
                   end if 
+                  ! Save Lamb vector: u x w
+                  if (self % saveLambVector) then
+                     ! Attention: qm in APE is -(w x u)'
+                     U_x => e % storage % U_x(:,i,j,k)
+                     U_y => e % storage % U_y(:,i,j,k)
+                     U_z => e % storage % U_z(:,i,j,k)
+                     call ComputeVorticity(U_x, U_y, U_z, vorticity)
+                     velocity = e % storage % Q(INSRHOU:INSRHOW,i,j,k) / e % storage % Q(INSRHO,i,j,k)
+                     call vCross(velocity, vorticity, LambVector)
+                     data(limits(5)+1:limits(6),i,j,k) = data(limits(5)+1:limits(6),i,j,k) * ratio + LambVector * inv_nsamples_plus_1
+                  end if
                end do                  ; end do                   ; end do
 
                end associate

Solver/src/libs/monitors/StatisticsMonitor.f90

@@ -145,6 +145,7 @@ MODULE PhysicsStorage_NS
       SUBROUTINE ConstructPhysicsStorage_NS( controlVariables, Lref, timeref, success )
       USE FTValueDictionaryClass
       USE Physics_NSKeywordsModule
+      use mainKeywordsModule, only: saveLambVectorToSolutionKey
       use Utilities, only: toLower, almostEqual
 !
 !     ---------
@@ -273,6 +274,14 @@ SUBROUTINE ConstructPhysicsStorage_NS( controlVariables, Lref, timeref, success
          end if
       end if
 
+      ! If Lamb vector is requested, computeGradient should be enabled
+      if (controlVariables % logicalValueForKey(saveLambVectorToSolutionKey)) then
+         if (.not. computeGradients) then
+            print *, trim(COMPUTE_GRADIENTS_KEY), " should be enabled to compute the Lamb vector."
+            error stop
+         end if
+      end if
+
       dimensionless_ % gammaM2 = thermodynamics_ % gamma * POW2( dimensionless_ % Mach )
 !
 !     ********************