Adding sound speed to surface and volume output

doc/options.yaml

@@ -1336,6 +1336,7 @@ surfaceVariables:
       * ``ptloss`` (relative total pressure loss)
       * ``mach`` (mach number)
       * ``blank`` (cell iblank values used for visualization or other post-processing)
+      * ``soundspeed`` (speed of sound)
 
     In addition, the following variables are available for viscous flows:
 
@@ -1374,6 +1375,7 @@ volumeVariables:
       * ``resrhoe`` (total energy residual)
       * ``resturb`` (turbulence residuals)
       * ``blank`` (cell iblank values used for visualization or other post-processing)
+      * ``soundspeed`` (speed of sound)
 
 forcesAsTractions:
   desc: >
@@ -1612,7 +1614,7 @@ sepSensorKsPhi:
   desc: >
     This option is for KS based separation sensors, ``sepSensorKs`` and ``sepSensorKsArea`` cost function.
     This angle is used to compute the allowable flow deviation region with respect to the freestream projected direction on the desired surface.
-    For example, on a wing, the streamlines on the upper surface flows inboard direction . 
+    For example, on a wing, the streamlines on the upper surface flows inboard direction .
     And in some regions, it goes in the outboard directions.
     Therefore, we provide an allowable flow deviation from the projected freestream on the surface to compute the sensor.
     For instance, when the angle is ``90`` degrees, this approximately closes to the verge of separation.

src/inputParam/inputParamRoutines.F90

@@ -677,6 +677,7 @@ subroutine defaultSurfaceOut
         surfWriteVz = .true.
         surfWriteCp = .true.
         surfWriteMach = .true.
+        surfWriteSoundSpeed = .true.
 
         ! Set the values which depend on the equations to be solved.
 
@@ -742,6 +743,7 @@ subroutine defaultVolumeOut
         volWriteResRho = .true.
         volWriteResMom = .false.
         volWriteResRhoe = .false.
+        volWriteSoundSpeed = .false.
 
         ! Set the values which depend on the equations to be solved.
 
@@ -2368,6 +2370,7 @@ subroutine surfaceVariables(variables)
         surfWriteCavitation = .false.
         surfWriteAxisMoment = .false.
         surfWriteGC = .false.
+        surfWriteSoundSpeed = .false.
 
         ! Initialize nVarSpecified to 0. This serves as a test
         ! later on.
@@ -2511,6 +2514,10 @@ subroutine surfaceVariables(variables)
                 surfWriteGC = .True.
                 nVarSpecified = nVarSpecified + 1
 
+            case ("soundspeed")
+                surfWriteSoundSpeed = .True.
+                nVarSpecified = nVarSpecified + 1
+
             case default
                 pos = len_trim(keyword)
                 write (errorMessage, "(3a)") "Unknown surface output &
@@ -2587,6 +2594,7 @@ subroutine volumeVariables(variables)
         volWriteGC = .false.
         volWriteStatus = .false.
         volWriteIntermittency = .false.
+        volWriteSoundSpeed = .false.
 
         ! Initialize nVarSpecified to 0. This serves as a test
         ! later on.
@@ -2738,6 +2746,10 @@ subroutine volumeVariables(variables)
                 volWriteIntermittency = .true.
                 nVarSpecified = nVarSpecified + 1
 
+            case ("soundspeed")
+                volWriteSoundSpeed = .true.
+                nVarSpecified = nVarSpecified + 1
+
             case default
                 pos = len_trim(keyword)
                 write (errorMessage, "(3a)") "Unknown extra volume output &

src/modules/extraOutput.f90

@@ -22,6 +22,7 @@ module extraOutput
     logical :: surfWriteBlank, surfWriteSepSensor, surfWriteSepSensorKs, surfWriteSepSensorKsArea
     logical :: surfWriteCavitation, surfWriteGC, surfWriteAxisMoment
     logical :: surfWriteForceInDragDir, surfWriteForceInLiftDir
+    logical :: surfWriteSoundSpeed
 !
 !       The logical variables, which define the extra volume variables
 !       to be written.
@@ -37,6 +38,7 @@ module extraOutput
     logical :: volWriteResRhoE, volWriteResTurb, volWriteBlank
     logical :: volWriteShock, volWriteFilteredShock, volWriteGC, volWriteStatus
     logical :: volWriteIntermittency
+    logical :: volWriteSoundSpeed
 !
 !       The logical variables, which define the isosurface variables
 !       to be written.

src/output/outputMod.F90

@@ -171,6 +171,7 @@ subroutine numberOfSurfSolVariables(nSolVar)
         if (surfWriteSepSensorKsArea) nSolVar = nSolVar + 1
         if (surfWriteCavitation) nSolVar = nSolVar + 1
         if (surfWriteGC) nSolVar = nSolVar + 1
+        if (surfWriteSoundSpeed) nSolVar = nSolVar + 1
 
     end subroutine numberOfSurfSolVariables
 
@@ -225,6 +226,7 @@ subroutine numberOfVolSolVariables(nVolSolvar, nVolDiscrVar)
         if (volWriteGC) nVolSolvar = nVolSolvar + 1
         if (volWriteStatus) nVolSolvar = nVolSolvar + 1
         if (volWriteIntermittency) nVolDiscrVar = nVolDiscrVar + 1
