diff --git a/.github/workflows/CI_sequential 3.yml b/.github/workflows/CI_sequential 3.yml index d55a3ff454..87db96eb59 100644 --- a/.github/workflows/CI_sequential 3.yml +++ b/.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 diff --git a/Solver/configure b/Solver/configure index 7636636919..c74e800d47 100755 --- a/Solver/configure +++ b/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 \ diff --git a/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 b/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 index ba2b375ab2..4d861a24b2 100644 --- a/Solver/src/MultiphaseSolver/SpatialDiscretization.f90 +++ b/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 diff --git a/Solver/src/libs/mesh/StorageClass.f90 b/Solver/src/libs/mesh/StorageClass.f90 index bd6a1a2b4a..9e79b1db12 100644 --- a/Solver/src/libs/mesh/StorageClass.f90 +++ b/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 diff --git a/Solver/src/libs/physics/multiphase/PhysicsStorage_MU.f90 b/Solver/src/libs/physics/multiphase/PhysicsStorage_MU.f90 index 5a1f8f271d..9aacadb660 100644 --- a/Solver/src/libs/physics/multiphase/PhysicsStorage_MU.f90 +++ b/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 diff --git a/Solver/src/libs/physics/multiphase/RiemannSolvers_MU.f90 b/Solver/src/libs/physics/multiphase/RiemannSolvers_MU.f90 index e11ffc0bec..046887483d 100644 --- a/Solver/src/libs/physics/multiphase/RiemannSolvers_MU.f90 +++ b/Solver/src/libs/physics/multiphase/RiemannSolvers_MU.f90 @@ -33,7 +33,7 @@ module RiemannSolvers_MU public RiemannSolver, ExactRiemannSolver abstract interface - subroutine RiemannSolverFCN(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, fL,fR) + subroutine RiemannSolverFCN(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, fL,fR, invMa2L, invMa2R) use SMConstants use PhysicsStorage_MU real(kind=RP), intent(in) :: QLeft(1:NCONS) @@ -47,6 +47,8 @@ subroutine RiemannSolverFCN(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, f real(kind=RP), intent(in) :: t2(1:NDIM) real(kind=RP), intent(out) :: fL(1:NCONS) real(kind=RP), intent(out) :: fR(1:NCONS) + real(kind=RP), intent(in) :: invMa2L + real(kind=RP), intent(in) :: invMa2R end subroutine RiemannSolverFCN end interface @@ -138,7 +140,7 @@ end subroutine DescribeRiemannSolver ! !/////////////////////////////////////////////////////////////////////////////////////////// ! - subroutine CentralRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, fL, fR) + subroutine CentralRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, fL, fR, invMa2L, invMa2R) implicit none real(kind=RP), intent(in) :: QLeft(1:NCONS) real(kind=RP), intent(in) :: QRight(1:NCONS) @@ -149,13 +151,15 @@ subroutine CentralRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t real(kind=RP), intent(in) :: nHat(1:NDIM), t1(NDIM), t2(NDIM) real(kind=RP), intent(out) :: fL(1:NCONS) real(kind=RP), intent(out) :: fR(1:NCONS) + real(kind=RP), intent(in) :: invMa2L + real(kind=RP), intent(in) :: invMa2R ! ! --------------- ! Local variables ! --------------- ! - real(kind=RP) :: cL, uL, vL, wL, pL, invSqrtRhoL, invMa2L - real(kind=RP) :: cR, uR, vR, wR, pR, invSqrtRhoR, invMa2R + real(kind=RP) :: cL, uL, vL, wL, pL, invSqrtRhoL + real(kind=RP) :: cR, uR, vR, wR, pR, invSqrtRhoR ! @@ -190,9 +194,6 @@ subroutine CentralRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t fR(IMSQRHOW) = 0.5_RP*rhoR*uR*wR fR(IMP) = 0.0_RP -! Get the left and right face inv Mach^2 - invMa2L = dimensionless % invMa2(1) * min(max(cL,0.0_RP),1.0_RP) + dimensionless % invMa2(2) * (1.0_RP - min(max(cL,0.0_RP),1.0_RP)) - invMa2R = dimensionless % invMa2(1) * min(max(cR,0.0_RP),1.0_RP) + dimensionless % invMa2(2) * (1.0_RP - min(max(cR,0.0_RP),1.0_RP)) ! ! Perform the average and rotation ! -------------------------------- @@ -216,7 +217,7 @@ subroutine CentralRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t end subroutine CentralRiemannSolver - subroutine ExactRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, fL, fR) + subroutine ExactRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, fL, fR, invMa2L, invMa2R) implicit none real(kind=RP), intent(in) :: QLeft(1:NCONS) real(kind=RP), intent(in) :: QRight(1:NCONS) @@ -227,15 +228,17 @@ subroutine ExactRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, real(kind=RP), intent(in) :: nHat(1:NDIM), t1(NDIM), t2(NDIM) real(kind=RP), intent(out) :: fL(1:NCONS) real(kind=RP), intent(out) :: fR(1:NCONS) + real(kind=RP), intent(in) :: invMa2L + real(kind=RP), intent(in) :: invMa2R ! ! --------------- ! Local variables ! --------------- ! - real(kind=RP) :: cL,uL, vL, wL, pL, invRhoL, invSqrtRhoL, lambdaMinusL, lambdaPlusL, invMa2L - real(kind=RP) :: cR,uR, vR, wR, pR, invRhoR, invSqrtRhoR, lambdaMinusR, lambdaPlusR, invMa2R + real(kind=RP) :: cL,uL, vL, wL, pL, invRhoL, invSqrtRhoL, lambdaMinusL, lambdaPlusL + real(kind=RP) :: cR,uR, vR, wR, pR, invRhoR, invSqrtRhoR, lambdaMinusR, lambdaPlusR real(kind=RP) :: rhoStarL, rhoStarR, uStar, pStar, rhoStar, vStar, wStar, cuStar, halfRhouStar - real(kind=RP) :: QLRot(NCONS), QRRot(NCONS) + real(kind=RP) :: QLRot(NCONS), QRRot(NCONS), invMa2Avg, halfRhouL, halfRhouR real(kind=RP) :: lambda_mu = 0.0_RP ! @@ -257,9 +260,7 @@ subroutine ExactRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, wR = invSqrtRhoR * (QRight(IMSQRHOU) * t2(1) + QRight(IMSQRHOV) * t2(2) + QRight(IMSQRHOW) * t2(3)) pR = QRight(IMP) -! Get the left and right face inv Mach^2 - invMa2L = dimensionless % invMa2(1) * min(max(cL,0.0_RP),1.0_RP) + dimensionless % invMa2(2) * (1.0_RP - min(max(cL,0.0_RP),1.0_RP)) - invMa2R = dimensionless % invMa2(1) * min(max(cR,0.0_RP),1.0_RP) + dimensionless % invMa2(2) * (1.0_RP - min(max(cR,0.0_RP),1.0_RP)) + invMa2Avg = 0.5_RP * (invMa2L/rhoL + invMa2R/rhoR) ! This is c^2 (normalized with Vref) ! ! Compute the Star Region ! ----------------------- @@ -288,15 +289,17 @@ subroutine ExactRiemannSolver(QLeft, QRight, rhoL, rhoR, muL, muR, nHat, t1, t2, cuStar = 0.5_RP*(cL*uL + cR*uR) halfRhouStar = 0.5_RP*rhoStar*uStar + halfRhouL = 0.5_RP*rhoL*uL + halfRhouR = 0.5_RP*rhoR*uR ! ! - Add first the common (conservative) part !fL = [cuStar+lambda_mu*(muL-muR), rhoStar*uStar*uStar + pStar, rhoStar*uStar*vStar, rhoStar*uStar*wStar, dimensionless % invMa2 * uStar] - fL = [cuStar+lambda_mu*(muL-muR), rhoStar*uStar*uStar + pStar, rhoStar*uStar*vStar, rhoStar*uStar*wStar, 0.5*(invMa2L+invMa2R) * uStar] + fL = [cuStar+lambda_mu*(muL-muR), halfRhouStar*uStar + pStar, halfRhouStar*vStar, halfRhouStar*wStar, 0.0_RP] ! 0.5*(invMa2L+invMa2R) * uStar fR = fL ! ! - Add the non--conservative part - fL = fL + [0.0_RP, cL*0.5_RP*(muR-muL)-halfRhouStar*uL,-halfRhouStar*vL, -halfRhouStar*wL, -invMa2L*uL] - fR = fR + [0.0_RP, cR*0.5_RP*(muL-muR)-halfRhouStar*uR,-halfRhouStar*vR, -halfRhouStar*wR, -invMa2R*uR] + fL = fL + [0.0_RP, cL*0.5_RP*(muR-muL)+ 0.5_RP*halfRhouStar*(uStar-uL),0.5_RP*halfRhouStar*(vStar-vL), 0.5_RP*halfRhouStar*(wStar-wL), (invMa2L)*(uStar-uL)] + fR = fR + [0.0_RP, cR*0.5_RP*(muL-muR)+ 0.5_RP*halfRhouStar*(uStar-uR),0.5_RP*halfRhouStar*(vStar-vR), 0.5_RP*halfRhouStar*(wStar-wR), (invMa2R)*(uStar-uR)] ! ! ************************************************ ! Return momentum equations to the cartesian frame diff --git a/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 index e0ae61ff97..5aaf0d7ef1 100644 --- a/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/ActuatorLineInterpolation/SETUP/ProblemFile.f90 @@ -591,17 +591,17 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & ! -------------------------------------------------- ! #if defined(MULTIPHASE) - + real(kind=RP), parameter :: residuals(5) = [ 0.0000000000000000E+00_RP,& - 6.6814255992417415E+01_RP,& - 2.7001087612424890E-01_RP,& - 2.9121154361911261E-01_RP,& - 2.2437631252185120E+03_RP] + 6.6848885010291738E+01_RP,& + 0.2700344444108091E+00_RP,& + 0.2914829241166736E+00_RP,& + 2.2442342918587788E+03_RP] real(kind=RP), parameter :: source(5) = [ 0.0000000000000000E+00_RP,& - -5.9825590522814762E-03_RP,& - 1.9265667284469030E-05_RP,& - -1.9887356284781115E-05_RP,& + -5.9825397593181728E-03_RP,& + 0.0000192653654527E+00_RP,& + -1.9887137428830E-05_RP,& 0.0000000000000000E+00_RP] ! CALL initializeSharedAssertionsManager diff --git a/Solver/test/Multiphase/Convergence/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/Convergence/SETUP/ProblemFile.f90 index 73fd8b30ea..7bb0265624 100644 --- a/Solver/test/Multiphase/Convergence/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/Convergence/SETUP/ProblemFile.f90 @@ -518,11 +518,11 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & real(kind=RP), intent(in) :: elapsedTime real(kind=RP), intent(in) :: CPUTime real(kind=RP) :: x, y, z, c, locErr(5), phi, u, v, w, p, rho, rho0 - real(kind=RP), parameter :: saved_errors(5) = [1.6902480371161803E-06_RP, & - 1.9409444210673717E-05_RP, & + real(kind=RP), parameter :: saved_errors(5) = [1.6900393871161803E-06_RP, & + 1.9475704000673717E-05_RP, & 7.4217194511428453E-17_RP, & - 1.9836864099170690E-05_RP, & - 1.7108477257107393E-04_RP] + 1.9905579979170690E-05_RP, & + 1.7071338898107393E-04_RP] integer :: i, j,k, eID CHARACTER(LEN=29) :: testName = "Multiphase convergence" real(kind=RP) :: error(5) diff --git a/Solver/test/Multiphase/Convergence_cs/Convergence_cs.control b/Solver/test/Multiphase/Convergence_cs/Convergence_cs.control new file mode 100644 index 0000000000..1ff5809e34 --- /dev/null +++ b/Solver/test/Multiphase/Convergence_cs/Convergence_cs.control @@ -0,0 +1,92 @@ +! +! ******************* +! Sample control file +! ******************* +! +!-------------------------- Configuration:- + Mesh file name = MESH/4x4.msh + Solution file name = RESULTS/N4.hsol + Save gradients to solution = .false. + Restart = .false. + Restart file name = RESTART_FILE_NAME.hsol + +!-------------------- Physical parameters:- + Fluid 1 Density (kg/m^3) = 1.0 + Fluid 2 Density (kg/m^3) = 2.0 + Fluid 1 Viscosity (Pa.s) = 0.001 + Fluid 2 Viscosity (Pa.s) = 0.001 + ! Gravity direction = [x,y,z] + Velocity direction = [1,0,0] + Compute gradients = .true. + +interface width (m) = 0.70711 +interface tension (N/m) = 6.236E-3 +chemical characteristic time (s) = 1000.0 +!artificial sound speed square (m/s) = 1000.0 +fluid 1 sound speed square (m/s) = 1000.0 +fluid 2 sound speed square (m/s) = 2000.0 + +!------------------------- Discretization:- + Polynomial order = 5 + Discretization nodes = Gauss-Lobatto + Riemann solver = exact ! Standard Roe/Low dissipation Roe + Viscous discretization = BR1 ! IP/BR2 + Cahn-Hilliard discretization = BR1 ! IP/BR2 +#define Jacobian + type = 1 +#end + + +!----------------------- Time integration:- + Time integration = explicit +!Explicit method = RK5 +! CFL = 0.5 +! dCFL = 100.0 +dt = 1.0e-5 + Compute time derivative after timestep = .true. + Number of time steps = 100000 + Output interval = 100 + Convergence tolerance = 1.0e-10 + Simulation type = time-accurate + Final time = 0.1 + Autosave mode = Iteration ! Time + Autosave interval = 1000000 ! 1.0 + +!-------------------- Boundary conditions:- +#define boundary Front + type = periodic + coupled boundary = Back +#end + +#define boundary bottom + type = periodic + coupled boundary = top +#end + +#define boundary top + type = periodic + coupled boundary = bottom +#end + +#define boundary Back + type = periodic + coupled boundary = Front +#end + +#define boundary Left + type = periodic + coupled boundary = Right +#end + +#define boundary Right + type = periodic + coupled boundary = Left +#end + +#define volume monitor 1 + name = entr-rate + variable = entropy rate +#end + + + diff --git a/Solver/test/Multiphase/Convergence_cs/MESH/4x4.msh b/Solver/test/Multiphase/Convergence_cs/MESH/4x4.msh new file mode 100644 index 0000000000..91e17bbe00 --- /dev/null +++ b/Solver/test/Multiphase/Convergence_cs/MESH/4x4.msh @@ -0,0 +1,243 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +7 +2 1 "front" +2 2 "back" +2 3 "top" +2 4 "bottom" +2 5 "left" +2 6 "right" +3 7 "fluid" +$EndPhysicalNames +$Entities +8 12 6 1 +1 -1 -1 -1 0 +2 1 -1 -1 0 +3 -1 1 -1 0 +4 1 1 -1 0 +5 -1 -1 1 0 +6 1 -1 1 0 +10 1 1 1 0 +14 -1 1 1 0 +1 -1 -1 -1 1 -1 -1 0 2 1 -2 +2 -1 1 -1 1 1 -1 0 2 3 -4 +3 -1 -1 -1 -1 1 -1 0 2 1 -3 +4 1 -1 -1 1 1 -1 0 2 2 -4 +7 -1 -1 1 1 -1 1 0 2 5 -6 +8 1 -1 1 1 1 1 0 2 6 -10 +9 -1 1 1 1 1 1 0 2 10 -14 +10 -1 -1 1 -1 1 1 0 2 14 -5 +12 -1 -1 -1 -1 -1 1 0 2 1 -5 +13 1 -1 -1 1 -1 1 0 2 2 -6 +17 1 1 -1 1 1 1 0 2 4 -10 +21 -1 1 -1 -1 1 1 0 2 3 -14 +5 -1 -1 -1 1 1 -1 1 6 4 1 4 -2 -3 +14 -1 -1 -1 1 -1 1 1 4 4 1 13 -7 -12 +18 1 -1 -1 1 1 1 1 2 4 4 17 -8 -13 +22 -1 1 -1 1 1 1 1 3 4 -2 21 -9 -17 +26 -1 -1 -1 -1 1 1 1 1 4 -3 12 -10 -21 +27 -1 -1 1 1 1 1 1 5 4 7 8 9 10 +1 -1 -1 -1 1 1 1 1 7 6 -5 27 14 18 -22 -26 +$EndEntities +$Nodes +23 50 1 50 +0 1 0 1 +1 +-1 -1 -1 +0 2 0 1 +2 +1 -1 -1 +0 3 0 1 +3 +-1 1 -1 +0 4 0 1 +4 +1 1 -1 +0 5 0 1 +5 +-1 -1 1 +0 6 0 1 +6 +1 -1 1 +0 10 0 1 +7 +1 1 1 +0 14 0 1 +8 +-1 1 1 +1 1 0 3 +9 +10 +11 +-0.5 -1 -1 +0 -1 -1 +0.5 -1 -1 +1 2 0 3 +12 +13 +14 +-0.5 1 -1 +0 1 -1 +0.5 1 -1 +1 3 0 3 +15 +16 +17 +-1 -0.5 -1 +-1 0 -1 +-1 0.5 -1 +1 4 0 3 +18 +19 +20 +1 -0.5 -1 +1 0 -1 +1 0.5 -1 +1 7 0 3 +21 +22 +23 +-0.5 -1 1 +0 -1 1 +0.5 -1 1 +1 8 0 3 +24 +25 +26 +1 -0.5 1 +1 0 1 +1 0.5 1 +1 9 0 3 +27 +28 +29 +0.5 1 1 +0 1 1 +-0.5 1 1 +1 10 0 3 +30 +31 +32 +-1 0.5 1 +-1 0 1 +-1 -0.5 1 +2 5 0 9 +33 +34 +35 +36 +37 +38 +39 +40 +41 +-0.5 -0.5 -1 +-0.5 0 -1 +-0.5 0.5 -1 +0 -0.5 -1 +0 0 -1 +0 0.5 -1 +0.5 -0.5 -1 +0.5 0 -1 +0.5 0.5 -1 +2 14 0 0 +2 18 0 0 +2 22 0 0 +2 26 0 0 +2 27 0 9 +42 +43 +44 +45 +46 +47 +48 +49 +50 +-0.5 -0.5 1 +-0.5 0 1 +-0.5 0.5 1 +0 -0.5 1 +0 0 1 +0 0.5 1 +0.5 -0.5 1 +0.5 0 1 +0.5 0.5 1 +3 1 0 0 +$EndNodes +$Elements +7 64 1 64 +2 5 3 16 +1 1 9 33 15 +2 15 33 34 16 +3 16 34 35 17 +4 17 35 12 3 +5 9 10 36 33 +6 33 36 37 34 +7 34 37 38 35 +8 35 38 13 12 +9 10 11 39 36 +10 36 39 40 37 +11 37 40 41 38 +12 38 41 14 13 +13 11 2 18 39 +14 39 18 19 40 +15 40 19 20 41 +16 41 20 4 14 +2 14 3 4 +17 1 9 21 5 +18 9 10 22 21 +19 10 11 23 22 +20 11 2 6 23 +2 18 3 4 +21 2 18 24 6 +22 18 19 25 24 +23 19 20 26 25 +24 20 4 7 26 +2 22 3 4 +25 3 8 29 12 +26 12 29 28 13 +27 13 28 27 14 +28 14 27 7 4 +2 26 3 4 +29 1 5 32 15 +30 15 32 31 16 +31 16 31 30 17 +32 17 30 8 3 +2 27 3 16 +33 5 21 42 32 +34 32 42 43 31 +35 31 43 44 30 +36 30 44 29 8 +37 21 22 45 42 +38 42 45 46 43 +39 43 46 47 44 +40 44 47 28 29 +41 22 23 48 45 +42 45 48 49 46 +43 46 49 50 47 +44 47 50 27 28 +45 23 6 24 48 +46 48 24 25 49 +47 49 25 26 50 +48 50 26 7 27 +3 1 5 16 +49 1 9 33 15 5 21 42 32 +50 15 33 34 16 32 42 43 31 +51 16 34 35 17 31 43 44 30 +52 17 35 12 3 30 44 29 8 +53 9 10 36 33 21 22 45 42 +54 33 36 37 34 42 45 46 43 +55 34 37 38 35 43 46 47 44 +56 35 38 13 12 44 47 28 29 +57 10 11 39 36 22 23 48 45 +58 36 39 40 37 45 48 49 46 +59 37 40 41 38 46 49 50 47 +60 38 41 14 13 47 50 27 28 +61 11 2 18 39 23 6 24 48 +62 39 18 19 40 48 24 25 49 +63 40 19 20 41 49 25 26 50 +64 41 20 4 14 50 26 7 27 +$EndElements diff --git a/Solver/test/Multiphase/Convergence_cs/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/Convergence_cs/SETUP/ProblemFile.f90 new file mode 100644 index 0000000000..114fa7775f --- /dev/null +++ b/Solver/test/Multiphase/Convergence_cs/SETUP/ProblemFile.f90 @@ -0,0 +1,669 @@ +! +!//////////////////////////////////////////////////////////////////////// +! +! The Problem File contains user defined procedures +! that are used to "personalize" i.e. define a specific +! problem to be solved. These procedures include initial conditions, +! exact solutions (e.g. for tests), etc. and allow modifications +! without having to modify the main code. +! +! The procedures, *even if empty* that must be defined are +! +! UserDefinedSetUp +! UserDefinedInitialCondition(mesh) +! UserDefinedPeriodicOperation(mesh) +! UserDefinedFinalize(mesh) +! UserDefinedTermination +! +!//////////////////////////////////////////////////////////////////////// +! +#include "Includes.h" +module ProblemFileFunctions + implicit none + + abstract interface + subroutine UserDefinedStartup_f + end subroutine UserDefinedStartup_f + + SUBROUTINE UserDefinedFinalSetup_f(mesh & +#ifdef FLOW + , thermodynamics_ & + , dimensionless_ & + , refValues_ & +#endif +#ifdef CAHNHILLIARD + , multiphase_ & +#endif + ) + USE HexMeshClass + use FluidData + IMPLICIT NONE + CLASS(HexMesh) :: mesh +#ifdef FLOW + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#endif +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif + END SUBROUTINE UserDefinedFinalSetup_f + + subroutine UserDefinedInitialCondition_f(mesh & +#ifdef FLOW + , thermodynamics_ & + , dimensionless_ & + , refValues_ & +#endif +#ifdef CAHNHILLIARD + , multiphase_ & +#endif + ) + use smconstants + use physicsstorage + use hexmeshclass + use fluiddata + implicit none + class(hexmesh) :: mesh +#ifdef FLOW + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#endif +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif + end subroutine UserDefinedInitialCondition_f +#ifdef FLOW + subroutine UserDefinedState_f(x, t, nHat, Q, thermodynamics_, dimensionless_, refValues_) +! +! ------------------------------------------------- +! Used to define an user defined boundary condition +! ------------------------------------------------- +! + use SMConstants + use PhysicsStorage + use FluidData + implicit none + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: t + real(kind=RP), intent(in) :: nHat(NDIM) + real(kind=RP), intent(inout) :: Q(NCONS) + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ + end subroutine UserDefinedState_f + + subroutine UserDefinedNeumann_f(x, t, nHat, U_x, U_y, U_z) + use SMConstants + use PhysicsStorage + use FluidData + implicit none + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: t + real(kind=RP), intent(in) :: nHat(NDIM) + real(kind=RP), intent(inout) :: U_x(NGRAD) + real(kind=RP), intent(inout) :: U_y(NGRAD) + real(kind=RP), intent(inout) :: U_z(NGRAD) + end subroutine UserDefinedNeumann_f +#endif +! +!//////////////////////////////////////////////////////////////////////// +! + SUBROUTINE UserDefinedPeriodicOperation_f(mesh, time, dt, Monitors) + use SMConstants + USE HexMeshClass + use MonitorsClass + IMPLICIT NONE + CLASS(HexMesh) :: mesh + REAL(KIND=RP) :: time + REAL(KIND=RP) :: dt + type(Monitor_t), intent(in) :: monitors + END SUBROUTINE UserDefinedPeriodicOperation_f +! +!//////////////////////////////////////////////////////////////////////// +! +#ifdef FLOW + subroutine UserDefinedSourceTermNS_f(x, Q, time, S, thermodynamics_, dimensionless_, refValues_& +#ifdef CAHNHILLIARD +,multiphase_ & +#endif +) + use SMConstants + USE HexMeshClass + use FluidData + use PhysicsStorage + IMPLICIT NONE + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: Q(NCONS) + real(kind=RP), intent(in) :: time + real(kind=RP), intent(inout) :: S(NCONS) + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif + end subroutine UserDefinedSourceTermNS_f +#endif +! +!//////////////////////////////////////////////////////////////////////// +! + SUBROUTINE UserDefinedFinalize_f(mesh, time, iter, maxResidual & +#ifdef FLOW + , thermodynamics_ & + , dimensionless_ & + , refValues_ & +#endif +#ifdef CAHNHILLIARD + , multiphase_ & +#endif + , monitors, & + elapsedTime, & + CPUTime ) + use SMConstants + USE HexMeshClass + use FluidData + use MonitorsClass + IMPLICIT NONE + CLASS(HexMesh) :: mesh + REAL(KIND=RP) :: time + integer :: iter + real(kind=RP) :: maxResidual +#ifdef FLOW + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#endif +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif + type(Monitor_t), intent(in) :: monitors + real(kind=RP), intent(in) :: elapsedTime + real(kind=RP), intent(in) :: CPUTime + END SUBROUTINE UserDefinedFinalize_f + + SUBROUTINE UserDefinedTermination_f + implicit none + END SUBROUTINE UserDefinedTermination_f + end interface + +end module ProblemFileFunctions + + SUBROUTINE UserDefinedStartup +! +! -------------------------------- +! Called before any other routines +! -------------------------------- +! + IMPLICIT NONE + END SUBROUTINE UserDefinedStartup +! +!//////////////////////////////////////////////////////////////////////// +! + SUBROUTINE UserDefinedFinalSetup(mesh & +#ifdef FLOW + , thermodynamics_ & + , dimensionless_ & + , refValues_ & +#endif +#ifdef CAHNHILLIARD + , multiphase_ & +#endif + ) +! +! ---------------------------------------------------------------------- +! Called after the mesh is read in to allow mesh related initializations +! or memory allocations. +! ---------------------------------------------------------------------- +! + USE HexMeshClass + use PhysicsStorage + use FluidData + IMPLICIT NONE + CLASS(HexMesh) :: mesh +#ifdef FLOW + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#endif +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif + END SUBROUTINE UserDefinedFinalSetup +! +!//////////////////////////////////////////////////////////////////////// +! + subroutine UserDefinedInitialCondition(mesh & +#ifdef FLOW + , thermodynamics_ & + , dimensionless_ & + , refValues_ & +#endif +#ifdef CAHNHILLIARD + , multiphase_ & +#endif + ) +! +! ------------------------------------------------ +! called to set the initial condition for the flow +! - by default it sets an uniform initial +! condition. +! ------------------------------------------------ +! + use smconstants + use physicsstorage + use hexmeshclass + use fluiddata + implicit none + class(hexmesh) :: mesh +#ifdef FLOW + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#endif +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif +! +! --------------- +! local variables +! --------------- +! + integer :: eid, i, j, k + real(kind=RP) :: c, u, v, w, p, phi + real(kind=RP) :: x, y, z +! +! ------------------------------------------------------ +! Incompressible Navier-Stokes default initial condition +! ------------------------------------------------------ +! +#if defined(MULTIPHASE) + do eID = 1, mesh % no_of_elements + associate( Nx => mesh % elements(eID) % Nxyz(1), & + ny => mesh % elemeNts(eID) % nxyz(2), & + Nz => mesh % elements(eID) % Nxyz(3) ) + do k = 0, Nz; do j = 0, Ny; do i = 0, Nx + x = mesh % elements(eID) % geom % x(IX,i,j,k) + y = mesh % elements(eID) % geom % x(IY,i,j,k) + z = mesh % elements(eID) % geom % x(IZ,i,j,k) + + phi = 0.0_RP + c = 0.5_RP + + u = 0.0_RP + v = 0.0_RP + w = 0.0_RP + + p = 2*sin(PI*x)*sin(PI*y) + + mesh % elements(eID) % storage % q(:,i,j,k) = [c,u,v,w,p] + end do; end do; end do + end associate + end do +#endif + + end subroutine UserDefinedInitialCondition +#ifdef FLOW + subroutine UserDefinedState1(x, t, nHat, Q, thermodynamics_, dimensionless_, refValues_) +! +! ------------------------------------------------- +! Used to define an user defined boundary condition +! ------------------------------------------------- +! + use SMConstants + use PhysicsStorage + use FluidData + implicit none + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: t + real(kind=RP), intent(in) :: nHat(NDIM) + real(kind=RP), intent(inout) :: Q(NCONS) + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ + end subroutine UserDefinedState1 + + subroutine UserDefinedGradVars1(x, t, nHat, Q, U, thermodynamics_, dimensionless_, refValues_) + use SMConstants + use PhysicsStorage + use FluidData + implicit none + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: t + real(kind=RP), intent(in) :: nHat(NDIM) + real(kind=RP), intent(in) :: Q(NCONS) + real(kind=RP), intent(inout) :: U(NGRAD) + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ + end subroutine UserDefinedGradVars1 + + subroutine UserDefinedNeumann1(x, t, nHat, Q, U_x, U_y, U_z, flux, thermodynamics_, dimensionless_, refValues_) + use SMConstants + use PhysicsStorage + use FluidData + implicit none + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: t + real(kind=RP), intent(in) :: nHat(NDIM) + real(kind=RP), intent(in) :: Q(NCONS) + real(kind=RP), intent(in) :: U_x(NGRAD) + real(kind=RP), intent(in) :: U_y(NGRAD) + real(kind=RP), intent(in) :: U_z(NGRAD) + real(kind=RP), intent(inout) :: flux(NCONS) + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ + end subroutine UserDefinedNeumann1 +#endif +! +!//////////////////////////////////////////////////////////////////////// +! + SUBROUTINE UserDefinedPeriodicOperation(mesh, time, dt, Monitors) +! +! ---------------------------------------------------------- +! Called before every time-step to allow periodic operations +! to be performed +! ---------------------------------------------------------- +! + use SMConstants + USE HexMeshClass + use MonitorsClass + IMPLICIT NONE + CLASS(HexMesh) :: mesh + REAL(KIND=RP) :: time + REAL(KIND=RP) :: dt + type(Monitor_t), intent(in) :: monitors + + END SUBROUTINE UserDefinedPeriodicOperation +! +!//////////////////////////////////////////////////////////////////////// +! +#ifdef FLOW + subroutine UserDefinedSourceTermNS(x, Q, time, S, thermodynamics_, dimensionless_, refValues_& +#ifdef CAHNHILLIARD +, multiphase_ & +#endif +) +! +! -------------------------------------------- +! Called to apply source terms to the equation +! -------------------------------------------- +! + use SMConstants + USE HexMeshClass + use PhysicsStorage + use FluidData + IMPLICIT NONE + real(kind=RP), intent(in) :: x(NDIM) + real(kind=RP), intent(in) :: Q(NCONS) + real(kind=RP), intent(in) :: time + real(kind=RP), intent(inout) :: S(NCONS) + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +! +! --------------- +! Local variables +! --------------- +! + real(kind=RP) :: rho0, eta = 0.001_RP + real(kind=RP) :: rho1, rho2, cs1, cs2 + + rho1 = dimensionless_ % rho(2) + rho2 = dimensionless_ % rho(1) + cs1 = sqrt(thermodynamics_ % c02(2)) + cs2 = sqrt(thermodynamics_ % c02(1)) + + S(IMC) = 0.5_RP*cos(time)*cos(PI*x(IX))*cos(PI*x(IY)) & !diff(c,t) + +PI* & + (4.0_RP*sin(time)*(cos(PI*x(IX)))**2*(cos(PI*x(IY)))**2 & + -1.0_RP*sin(time)*(cos(PI*x(IX)))**2 & + -1.0_RP*sin(time)*(cos(PI*x(IY)))**2 & + +2.0_RP*cos(PI*x(IX))*cos(PI*x(IY)) & + )*(sin(time)) & !+diff(c*u,x)+diff(c*v,y) + +(PI**2)*(multiphase_ % M0*multiphase_ % sigma/multiphase_ % eps) & + *( & + 3.0_RP*(PI*multiphase_ % eps)**2 & + +108.0_RP*(sin(time)*sin (PI*x(IX))*sin(PI*x(IY)))**2 & + -72.0_RP*(sin(time)*sin(PI*x(IX)))**2 & + -72.0_RP*(sin(time)*sin(PI*x(IY)))**2 & + +36.0_RP*(sin(time)**2) & + -12.0_RP & + ) & + *sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) !-M0*(diff(μ,x,2)+diff(μ,y,2)) + + S(IMSQRHOU) = 1.0_RP*( & + -1.5_RP*rho1*sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + rho1 & + +1.5_RP*rho2*sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + rho2 & + ) & + *sin(PI*x(IX))*cos(time)*cos(PI*x(IY)) & !sqrt(ρ)*diff(sqrt(ρ)*u,t) + + PI * ( & + - 6.0_RP * rho1 * sin(time) * (cos(PI*x(IX)))**2 * (cos(PI*x(IY)))**3 & + + 2.0_RP * rho1 * sin(time) * (cos(PI*x(IX)))**2 * cos(PI*x(IY)) & + + 1.0_RP * rho1 * sin(time) * cos(PI*x(IY))**3 & + + 4.0_RP * rho1 * cos(PI*x(IX)) * (cos(PI*x(IY)))**2 & + - 1.0_RP * rho1 * cos(PI*x(IX)) & + + 6.0_RP * rho2 * sin(time) * (cos(PI*x(IX)))**2 * (cos(PI*x(IY)))**3 & + - 2.0_RP * rho2 * sin(time) * (cos(PI*x(IX)))**2 * cos(PI*x(IY)) & + - 1.0_RP * rho2 * sin(time) * cos(PI*x(IY))**3 & + + 4.0_RP * rho2 * cos(PI*x(IX)) * (cos(PI*x(IY)))**2 & + - 1.0_RP * rho2 * cos(PI*x(IX)) & + ) * (sin(time))**2 * sin(PI*x(IX)) & !0.5*(diff(ρ*u*u,x)+diff(ρ*u*v,y)) + -PI*( & + rho1*(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) - 1.0_RP) & + -rho2*(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + 1.0_RP) & + ) & + *(sin(time))**2*sin(PI*x(IX))*cos(PI*x(IX))*cos(2*PI*x(IY)) & !0.5*ρ*(u*diff(u,x)+v*diff(u,y)) + -0.5_RP*PI*multiphase_ % sigma / multiphase_ % eps & + *(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + 1.0_RP) & + *(1.5_RP*(PI*multiphase_ % eps)**2 + 18.0_RP*(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)))**2 - 6.0_RP) & + *sin(time)*sin(PI*x(IX))*cos(PI*x(IY)) & !+ c *diff(μ,x) + + 2.0_RP*PI*sin(PI*x(IY))*cos(time)*cos(PI*x(IX)) & !+diff(p,x) + + 8.0_RP*(PI**2) *eta*sin(time)*sin(PI*x(IX))*cos(PI*x(IY)) !- diff(η*(diff(u,x) + diff(u,x)),x) - diff(η*(diff(u,y) + diff(v,x)),y) + + S(IMSQRHOV) = 1.0_RP*( & + -1.5_RP*rho1*sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + rho1 & + +1.5_RP*rho2*sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + rho2 & + ) & + *sin(PI*x(IY))*cos(time)*cos(PI*x(IX)) & !sqrt(ρ)*diff(sqrt(ρ)*v,t) + +PI * ( & + -6.0_RP * rho1 * sin(time) * (cos(PI*x(IX)))**3 * (cos(PI*x(IY)))**2 & + + 1.0_RP * rho1 * sin(time) * (cos(PI*x(IX)))**3 & + + 2.0_RP * rho1 * sin(time) * cos(PI*x(IX)) * (cos(PI*x(IY)))**2 & + + 4.0_RP * rho1 * (cos(PI*x(IX)))**2 * cos(PI*x(IY)) & + - 1.0_RP * rho1 * (cos(PI*x(IY))) & + + 6.0_RP * rho2 * sin(time) * (cos(PI*x(IX)))**3 * (cos(PI*x(IY)))**2 & + - 1.0_RP * rho2 * sin(time) * (cos(PI*x(IX)))**3 & + - 2.0_RP * rho2 * sin(time) * cos(PI*x(IX)) * (cos(PI*x(IY)))**2 & + + 4.0_RP * rho2 * (cos(PI*x(IX)))**2 * cos(PI*x(IY)) & + - 1.0_RP * rho2 * (cos(PI*x(IY))) ) & + * (sin(time))**2 * sin(PI*x(IY)) & !0.5*(diff(ρ*v*u,x)+diff(ρ*v*v,y)) + -PI*( & + rho1*(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) - 1.0_RP) & + -rho2*(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + 1.0_RP) & + ) & + *(sin(time))**2*sin(PI*x(IY))*cos(2.0_RP*PI*x(IX))*cos(PI*x(IY)) & !0.5*ρ*(u*diff(v,x)+v*diff(v,y)) + -0.5_RP*PI*multiphase_ % sigma / multiphase_ % eps & + *(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)) + 1.0_RP) & + *(1.5_RP*(PI*multiphase_ % eps)**2 + 18.0_RP*(sin(time)*cos(PI*x(IX))*cos(PI*x(IY)))**2 - 6.0_RP) & + *sin(time)*sin(PI*x(IY))*cos(PI*x(IX)) & !+ c *diff(μ,y) + + 2.0_RP*PI*sin(PI*x(IX))*cos(time)*cos(PI*x(IY)) & !+diff(p,y) + + 8.0_RP *(PI**2)*eta*sin(time)*sin(PI*x(IY))*cos(PI*x(IX)) !- diff(η*(diff(v,x) + diff(u,y)),x) - diff(η*(diff(v,y) + diff(v,y)),y) + + + S(IMSQRHOW) = 0.0_RP + + S(IMP) = sin(PI*x(IX))*sin(PI*x(IY))*sin(time)*(-2.0_RP) & !diff(p,t) + -0.5_RP * PI * ( & + cs1 * (sin(time) * cos(PI*x(IX)) * cos(PI*x(IY)) - 1.0_RP) & + - cs2 * (sin(time) * cos(PI*x(IX)) * cos(PI*x(IY)) + 1.0_RP) )**2 & + * ( & + rho1 * (sin(time) * cos(PI*x(IX)) * cos(PI*x(IY)) - 1.0_RP) & + - rho2 * (sin(time) * cos(PI*x(IX)) * cos(PI*x(IY)) + 1.0_RP) ) & + * sin(time) * cos(PI*x(IX)) * cos(PI*x(IY)) !+ρ*cs^2*(diff(u,x)+diff(v,y)) +#endif + + end subroutine UserDefinedSourceTermNS +#endif +! +!//////////////////////////////////////////////////////////////////////// +! + SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & +#ifdef FLOW + , thermodynamics_ & + , dimensionless_ & + , refValues_ & +#endif +#ifdef CAHNHILLIARD + , multiphase_ & +#endif + , monitors, & + elapsedTime, & + CPUTime ) +! +! -------------------------------------------------------- +! Called after the solution computed to allow, for example +! error tests to be performed +! -------------------------------------------------------- +! + use SMConstants + USE HexMeshClass + use PhysicsStorage + use FluidData + use MonitorsClass + use FTAssertions + IMPLICIT NONE + CLASS(HexMesh) :: mesh + REAL(KIND=RP) :: time + integer :: iter + real(kind=RP) :: maxResidual +#ifdef FLOW + type(Thermodynamics_t), intent(in) :: thermodynamics_ + type(Dimensionless_t), intent(in) :: dimensionless_ + type(RefValues_t), intent(in) :: refValues_ +#endif +#ifdef CAHNHILLIARD + type(Multiphase_t), intent(in) :: multiphase_ +#endif + type(Monitor_t), intent(in) :: monitors + real(kind=RP), intent(in) :: elapsedTime + real(kind=RP), intent(in) :: CPUTime + real(kind=RP) :: x, y, z, c, locErr(5), phi, u, v, w, p, rho, rho0 + real(kind=RP), parameter :: saved_errors(5) = [4.877120855395489E-6_RP, & + 6.8958937066338505E-05_RP, & + 6.8958937066344671E-05_RP, & + 7.4217194511428453E-17_RP, & + 6.1637135702849238E-04_RP] + integer :: i, j,k, eID + CHARACTER(LEN=29) :: testName = "Multiphase convergence non constant sound speed" + real(kind=RP) :: error(5) + TYPE(FTAssertionsManager), POINTER :: sharedManager + LOGICAL :: success + real(kind=RP), parameter :: w_LGL(0:5) = [0.066666666666667_RP, & + 0.378474956297847_RP, & + 0.554858377035486_RP, & + 0.554858377035486_RP, & + 0.378474956297847_RP, & + 0.066666666666667_RP ] +! +! ********************************* +! Check the L-inf norm of the error +! ********************************* +! +#ifdef MULTIPHASE + rho0 = dimensionless_ % rho(2) +#endif + error = 0.0_RP + do eID = 1, mesh % no_of_elements + associate(e => mesh % elements(eID) ) + do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1) + x = e % geom % x(IX,i,j,k) + y = e % geom % x(IY,i,j,k) + z = e % geom % x(IZ,i,j,k) + + phi = cos(PI*x)*cos(PI*y)*sin(time) + c = 0.5_RP*(1.0_RP + phi) + + rho = c + rho0*(1.0_RP-c) + + u = 2*sin(PI*x)*cos(PI*y)*sin(time) + v = 2*cos(PI*x)*sin(PI*y)*sin(time) + w = 0.0_RP + p = 2*sin(PI*x)*sin(PI*y)*cos(time) + + locErr = e % storage % Q(:,i,j,k) - [c,sqrt(rho)*u,sqrt(rho)*v,sqrt(rho)*w,p] + + error = error + e % geom % jacobian(i,j,k)*locErr**2*w_LGL(i)*w_LGL(j)*2.0_RP + + end do ; end do ; end do + end associate + end do + + error = sqrt(error) + + CALL initializeSharedAssertionsManager + sharedManager => sharedAssertionsManager() + + CALL FTAssertEqual(expectedValue = saved_errors(1)+1.0_RP, & + actualValue = error(1)+1.0_RP, & + tol = 1.d-11, & + msg = "Concentration error") + + + CALL FTAssertEqual(expectedValue = saved_errors(2)+1.0_RP, & + actualValue = error(2)+1.0_RP, & + tol = 1.d-11, & + msg = "X-Momentum error") + + CALL FTAssertEqual(expectedValue = saved_errors(3)+1.0_RP, & + actualValue = error(3)+1.0_RP, & + tol = 1.d-11, & + msg = "Y-Momentum error") + + CALL FTAssertEqual(expectedValue = saved_errors(4)+1.0_RP, & + actualValue = error(4)+1.0_RP, & + tol = 1.d-11, & + msg = "Z-Momentum error") + + CALL FTAssertEqual(expectedValue = saved_errors(5)+1.0_RP, & + actualValue = error(5)+1.0_RP, & + tol = 1.d-11, & + msg = "Pressure error") + + CALL sharedManager % summarizeAssertions(title = testName,iUnit = 6) + + IF ( sharedManager % numberOfAssertionFailures() == 0 ) THEN + WRITE(6,*) testName, " ... Passed" + WRITE(6,*) "This test case has no expected solution yet, only checks the residual after 100 iterations." + ELSE + WRITE(6,*) testName, " ... Failed" + WRITE(6,*) "NOTE: Failure is expected when the max eigenvalue procedure is changed." + WRITE(6,*) " If that is done, re-compute the expected values and modify this procedure" + error stop 99 + END IF + WRITE(6,*) + + CALL finalizeSharedAssertionsManager + CALL detachSharedAssertionsManager + + + + + END SUBROUTINE UserDefinedFinalize +! +!//////////////////////////////////////////////////////////////////////// +! + SUBROUTINE UserDefinedTermination +! +! ----------------------------------------------- +! Called at the the end of the main driver after +! everything else is done. +! ----------------------------------------------- +! + IMPLICIT NONE + END SUBROUTINE UserDefinedTermination + \ No newline at end of file diff --git a/Solver/test/Multiphase/Convergence_cs/convergence_cs.control b/Solver/test/Multiphase/Convergence_cs/convergence_cs.control new file mode 100644 index 0000000000..1ff5809e34 --- /dev/null +++ b/Solver/test/Multiphase/Convergence_cs/convergence_cs.control @@ -0,0 +1,92 @@ +! +! ******************* +! Sample control file +! ******************* +! +!-------------------------- Configuration:- + Mesh file name = MESH/4x4.msh + Solution file name = RESULTS/N4.hsol + Save gradients to solution = .false. + Restart = .false. + Restart file name = RESTART_FILE_NAME.hsol + +!-------------------- Physical parameters:- + Fluid 1 Density (kg/m^3) = 1.0 + Fluid 2 Density (kg/m^3) = 2.0 + Fluid 1 Viscosity (Pa.s) = 0.001 + Fluid 2 Viscosity (Pa.s) = 0.001 + ! Gravity direction = [x,y,z] + Velocity direction = [1,0,0] + Compute gradients = .true. + +interface width (m) = 0.70711 +interface tension (N/m) = 6.236E-3 +chemical characteristic time (s) = 1000.0 +!artificial sound speed square (m/s) = 1000.0 +fluid 1 sound speed square (m/s) = 1000.0 +fluid 2 sound speed square (m/s) = 2000.0 + +!------------------------- Discretization:- + Polynomial order = 5 + Discretization nodes = Gauss-Lobatto + Riemann solver = exact ! Standard Roe/Low dissipation Roe + Viscous discretization = BR1 ! IP/BR2 + Cahn-Hilliard discretization = BR1 ! IP/BR2 +#define Jacobian + type = 1 +#end + + +!----------------------- Time integration:- + Time integration = explicit +!Explicit method = RK5 +! CFL = 0.5 +! dCFL = 100.0 +dt = 1.0e-5 + Compute time derivative after timestep = .true. + Number of time steps = 100000 + Output interval = 100 + Convergence tolerance = 1.0e-10 + Simulation type = time-accurate + Final time = 0.1 + Autosave mode = Iteration ! Time + Autosave interval = 1000000 ! 1.0 + +!-------------------- Boundary conditions:- +#define boundary Front + type = periodic + coupled boundary = Back +#end + +#define boundary bottom + type = periodic + coupled boundary = top +#end + +#define boundary top + type = periodic + coupled boundary = bottom +#end + +#define boundary Back + type = periodic + coupled boundary = Front +#end + +#define boundary Left + type = periodic + coupled boundary = Right +#end + +#define boundary Right + type = periodic + coupled boundary = Left +#end + +#define volume monitor 1 + name = entr-rate + variable = entropy rate +#end + + + diff --git a/Solver/test/Multiphase/MixedRK/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/MixedRK/SETUP/ProblemFile.f90 index b6f8daf5a3..7b2506e0e1 100644 --- a/Solver/test/Multiphase/MixedRK/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/MixedRK/SETUP/ProblemFile.f90 @@ -542,38 +542,38 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & CHARACTER(LEN=29) :: testName = "Multiphase:: MixedRK" TYPE(FTAssertionsManager), POINTER :: sharedManager LOGICAL :: success - real(kind=RP), parameter :: residuals_saved(5) = [2.8665833841075737E-06_RP, & - 1.9361250101172587E-03_RP, & - 5.2814258414213508E-11_RP, & + real(kind=RP), parameter :: residuals_saved(5) = [2.8665833841448300E-06_RP, & + 1.9361250102237703E-03_RP, & + 5.9347820202367154E-11_RP, & 2.3047237097072765E-17_RP, & - 7.1908221779661385E-01_RP] + 7.1908221803287242E-01_RP] CALL initializeSharedAssertionsManager sharedManager => sharedAssertionsManager() - CALL FTAssertEqual(expectedValue = residuals_saved(1)+1.0_RP, & - actualValue = monitors % residuals % values(1,1)+1.0_RP, & + CALL FTAssertEqual(expectedValue = residuals_saved(1), & + actualValue = monitors % residuals % values(1,1), & tol = 1.d-11, & msg = "Continuity Residual") - CALL FTAssertEqual(expectedValue = residuals_saved(2)+1.0_RP, & - actualValue = monitors % residuals % values(2,1)+1.0_RP, & + CALL FTAssertEqual(expectedValue = residuals_saved(2), & + actualValue = monitors % residuals % values(2,1), & tol = 1.d-11, & msg = "X-Momentum Residual") - CALL FTAssertEqual(expectedValue = residuals_saved(3)+1.0_RP, & - actualValue = monitors % residuals % values(3,1)+1.0_RP, & + CALL FTAssertEqual(expectedValue = residuals_saved(3), & + actualValue = monitors % residuals % values(3,1), & tol = 1.d-11, & msg = "Y-Momentum Residual") - CALL FTAssertEqual(expectedValue = residuals_saved(4)+1.0_RP, & - actualValue = monitors % residuals % values(4,1)+1.0_RP, & + CALL FTAssertEqual(expectedValue = residuals_saved(4), & + actualValue = monitors % residuals % values(4,1), & tol = 1.d-11, & msg = "Z-Momentum Residual") - CALL FTAssertEqual(expectedValue = residuals_saved(5)+1.0_RP, & - actualValue = monitors % residuals % values(5,1)+1.0_RP, & + CALL FTAssertEqual(expectedValue = residuals_saved(5), & + actualValue = monitors % residuals % values(5,1), & tol = 1.d-11, & msg = "Div Residual") diff --git a/Solver/test/Multiphase/Monopole_pAdaptationRL/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/Monopole_pAdaptationRL/SETUP/ProblemFile.f90 index 0b49ee263d..d82fc68339 100644 --- a/Solver/test/Multiphase/Monopole_pAdaptationRL/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/Monopole_pAdaptationRL/SETUP/ProblemFile.f90 @@ -682,7 +682,7 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & ! Local variables ! --------------- ! - CHARACTER(LEN=29) :: testName = "Re 100 Cylinder acoustics p-Adaptation" + CHARACTER(LEN=29) :: testName = "Monopole P-adaptation RL" REAL(KIND=RP) :: maxError REAL(KIND=RP), ALLOCATABLE :: QExpected(:,:,:,:) INTEGER :: eID @@ -708,11 +708,11 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & ! ----------------------------------------------------------------- ! #if defined(MULTIPHASE) - real(kind=RP), parameter :: residuals(5) = [ 6.7612332904056904E+00_RP, & - 1.3126856239198450E+00_RP, & - 1.3630698469480864E+00_RP, & + real(kind=RP), parameter :: residuals(5) = [ 6.7612400583858490E+00_RP, & + 1.7330535798928657E+00_RP, & + 1.4330968454689006E+00_RP, & 4.5116531519710653E-15_RP, & - 7.2809447388229482E+03_RP] + 5.1783390535161334E+03_RP] CALL initializeSharedAssertionsManager sharedManager => sharedAssertionsManager() diff --git a/Solver/test/Multiphase/Pipe/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/Pipe/SETUP/ProblemFile.f90 index 7b747a01c8..72ad886e80 100644 --- a/Solver/test/Multiphase/Pipe/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/Pipe/SETUP/ProblemFile.f90 @@ -559,11 +559,11 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & TYPE(FTAssertionsManager), POINTER :: sharedManager LOGICAL :: success real(kind=RP), parameter :: residuals_saved(5) = [115.279163488786_RP, & - 218.587148493799_RP, & - 147.034918533523_RP, & - 865.653084759232_RP, & - 72786.6199585908_RP ] - real(kind=RP), parameter :: entropyRate_saved = -49.577609575581214_RP + 221.87180239169732_RP, & + 148.44094157191935_RP, & + 855.40356932873965_RP, & + 72790.732640089103_RP ] + real(kind=RP), parameter :: entropyRate_saved = -49.587092808563099_RP integer :: i CALL initializeSharedAssertionsManager diff --git a/Solver/test/Multiphase/RisingBubble/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/RisingBubble/SETUP/ProblemFile.f90 index 45fff202a4..dbdab6b0d4 100644 --- a/Solver/test/Multiphase/RisingBubble/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/RisingBubble/SETUP/ProblemFile.f90 @@ -488,12 +488,12 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & real(kind=RP), parameter :: area_saved = 2.0280805425949214E-01_RP real(kind=RP), parameter :: xcog_saved = 1.0140441754564830E-01_RP real(kind=RP), parameter :: risevel_saved = 3.8407823048394762E-04_RP - real(kind=RP), parameter :: residuals_saved(5) = [7.2798764748471467E-01_RP, & - 4.1495228250435803E+00_RP, & + real(kind=RP), parameter :: residuals_saved(5) = [7.2798764883600600E-01_RP, & + 4.1495193972979930E+00_RP, & 8.9537579555904704E-14_RP, & - 3.2663420227766222E+00_RP, & - 1.5783130449534988E+02_RP] - real(kind=RP), parameter :: entropyRate_saved = -6.5422251656711682E-03_RP + 3.2663393470370040E+00_RP, & + 1.5783130280842403E+02_RP] + real(kind=RP), parameter :: entropyRate_saved = -6.5422252873259532E-03_RP CALL initializeSharedAssertionsManager sharedManager => sharedAssertionsManager() diff --git a/Solver/test/Multiphase/RisingBubbleVreman/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/RisingBubbleVreman/SETUP/ProblemFile.f90 index 0b85b4fa92..9ffadfa9a6 100644 --- a/Solver/test/Multiphase/RisingBubbleVreman/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/RisingBubbleVreman/SETUP/ProblemFile.f90 @@ -488,12 +488,12 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & real(kind=RP), parameter :: area_saved = 2.0280805425949214E-01_RP real(kind=RP), parameter :: xcog_saved = 1.014044175764861E-01_RP real(kind=RP), parameter :: risevel_saved = 3.840682144321292E-04_RP - real(kind=RP), parameter :: residuals_saved(5) = [7.2798703362831130E-01_RP, & - 4.1415375876201672E+00_RP, & + real(kind=RP), parameter :: residuals_saved(5) = [7.2798703496061990E-01_RP, & + 4.1415341863462345E+00_RP, & 9.7885082257257918E-14_RP, & - 3.2583693013291191E+00_RP, & - 1.5785032659456590E+02_RP] - real(kind=RP), parameter :: entropyRate_saved = -6.5422250365241488E-03_RP + 3.2583666451827610E+00_RP, & + 1.5785032494113227E+02_RP] + real(kind=RP), parameter :: entropyRate_saved = -6.542225158039039E-03_RP CALL initializeSharedAssertionsManager sharedManager => sharedAssertionsManager() diff --git a/Solver/test/Multiphase/Snell/SETUP/ProblemFile.f90 b/Solver/test/Multiphase/Snell/SETUP/ProblemFile.f90 index 55fa3fc689..10db552ce9 100644 --- a/Solver/test/Multiphase/Snell/SETUP/ProblemFile.f90 +++ b/Solver/test/Multiphase/Snell/SETUP/ProblemFile.f90 @@ -546,7 +546,7 @@ SUBROUTINE UserDefinedFinalize(mesh, time, iter, maxResidual & 1.9361250101172587E-03_RP, & 5.2814258414213508E-11_RP, & 2.3047237097072765E-17_RP, & - 7.1908221779661385E-01_RP] + 7.1908221803287240E-01_RP] CALL initializeSharedAssertionsManager