fix exact rieman solver MU

.github/workflows/CI_sequential 3.yml

@@ -192,9 +192,26 @@ jobs:
           source /opt/intel/oneapi/setvars.sh || true
           ./horses3d.mu Convergence.control
         if: ('!cancelled()' && matrix.compiler == 'ifort')
+#
+# 2) Convergence_cs
+# ----------------------------
+# This test case is skipped with gfortran as the result is slightly different
 
+      - name: Build MultiphaseConvergence_cs
+        working-directory: ./Solver/test/Multiphase/Convergence_cs/SETUP
+        run: |
+          source /opt/intel/oneapi/setvars.sh || true
+          make MODE=${{matrix.mode}} COMPILER=${{matrix.compiler}} COMM=${{matrix.comm}} ENABLE_THREADS=${{matrix.enable_threads}} WITH_MKL=${{matrix.mkl}}
+        if: ('!cancelled()' && matrix.compiler == 'ifort')
+
+      - name: Run MultiphaseConvergence_cs
+        working-directory: ./Solver/test/Multiphase/Convergence_cs
+        run: |
+          source /opt/intel/oneapi/setvars.sh || true
+          ./horses3d.mu Convergence_cs.control
+        if: ('!cancelled()' && matrix.compiler == 'ifort')
 #
-# 2) RisingBubble
+# 3) RisingBubble
 # ----------------------------
 # This test case is skipped with gfortran as the result is slightly different
 
@@ -213,7 +230,7 @@ jobs:
         if: ('!cancelled()' && matrix.compiler == 'ifort')
 
 #
-# 3) Pipe
+# 4) Pipe
 # ----------------------------
 
       - name: Build MultiphasePipe
@@ -231,7 +248,7 @@ jobs:
         if: '!cancelled()'
 
 #
-# 4) Entropy conserving test
+# 5) Entropy conserving test
 # ----------------------------
 
       - name: Build MultiphaseEntropyConservingTest
@@ -249,7 +266,7 @@ jobs:
         if: '!cancelled()'
 
 #
-# 5) RisingBubbleVreman
+# 6) RisingBubbleVreman
 # ----------------------------
 # This test case is skipped with gfortran as the result is slightly different
 
@@ -268,7 +285,7 @@ jobs:
         if: ('!cancelled()')
 
 #
-# 6) Mu AL
+# 7) Mu AL
 # -------------------------------
 
       - name: Build Mu-AL
@@ -286,7 +303,7 @@ jobs:
         if: '!cancelled()'
 
 #
-# 7) Multiphase monopole with acoustics, sponge, and Reinforcement Learning p-adaptation
+# 8) Multiphase monopole with acoustics, sponge, and Reinforcement Learning p-adaptation
 # --------------------------------------------
       - name: Build Multiphase_Monopole_pAdaptationRL
         working-directory: ./Solver/test/Multiphase/Monopole_pAdaptationRL/SETUP
@@ -303,7 +320,7 @@ jobs:
         if: '!cancelled()'
 
 #
-# 8) Snell
+# 9) Snell
 # ----------------------------
       - name: Build MultiphaseSnell
         working-directory: ./Solver/test/Multiphase/Snell/SETUP
@@ -320,7 +337,7 @@ jobs:
         if: ('!cancelled()')
 
 #
-# 9) Mixed RK
+# 10) Mixed RK
 # ----------------------------
       - name: Build MultiphaseMixedRK
         working-directory: ./Solver/test/Multiphase/MixedRK/SETUP

Solver/configure

@@ -88,6 +88,7 @@ TEST_CASES="./Euler/BoxAroundCircle \
             ./IncompressibleNS/ActuatorLineInterpolation \
             ./CahnHilliard/TJokisaari \
             ./Multiphase/Convergence \
+            ./Multiphase/Convergence_cs \
             ./Multiphase/EntropyConservingTest \
             ./Multiphase/RisingBubble \
             ./Multiphase/RisingBubbleVreman \

Solver/src/MultiphaseSolver/SpatialDiscretization.f90

