Skip to content

Commit 541ae5a

Browse files
AbbBalloutAbbBallout
committed
fix exact rieman solver MU
1 parent fd42ad3 commit 541ae5a

17 files changed

Lines changed: 1178 additions & 80 deletions

File tree

‎.github/workflows/CI_sequential 3.yml‎

Lines changed: 25 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -192,9 +192,26 @@ jobs:
192192
source /opt/intel/oneapi/setvars.sh || true
193193
./horses3d.mu Convergence.control
194194
if: ('!cancelled()' && matrix.compiler == 'ifort')
195+
#
196+
# 2) Convergence_cs
197+
# ----------------------------
198+
# This test case is skipped with gfortran as the result is slightly different
195199

200+
- name: Build MultiphaseConvergence_cs
201+
working-directory: ./Solver/test/Multiphase/Convergence_cs/SETUP
202+
run: |
203+
source /opt/intel/oneapi/setvars.sh || true
204+
make MODE=${{matrix.mode}} COMPILER=${{matrix.compiler}} COMM=${{matrix.comm}} ENABLE_THREADS=${{matrix.enable_threads}} WITH_MKL=${{matrix.mkl}}
205+
if: ('!cancelled()' && matrix.compiler == 'ifort')
206+
207+
- name: Run MultiphaseConvergence_cs
208+
working-directory: ./Solver/test/Multiphase/Convergence_cs
209+
run: |
210+
source /opt/intel/oneapi/setvars.sh || true
211+
./horses3d.mu Convergence_cs.control
212+
if: ('!cancelled()' && matrix.compiler == 'ifort')
196213
#
197-
# 2) RisingBubble
214+
# 3) RisingBubble
198215
# ----------------------------
199216
# This test case is skipped with gfortran as the result is slightly different
200217

@@ -213,7 +230,7 @@ jobs:
213230
if: ('!cancelled()' && matrix.compiler == 'ifort')
214231

215232
#
216-
# 3) Pipe
233+
# 4) Pipe
217234
# ----------------------------
218235

219236
- name: Build MultiphasePipe
@@ -231,7 +248,7 @@ jobs:
231248
if: '!cancelled()'
232249

233250
#
234-
# 4) Entropy conserving test
251+
# 5) Entropy conserving test
235252
# ----------------------------
236253

237254
- name: Build MultiphaseEntropyConservingTest
@@ -249,7 +266,7 @@ jobs:
249266
if: '!cancelled()'
250267

251268
#
252-
# 5) RisingBubbleVreman
269+
# 6) RisingBubbleVreman
253270
# ----------------------------
254271
# This test case is skipped with gfortran as the result is slightly different
255272

@@ -268,7 +285,7 @@ jobs:
268285
if: ('!cancelled()')
269286

270287
#
271-
# 6) Mu AL
288+
# 7) Mu AL
272289
# -------------------------------
273290

274291
- name: Build Mu-AL
@@ -286,7 +303,7 @@ jobs:
286303
if: '!cancelled()'
287304

288305
#
289-
# 7) Multiphase monopole with acoustics, sponge, and Reinforcement Learning p-adaptation
306+
# 8) Multiphase monopole with acoustics, sponge, and Reinforcement Learning p-adaptation
290307
# --------------------------------------------
291308
- name: Build Multiphase_Monopole_pAdaptationRL
292309
working-directory: ./Solver/test/Multiphase/Monopole_pAdaptationRL/SETUP
@@ -303,7 +320,7 @@ jobs:
303320
if: '!cancelled()'
304321

305322
#
306-
# 8) Snell
323+
# 9) Snell
307324
# ----------------------------
308325
- name: Build MultiphaseSnell
309326
working-directory: ./Solver/test/Multiphase/Snell/SETUP
@@ -320,7 +337,7 @@ jobs:
320337
if: ('!cancelled()')
321338

322339
#
323-
# 9) Mixed RK
340+
# 10) Mixed RK
324341
# ----------------------------
325342
- name: Build MultiphaseMixedRK
326343
working-directory: ./Solver/test/Multiphase/MixedRK/SETUP

