@@ -58,6 +58,7 @@ end subroutine computeBoundaryFluxF
5858 real (kind= RP), protected :: IMEX_S0 = 0.0_RP
5959 real (kind= RP), protected :: IMEX_K0 = 1.0_RP
6060 logical :: use_non_constant_speed_of_sound = .false.
61+ integer :: speed_of_sound_model = 1
6162!
6263! ========
6364 CONTAINS
@@ -192,6 +193,11 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh)
192193 end select
193194
194195 use_non_constant_speed_of_sound = controlVariables % ContainsKey(FLUID1_COMPRESSIBILITY_KEY)
196+ if ( controlVariables % ContainsKey(" speed of sound profile" ) .and. (trim (controlVariables % stringValueForKey(' speed of sound profile' , requestedLength = LINE_LENGTH)) == ' linear quadratic' )) then
197+ speed_of_sound_model = 2
198+ else
199+ speed_of_sound_model = 1
200+ end if
195201
196202 call CHDiscretization % Construct(controlVariables, ELLIPTIC_CH)
197203 call CHDiscretization % Describe
@@ -200,6 +206,12 @@ subroutine Initialize_SpaceAndTimeMethods(controlVariables, mesh)
200206
201207 if (use_non_constant_speed_of_sound) then
202208 write (STD_OUT,' (A)' ) " Implementing artificial compressibility with a non-constant speed of sound in each phase"
209+ if (speed_of_sound_model.eq. 2 ) then
210+ write (STD_OUT,' (A)' ) " Speed of sound profile: linear quadratic along interface"
211+ else
212+ write (STD_OUT,' (A)' ) " Speed of sound profile: linear along interface"
213+ end if
214+
203215 else
204216 write (STD_OUT,' (A)' ) " Implementing artificial compressibility with a constant ACM factor in each phase"
205217 endif
@@ -267,6 +279,24 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
267279!
268280 select case (mode)
269281 case (CTD_IGNORE_MODE,CTD_IMEX_EXPLICIT)
282+
283+ !
284+ ! --------------------------------------------
285+ ! Prolong Cahn-Hilliard and iNS to faces
286+ ! --------------------------------------------
287+ !
288+ call mesh % ProlongSolutionToFaces(NCONS, element_mask= element_mask, Level= locLevel)
289+ !
290+ ! ----------------
291+ ! Update MPI Faces
292+ ! ----------------
293+ !
294+ #ifdef _HAS_MPI_
295+ ! $omp single
296+ call mesh % UpdateMPIFacesSolution(NCONS)
297+ ! $omp end single
298+ #endif
299+
270300! $omp do schedule(runtime) private(eID)
271301 do lID = 1 , MLIter(locLevel,1 )
272302 eID = MLIter_eID(lID)
@@ -278,6 +308,23 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
278308 endif
279309 end do
280310! $omp end do
311+
312+ #ifdef _HAS_MPI_
313+ ! $omp single
314+ call mesh % GatherMPIFacesSolution(NCONS)
315+ ! $omp end single
316+ #endif
317+
318+ ! Not optimized for MLRK and MixedRK but okey
319+ ! $omp do schedule(runtime)
320+ do fID = 1 , size (mesh % faces)
321+ associate(f = > mesh % faces(fID))
322+ f % storage(1 ) % c(1 ,:,:) = f % storage(1 ) % QNS(IMC,:,:)
323+ f % storage(2 ) % c(1 ,:,:) = f % storage(2 ) % QNS(IMC,:,:)
324+ end associate
325+ end do
326+ ! $omp end do
327+
281328 end select
282329!
283330! -------------------------------
@@ -296,12 +343,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
296343 end select
297344! $omp end single
298345!
299- ! --------------------------------------------
300- ! Prolong Cahn-Hilliard concentration to faces
301- ! --------------------------------------------
302- !
303- call mesh % ProlongSolutionToFaces(NCOMP, element_mask= element_mask, Level= locLevel)
304- !
305346! ----------------
306347! Update MPI Faces
307348! ----------------
@@ -453,20 +494,13 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
453494 call SetBoundaryConditionsEqn(NS_BC)
454495! $omp end single
455496!
456- ! -------------------------
457- ! Prolong solution to faces
458- ! -------------------------
459- !
460- call mesh % ProlongSolutionToFaces(NCONS, element_mask= element_mask, Level= locLevel)
461- !
462497! ----------------
463498! Update MPI Faces
464499! ----------------
465500!
466501#ifdef _HAS_MPI_
467502! $omp single
468503 call mesh % UpdateMPIFacesSolution(NCONS)
469- call mesh % GatherMPIFacesSolution(NCONS)
470504! $omp end single
471505#endif
472506!
@@ -486,8 +520,14 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
486520 mesh % elements(eID) % storage % rho = min (max (mesh % elements(eID) % storage % rho, dimensionless % rho_min),dimensionless % rho_max)
487521
488522 if (use_non_constant_speed_of_sound ) then
489- 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
490- mesh % elements(eID) % storage % invMa2 = mesh % elements(eID) % storage % invMa2* mesh % elements(eID) % storage % rho
523+ if (speed_of_sound_model.eq. 1 ) then
524+ ! Linear profile
525+ 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
526+ mesh % elements(eID) % storage % invMa2 = mesh % elements(eID) % storage % invMa2* mesh % elements(eID) % storage % rho
527+ else
528+ ! Linear quadratic profile
529+ mesh % elements(eID) % storage % invMa2 = dimensionless % invMa2(2 ) + (dimensionless % invMa2(1 ) - dimensionless % invMa2(2 )) * min (max (mesh % elements(eID) % storage % Q(IMC,:,:,:), 0.0_RP ), 1.0_RP )
530+ end if
491531 else
492532 mesh % elements(eID) % storage % invMa2 = dimensionless % invMa2(1 )
493533 endif
@@ -504,7 +544,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
504544 if (present (element_mask)) compute_element = face_mask(fID)
505545
506546 if (compute_element) then
507-
508547 mesh % faces(fID) % storage(1 ) % rho = dimensionless % rho(2 ) + (dimensionless % rho(1 )- dimensionless % rho(2 ))* mesh % faces(fID) % storage(1 ) % Q(IMC,:,:)
509548 mesh % faces(fID) % storage(2 ) % rho = dimensionless % rho(2 ) + (dimensionless % rho(1 )- dimensionless % rho(2 ))* mesh % faces(fID) % storage(2 ) % Q(IMC,:,:)
510549
@@ -513,11 +552,19 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
513552
514553
515554 if (use_non_constant_speed_of_sound ) then
516- 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
517- 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
555+ if (speed_of_sound_model.eq. 1 ) then
556+ ! Linear profile
557+ 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
558+ 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
518559
519- mesh % faces(fID) % storage(1 ) % invMa2 = mesh % faces(fID) % storage(1 ) % invMa2* mesh % faces(fID) % storage(1 ) % rho
520- mesh % faces(fID) % storage(2 ) % invMa2 = mesh % faces(fID) % storage(2 ) % invMa2* mesh % faces(fID) % storage(2 ) % rho
560+ mesh % faces(fID) % storage(1 ) % invMa2 = mesh % faces(fID) % storage(1 ) % invMa2* mesh % faces(fID) % storage(1 ) % rho
561+ mesh % faces(fID) % storage(2 ) % invMa2 = mesh % faces(fID) % storage(2 ) % invMa2* mesh % faces(fID) % storage(2 ) % rho
562+
563+ else
564+ ! Linear quadratic profile
565+ mesh % faces(fID) % storage(1 ) % invMa2 = dimensionless % invMa2(2 ) + (dimensionless % invMa2(1 ) - dimensionless % invMa2(2 )) * min (max (mesh % faces(fID) % storage(1 ) % Q(IMC,:,:), 0.0_RP ), 1.0_RP )
566+ mesh % faces(fID) % storage(2 ) % invMa2 = dimensionless % invMa2(2 ) + (dimensionless % invMa2(1 ) - dimensionless % invMa2(2 )) * min (max (mesh % faces(fID) % storage(2 ) % Q(IMC,:,:), 0.0_RP ), 1.0_RP )
567+ end if
521568
522569 else
523570 mesh % faces(fID) % storage(1 ) % invMa2 = dimensionless % invMa2(1 )
@@ -534,16 +581,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
534581!
535582 call ViscousDiscretization % ComputeLocalGradients( NCONS, NCONS, mesh , time , mGradientVariables, Level = locLevel)
536583!
537- ! --------------------
538- ! Update MPI Gradients
539- ! --------------------
540- !
541- #ifdef _HAS_MPI_
542- ! $omp single
543- call mesh % UpdateMPIFacesGradients(NCONS)
544- ! $omp end single
545- #endif
546- !
547584! -------------------------------------
548585! Add the Non-Conservative term to QDot
549586! -------------------------------------
@@ -593,7 +630,6 @@ SUBROUTINE ComputeTimeDerivative( mesh, particles, time, mode, HO_Elements, elem
593630
594631#ifdef _HAS_MPI_
595632! $omp single
596- ! Not sure about the position of this w.r.t the MPI directly above
597633 call mesh % UpdateMPIFacesGradients(NCONS)
598634 call mesh % GatherMPIFacesGradients(NCONS)
599635! $omp end single
0 commit comments