@@ -52,9 +52,12 @@ end subroutine computeBoundaryFluxF
 
       character(len=LINE_LENGTH), parameter  :: viscousDiscretizationKey = "viscous discretization"
       character(len=LINE_LENGTH), parameter     :: CHDiscretizationKey = "cahn-hilliard discretization"
+      character(len=LINE_LENGTH), parameter  :: FLUID1_COMPRESSIBILITY_KEY = "fluid 1 sound speed square (m/s)"
+
 
       real(kind=RP), protected :: IMEX_S0 = 0.0_RP 
       real(kind=RP), protected :: IMEX_K0 = 1.0_RP
+      logical :: use_non_constant_speed_of_sound = .false.
 !
 !     ========      
       CONTAINS 
@@ -188,6 +191,13 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh)
 
             end select
 
+            use_non_constant_speed_of_sound = controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY)
+            if(use_non_constant_speed_of_sound) then
+               write(STD_OUT,'(A)') "  Implementing artificial compressibility with a non-constant speed of sound in each phase"
+            else
+               write(STD_OUT,'(A)') "  Implementing artificial compressibility with a constant speed of sound in each phase"
+            endif
+
             call CHDiscretization % Construct(controlVariables, ELLIPTIC_CH)
             call CHDiscretization % Describe
 
@@ -448,7 +458,7 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
 #endif
 !
 !        -------------------------------------
-!        Get the density in faces and elements
+!        Get the density and invMa2 in faces and elements
 !        -------------------------------------
 !
 !$omp do schedule(runtime)
@@ -459,9 +469,17 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
             if (compute_element) then
 
                mesh % elements(eID) % storage % rho = dimensionless % rho(2) + (dimensionless % rho(1)-dimensionless % rho(2))*mesh % elements(eID) % storage % Q(IMC,:,:,:)
-
                mesh % elements(eID) % storage % rho = min(max(mesh % elements(eID) % storage % rho, dimensionless % rho_min),dimensionless % rho_max)
+
+               if (use_non_constant_speed_of_sound ) then
+                  mesh % elements(eID) % storage % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % elements(eID) % storage % Q(IMC,:,:,:),0.0_RP),1.0_RP)))**2 
+                  mesh % elements(eID) % storage % invMa2 = mesh % elements(eID) % storage % invMa2*mesh % elements(eID) % storage % rho 
+               else
+                  mesh % elements(eID) % storage % invMa2 = dimensionless % invMa2(1)
+               endif
+
             endif 
+
          end do
 !$omp end do nowait
 
@@ -477,6 +495,20 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
 
                mesh % faces(fID) % storage(1) % rho = min(max(mesh % faces(fID) % storage(1) % rho, dimensionless % rho_min),dimensionless % rho_max)
                mesh % faces(fID) % storage(2) % rho = min(max(mesh % faces(fID) % storage(2) % rho, dimensionless % rho_min),dimensionless % rho_max)
+
+
+               if (use_non_constant_speed_of_sound ) then
+                  mesh % faces(fID) % storage(1) % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % faces(fID) % storage(1) % Q(IMC,:,:),0.0_RP),1.0_RP)))**2
+                  mesh % faces(fID) % storage(2) % invMa2 = (sqrt(dimensionless % invMa2(1)/dimensionless % rho(1)) * min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:),0.0_RP),1.0_RP) + sqrt(dimensionless % invMa2(2)/dimensionless % rho(2)) * (1.0_RP - min(max(mesh % faces(fID) % storage(2) % Q(IMC,:,:),0.0_RP),1.0_RP)))**2 
+
+                  mesh % faces(fID) % storage(1) % invMa2 = mesh % faces(fID) % storage(1) % invMa2*mesh % faces(fID) % storage(1) % rho
+                  mesh % faces(fID) % storage(2) % invMa2 = mesh % faces(fID) % storage(2) % invMa2*mesh % faces(fID) % storage(2) % rho
+
+               else
+                  mesh % faces(fID) % storage(1) % invMa2 = dimensionless % invMa2(1)
+                  mesh % faces(fID) % storage(2) % invMa2 = dimensionless % invMa2(2)
+               endif
+
             endif
          end do
 !$omp end do