‎Solver/configure‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@ TEST_CASES="./Euler/BoxAroundCircle \
8888
./IncompressibleNS/ActuatorLineInterpolation \
8989
./CahnHilliard/TJokisaari \
9090
./Multiphase/Convergence \
91+
./Multiphase/Convergence_cs \
9192
./Multiphase/EntropyConservingTest \
9293
./Multiphase/RisingBubble \
9394
./Multiphase/RisingBubbleVreman \

‎Solver/src/MultiphaseSolver/SpatialDiscretization.f90‎

Lines changed: 46 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -52,9 +52,12 @@ end subroutine computeBoundaryFluxF
5252

5353
character(len=LINE_LENGTH), parameter :: viscousDiscretizationKey = "viscous discretization"
5454
character(len=LINE_LENGTH), parameter :: CHDiscretizationKey = "cahn-hilliard discretization"
55+
character(len=LINE_LENGTH), parameter :: FLUID1_COMPRESSIBILITY_KEY = "fluid 1 sound speed square (m/s)"
56+
5557

5658
real(kind=RP), protected :: IMEX_S0 = 0.0_RP
5759
real(kind=RP), protected :: IMEX_K0 = 1.0_RP
60+
logical :: use_non_constant_speed_of_sound = .false.
5861
!
5962
! ========
6063
CONTAINS
@@ -188,6 +191,13 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh)
188191

189192
end select
190193

194+
use_non_constant_speed_of_sound = controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY)
195+
if(use_non_constant_speed_of_sound) then
196+
write(STD_OUT,'(A)') " Implementing artificial compressibility with a non-constant speed of sound in each phase"
197+
else
198+
write(STD_OUT,'(A)') " Implementing artificial compressibility with a constant speed of sound in each phase"
199+
endif
200+
191201
call CHDiscretization % Construct(controlVariables, ELLIPTIC_CH)
192202
call CHDiscretization % Describe
193203

@@ -448,7 +458,7 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
448458
#endif
449459
!
450460
! -------------------------------------
451-
! Get the density in faces and elements
461+
! Get the density and invMa2 in faces and elements
452462
! -------------------------------------
453463
!
454464
!$omp do schedule(runtime)
@@ -459,9 +469,17 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
459469
if (compute_element) then
460470

461471
mesh % elements(eID) % storage % rho = dimensionless % rho(2) + (dimensionless % rho(1)-dimensionless % rho(2))*mesh % elements(eID) % storage % Q(IMC,:,:,:)
462-
463472
mesh % elements(eID) % storage % rho = min(max(mesh % elements(eID) % storage % rho, dimensionless % rho_min),dimensionless % rho_max)
473+
474+
if (use_non_constant_speed_of_sound ) then
475+
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
476+
mesh % elements(eID) % storage % invMa2 = mesh % elements(eID) % storage % invMa2*mesh % elements(eID) % storage % rho
477+
else
478+
mesh % elements(eID) % storage % invMa2 = dimensionless % invMa2(1)
479+
endif
480+
464481
endif
482+
465483
end do
466484
!$omp end do nowait
467485

@@ -477,6 +495,20 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
477495