+        if (volWriteSoundSpeed) nVolDiscrVar = nVolDiscrVar + 1
 
         ! Check the discrete variables.
 
@@ -571,6 +573,11 @@ subroutine volSolNames(solNames)
             solNames(nn) = cgnsIntermittency
         end if
 
+        if (volWriteSoundSpeed) then
+            nn = nn + 1
+            solNames(nn) = cgnsSoundSpeed
+        end if
+
     end subroutine volSolNames
 
     subroutine surfSolNames(solNames)
@@ -738,6 +745,11 @@ subroutine surfSolNames(solNames)
             solNames(nn) = cgnsGC
         end if
 
+        if (surfWriteSoundSpeed) then
+            nn = nn + 1
+            solNames(nn) = cgnsSoundSpeed
+        end if
+
     end subroutine surfSolNames
 
     subroutine storeSolInBuffer(buffer, copyInBuffer, solName, &
@@ -1370,6 +1382,15 @@ subroutine storeSolInBuffer(buffer, copyInBuffer, solName, &
                 end do
             end do
 
+        case (cgnsSoundSpeed)
+            do k = kBeg, kEnd
+                do j = jBeg, jEnd
+                    do i = iBeg, iEnd
+                        wIO(i, j, k, 1) = sqrt(real(aa(i, j, k)))
+                    end do
+                end do
+            end do
+
         case default
             call terminate("storeSolInBuffer", &
                            "This should not happen")
@@ -1458,6 +1479,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
         real(kind=realType), dimension(:, :), pointer :: pp1, pp2
 
         real(kind=realType), dimension(:, :), pointer :: gamma1, gamma2
+        real(kind=realType), dimension(:, :), pointer :: aa1, aa2
         real(kind=realType), dimension(:, :), pointer :: rlv1, rlv2
         real(kind=realType), dimension(:, :), pointer :: dd2Wall
 
@@ -1571,6 +1593,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             pp1 => p(1, 1:, 1:); pp2 => p(2, 1:, 1:)
             gamma1 => gamma(1, 1:, 1:); gamma2 => gamma(2, 1:, 1:)
+            aa1 => aa(1, 1:, 1:); aa2 => aa(2, 1:, 1:)
 
             if (blockIsMoving) then
                 ss1 => s(1, 1:, 1:, :); ss2 => s(2, 1:, 1:, :)
@@ -1604,6 +1627,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             pp1 => p(ie, 1:, 1:); pp2 => p(il, 1:, 1:)
             gamma1 => gamma(ie, 1:, 1:); gamma2 => gamma(il, 1:, 1:)
+            aa1 => aa(ie, 1:, 1:); aa2 => aa(il, 1:, 1:)
 
             if (blockIsMoving) then
                 ss1 => s(ie - 1, 1:, 1:, :); ss2 => s(ie, 1:, 1:, :)
@@ -1638,6 +1662,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             pp1 => p(1:, 1, 1:); pp2 => p(1:, 2, 1:)
             gamma1 => gamma(1:, 1, 1:); gamma2 => gamma(1:, 2, 1:)
+            aa1 => aa(1:, 1, 1:); aa2 => aa(1:, 2, 1:)
 
             if (blockIsMoving) then
                 ss1 => s(1:, 1, 1:, :); ss2 => s(1:, 2, 1:, :)
@@ -1672,6 +1697,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             pp1 => p(1:, je, 1:); pp2 => p(1:, jl, 1:)
             gamma1 => gamma(1:, je, 1:); gamma2 => gamma(1:, jl, 1:)
+            aa1 => aa(1:, je, 1:); aa2 => aa(1:, jl, 1:)
 
             if (blockIsMoving) then
                 ss1 => s(1:, je - 1, 1:, :); ss2 => s(1:, je, 1:, :)
@@ -1706,6 +1732,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             pp1 => p(1:, 1:, 1); pp2 => p(1:, 1:, 2)
             gamma1 => gamma(1:, 1:, 1); gamma2 => gamma(1:, 1:, 2)
+            aa1 => aa(1:, 1:, 1); aa2 => aa(1:, 1:, 2)
 
             if (blockIsMoving) then
                 ss1 => s(1:, 1:, 1, :); ss2 => s(1:, 1:, 2, :)
@@ -1740,6 +1767,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             pp1 => p(1:, 1:, ke); pp2 => p(1:, 1:, kl)
             gamma1 => gamma(1:, 1:, ke); gamma2 => gamma(1:, 1:, kl)
+            aa1 => aa(1:, 1:, ke); aa2 => aa(1:, 1:, kl)
 
             if (blockIsMoving) then
                 ss1 => s(1:, 1:, ke - 1, :); ss2 => s(1:, 1:, ke, :)
@@ -2371,6 +2399,17 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
                 !print*, sensor
             end do
         end do
+
+        case (cgnsSoundSpeed)
+        do j = rangeFace(2, 1), rangeFace(2, 2)
+            do i = rangeFace(1, 1), rangeFace(1, 2)
+                nn=nn+1
+                ! Compute face sound speed
+                buffer(nn) = sqrt(half * (aa1(i, j) + aa2(i, j)))
+
+            end do
+        end do
+
         end select varName
 
     contains