@@ -511,7 +543,9 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
                associate(e => mesh % elements(eID))
                do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
                   sqrtRho = sqrt(e % storage % rho(i,j,k))
-                  invMa2 = dimensionless % invMa2(1) * min(max(e % storage % Q(IMC,i,j,k),0.0_RP),1.0_RP) + dimensionless % invMa2(2) * (1.0_RP - min(max(e % storage % Q(IMC,i,j,k),0.0_RP),1.0_RP)) 
+
+                  invMa2 = e % storage % invMa2(i,j,k)
+
                   e % storage % QDot(IMC,i,j,k)      = 0.0_RP
                   e % storage % QDot(IMSQRHOU,i,j,k) = -0.5_RP*sqrtRho*(  e % storage % Q(IMSQRHOU,i,j,k)*e % storage % U_x(IGU,i,j,k) & 
                                                                         + e % storage % Q(IMSQRHOV,i,j,k)*e % storage % U_y(IGU,i,j,k) &   
@@ -1094,7 +1128,9 @@ SUBROUTINE computeElementInterfaceFlux_MU(f)
                                   t1     = f % geom % t1(:,i,j), &
                                   t2     = f % geom % t2(:,i,j), &
                                   fL     = inv_fluxL(:,i,j), &
-                                  fR     = inv_fluxR(:,i,j) )
+                                  fR     = inv_fluxR(:,i,j), &
+                                  invMa2L= f % storage(1) % invMa2(i,j), &
+                                  invMa2R= f % storage(2) % invMa2(i,j))
 
 
 !
@@ -1190,7 +1226,9 @@ SUBROUTINE computeMPIFaceFlux_MU(f)
                                   t1     = f % geom % t1(:,i,j), &
                                   t2     = f % geom % t2(:,i,j), &
                                   fL     = inv_fluxL(:,i,j), &
-                                  fR     = inv_fluxR(:,i,j) )
+                                  fR     = inv_fluxR(:,i,j),&
+                                  invMa2L= f % storage(1) % invMa2(i,j), &
+                                  invMa2R= f % storage(2) % invMa2(i,j)) 
 
 
 !
@@ -1293,7 +1331,9 @@ SUBROUTINE computeBoundaryFlux_MU(f, time)
                                t1     = f % geom % t1(:,i,j), &
                                t2     = f % geom % t2(:,i,j), &
                                fL     = inv_fluxL, &
-                               fR     = inv_fluxR )
+                               fR     = inv_fluxR,&
+                               invMa2L= f % storage(1) % invMa2(i,j), &
+                               invMa2R= f % storage(2) % invMa2(i,j)) 
 
             fStar(:,i,j) = (inv_fluxL - visc_flux(:,i,j) ) * f % geom % jacobian(i,j)
          end do   

Solver/src/libs/mesh/StorageClass.f90

@@ -106,6 +106,9 @@ module StorageClass
       real(kind=RP), dimension(:,:,:,:),   allocatable :: mu_z  ! CHE chemical potential z-gradient
       real(kind=RP), dimension(:,:,:,:),   allocatable :: v     ! CHE flow field velocity
       real(kind=RP), dimension(:,:,:,:),   allocatable :: G_CH  ! CHE auxiliary storage
+#endif
+#ifdef MULTIPHASE
+      real(kind=RP),           allocatable :: invMa2(:,:,:)        !Storage for the density artificial compressibility factor
 #endif
       contains
          procedure   :: Assign              => ElementStorage_Assign
@@ -192,6 +195,9 @@ module StorageClass
 #ifdef ACOUSTIC
       real(kind=RP), dimension(:,:,:),     allocatable :: Qbase ! Base flow State vector
 #endif
+#ifdef MULTIPHASE
+      real(kind=RP), dimension(:,:),       allocatable :: invMa2
+#endif
 !
 !     Inviscid Jacobians
 !     ------------------
@@ -858,6 +864,8 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
                allocate(self % RKSteps(k) % hatK(1:NCONS,0:Nx, 0:Ny, 0:Nz))
             enddo
          end if