478496
mesh % faces(fID) % storage(1) % rho = min(max(mesh % faces(fID) % storage(1) % rho, dimensionless % rho_min),dimensionless % rho_max)
479497
mesh % faces(fID) % storage(2) % rho = min(max(mesh % faces(fID) % storage(2) % rho, dimensionless % rho_min),dimensionless % rho_max)
498+
499+
500+
if (use_non_constant_speed_of_sound ) then
501+
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
502+
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
503+
504+
mesh % faces(fID) % storage(1) % invMa2 = mesh % faces(fID) % storage(1) % invMa2*mesh % faces(fID) % storage(1) % rho
505+
mesh % faces(fID) % storage(2) % invMa2 = mesh % faces(fID) % storage(2) % invMa2*mesh % faces(fID) % storage(2) % rho
506+
507+
else
508+
mesh % faces(fID) % storage(1) % invMa2 = dimensionless % invMa2(1)
509+
mesh % faces(fID) % storage(2) % invMa2 = dimensionless % invMa2(2)
510+
endif
511+
480512
endif
481513
end do
482514
!$omp end do
@@ -511,7 +543,9 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
511543
associate(e => mesh % elements(eID))
512544
do k = 0, e % Nxyz(3) ; do j = 0, e % Nxyz(2) ; do i = 0, e % Nxyz(1)
513545
sqrtRho = sqrt(e % storage % rho(i,j,k))
514-
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))
546+
547+
invMa2 = e % storage % invMa2(i,j,k)
548+
515549
e % storage % QDot(IMC,i,j,k) = 0.0_RP
516550
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) &
517551
+ e % storage % Q(IMSQRHOV,i,j,k)*e % storage % U_y(IGU,i,j,k) &
@@ -1094,7 +1128,9 @@ SUBROUTINE computeElementInterfaceFlux_MU(f)
10941128
t1 = f % geom % t1(:,i,j), &
10951129
t2 = f % geom % t2(:,i,j), &
10961130
fL = inv_fluxL(:,i,j), &
1097-
fR = inv_fluxR(:,i,j) )
1131+
fR = inv_fluxR(:,i,j), &
1132+
invMa2L= f % storage(1) % invMa2(i,j), &
1133+
invMa2R= f % storage(2) % invMa2(i,j))
10981134

10991135

11001136
!
@@ -1190,7 +1226,9 @@ SUBROUTINE computeMPIFaceFlux_MU(f)
11901226
t1 = f % geom % t1(:,i,j), &
11911227
t2 = f % geom % t2(:,i,j), &
11921228
fL = inv_fluxL(:,i,j), &
1193-
fR = inv_fluxR(:,i,j) )
1229+
fR = inv_fluxR(:,i,j),&
1230+
invMa2L= f % storage(1) % invMa2(i,j), &
1231+
invMa2R= f % storage(2) % invMa2(i,j))
11941232

11951233

11961234
!
@@ -1293,7 +1331,9 @@ SUBROUTINE computeBoundaryFlux_MU(f, time)
12931331
t1 = f % geom % t1(:,i,j), &
12941332
t2 = f % geom % t2(:,i,j), &
12951333
fL = inv_fluxL, &
1296-
fR = inv_fluxR )
1334+
fR = inv_fluxR,&
1335+
invMa2L= f % storage(1) % invMa2(i,j), &
1336+
invMa2R= f % storage(2) % invMa2(i,j))
12971337

12981338
fStar(:,i,j) = (inv_fluxL - visc_flux(:,i,j) ) * f % geom % jacobian(i,j)
12991339
end do

‎Solver/src/libs/mesh/StorageClass.f90‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,9 @@ module StorageClass
106106
real(kind=RP), dimension(:,:,:,:), allocatable :: mu_z ! CHE chemical potential z-gradient
107107
real(kind=RP), dimension(:,:,:,:), allocatable :: v ! CHE flow field velocity
108108
real(kind=RP), dimension(:,:,:,:), allocatable :: G_CH ! CHE auxiliary storage
109+
#endif
110+
#ifdef MULTIPHASE
111+
real(kind=RP), allocatable :: invMa2(:,:,:) !Storage for the density artificial compressibility factor
109112
#endif
110113
contains
111114
procedure :: Assign => ElementStorage_Assign
@@ -192,6 +195,9 @@ module StorageClass
192195
#ifdef ACOUSTIC
193196
real(kind=RP), dimension(:,:,:), allocatable :: Qbase ! Base flow State vector
194197
#endif
198+
#ifdef MULTIPHASE
199+
real(kind=RP), dimension(:,:), allocatable :: invMa2
200+
#endif
195201
!
196202
! Inviscid Jacobians
197203
! ------------------
@@ -858,6 +864,8 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
858864
allocate(self % RKSteps(k) % hatK(1:NCONS,0:Nx, 0:Ny, 0:Nz))
859865
enddo
860866
end if
867+
868+
allocate ( self % invMa2 (0:Nx,0:Ny,0:Nz) )
861869
#endif
862870
!
863871
! -----------------
@@ -905,6 +913,9 @@ elemental subroutine ElementStorage_Construct(self, Nx, Ny, Nz, computeGradients
905913
self % v = 0.0_RP
906914
#endif
907915

916+
#ifdef MULTIPHASE
917+
self % invMa2 = 0.0_RP
918+
#endif
908919
self % first_sensed = huge(1)
909920
self % prev_sensor = 1.0_RP
910921
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)
11211132
safedeallocate(self % G_CH)
11221133
safedeallocate(self % v)
11231134
#endif
1135+
1136+
#ifdef MULTIPHASE
1137+
safedeallocate(self % invMa2)
1138+
#endif
1139+
11241140
safedeallocate(self % PrevQ)
11251141

11261142
end subroutine ElementStorage_Destruct
@@ -1390,6 +1406,10 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
13901406
! -----------------------------------------------------------------------
13911407
interfaceFluxMemorySize = max(interfaceFluxMemorySize, NCOMP*nDIM*product(Nf+1))
13921408
#endif
1409+
1410+
#ifdef MULTIPHASE
1411+
allocate( self % invMa2 (0:Nf(1),0:Nf(2)) )
1412+
#endif
13931413
!
13941414
! Reserve memory for the interface fluxes
13951415
! ---------------------------------------
@@ -1447,6 +1467,11 @@ pure subroutine FaceStorage_Construct(self, NDIM, Nf, Nel, computeGradients, ana
14471467
self % mu_z = 0.0_RP
14481468
self % v = 0.0_RP
14491469
#endif
1470+
1471+
#ifdef MULTIPHASE
1472+
self % invMa2 = 0.0_RP
1473+
#endif
1474+
14501475
self % constructed = .TRUE.
14511476
end subroutine FaceStorage_Construct
14521477
!
@@ -1539,6 +1564,11 @@ elemental subroutine FaceStorage_Destruct(self)
15391564
safedeallocate(self % mu_z)
15401565
safedeallocate(self % v)
15411566
#endif
1567+
1568+
#ifdef MULTIPHASE
1569+
safedeallocate(self % invMa2 )
1570+
#endif
1571+
15421572
safedeallocate(self % genericInterfaceFluxMemory)
15431573

15441574
self % Q => NULL()
@@ -1706,6 +1736,9 @@ elemental subroutine FaceStorage_Assign(to, from)
17061736
to % mu_y = from % mu_y
17071737
to % mu_z = from % mu_z
17081738
to % v = from % v
1739+
#endif
1740+
#ifdef MULTIPHASE
1741+
to % invMa2 = from % invMa2
17091742
#endif
17101743
call to % PointStorage
17111744
end subroutine FaceStorage_Assign

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

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -174,7 +174,7 @@ SUBROUTINE ConstructPhysicsStorage_MU( controlVariables, Lref, timeref, success
174174

175175
if ( controlVariables % ContainsKey(ARTIFICIAL_COMPRESSIBILITY_KEY) .AND. controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY) .OR. &
176176
controlVariables % ContainsKey(ARTIFICIAL_COMPRESSIBILITY_KEY) .AND. controlVariables % ContainsKey(FLUID2_COMPRESSIBILITY_KEY)) then
177-
error stop "Error: Either define a single artificial compressibility or the speed of sound in each fluids. Not both"
177+
error stop "Error: Either define a single artificial compressibility or the speed of sound in each fluid, not both"
178178
endif
179179

180180

@@ -291,7 +291,7 @@ SUBROUTINE ConstructPhysicsStorage_MU( controlVariables, Lref, timeref, success
291291
dimensionless_ % invMa2(1) = dimensionless_ % rho(1) * thermodynamics_ % c02(1) / POW2(refValues_ % V)
292292
endif
293293

294-
if(controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY)) then
294+
if(controlVariables % ContainsKey(FLUID2_COMPRESSIBILITY_KEY)) then
295295
dimensionless_ % Ma2(2) = POW2(refValues_ % V)/(dimensionless_ % rho(2) * thermodynamics_ % c02(2))
296296
dimensionless_ % invMa2(2) = dimensionless_ % rho(2) * thermodynamics_ % c02(2) / POW2(refValues_ % V)
297297
endif

0 commit comments

Comments
 (0)