+
+         allocate ( self % invMa2   (0:Nx,0:Ny,0:Nz) )
 #endif
 !
 !        -----------------
@@ -905,6 +913,9 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
          self % v     = 0.0_RP
 #endif
 
+#ifdef MULTIPHASE
+         self % invMa2     = 0.0_RP
+#endif
          self % first_sensed = huge(1)
          self % prev_sensor = 1.0_RP
          self % sensor = 1.0_RP  ! Activate the sensor by default (first time-step when SC is on)
@@ -1121,6 +1132,11 @@ elemental subroutine ElementStorage_Destruct(self)
          safedeallocate(self % G_CH)
          safedeallocate(self % v)
 #endif
+
+#ifdef MULTIPHASE
+         safedeallocate(self % invMa2)
+#endif
+
          safedeallocate(self % PrevQ)
 
       end subroutine ElementStorage_Destruct
@@ -1390,6 +1406,10 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
 !        -----------------------------------------------------------------------
          interfaceFluxMemorySize = max(interfaceFluxMemorySize, NCOMP*nDIM*product(Nf+1))
 #endif
+
+#ifdef MULTIPHASE
+         allocate( self % invMa2       (0:Nf(1),0:Nf(2)) )
+#endif
 !
 !        Reserve memory for the interface fluxes
 !        ---------------------------------------
@@ -1447,6 +1467,11 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
          self % mu_z  = 0.0_RP
          self % v     = 0.0_RP
 #endif
+
+#ifdef MULTIPHASE
+         self % invMa2    = 0.0_RP
+#endif 
+
          self % constructed = .TRUE.
       end subroutine FaceStorage_Construct
 !
@@ -1539,6 +1564,11 @@ elemental subroutine FaceStorage_Destruct(self)
          safedeallocate(self % mu_z)
          safedeallocate(self % v)
 #endif
+
+#ifdef MULTIPHASE
+         safedeallocate(self % invMa2 )
+#endif
+
          safedeallocate(self % genericInterfaceFluxMemory)
 
          self % Q      => NULL()
@@ -1706,6 +1736,9 @@ elemental subroutine FaceStorage_Assign(to, from)
          to % mu_y = from % mu_y
          to % mu_z = from % mu_z
          to % v = from % v
+#endif
+#ifdef MULTIPHASE
+         to % invMa2 = from % invMa2
 #endif
          call to % PointStorage
       end subroutine FaceStorage_Assign

Solver/src/libs/physics/multiphase/PhysicsStorage_MU.f90

@@ -174,7 +174,7 @@ SUBROUTINE ConstructPhysicsStorage_MU( controlVariables, Lref, timeref, success
 
       if ( controlVariables % ContainsKey(ARTIFICIAL_COMPRESSIBILITY_KEY) .AND.  controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY) .OR. & 
            controlVariables % ContainsKey(ARTIFICIAL_COMPRESSIBILITY_KEY) .AND.  controlVariables % ContainsKey(FLUID2_COMPRESSIBILITY_KEY)) then
-            error stop "Error: Either define a single artificial compressibility or the speed of sound in each fluids. Not both"
+            error stop "Error: Either define a single artificial compressibility or the speed of sound in each fluid, not both"
       endif
 
 
@@ -291,7 +291,7 @@ SUBROUTINE ConstructPhysicsStorage_MU( controlVariables, Lref, timeref, success
          dimensionless_ % invMa2(1) = dimensionless_ % rho(1) * thermodynamics_ % c02(1) / POW2(refValues_ % V)
       endif
 
-      if(controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY)) then 
+      if(controlVariables % ContainsKey(FLUID2_COMPRESSIBILITY_KEY)) then 
          dimensionless_ % Ma2(2) = POW2(refValues_ % V)/(dimensionless_ % rho(2) * thermodynamics_ % c02(2))
          dimensionless_ % invMa2(2) = dimensionless_ % rho(2) * thermodynamics_ % c02(2) / POW2(refValues_ % V)
       endif