-
Notifications
You must be signed in to change notification settings - Fork 65
/
m_riemann_solvers.fpp
4833 lines (3979 loc) · 236 KB
/
m_riemann_solvers.fpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!>
!! @file m_riemann_solvers.f90
!! @brief Contains module m_riemann_solvers
!> @brief This module features a database of approximate and exact Riemann
!! problem solvers for the Navier-Stokes system of equations, which
!! is supplemented by appropriate advection equations that are used
!! to capture the material interfaces. The closure of the system is
!! achieved by the stiffened gas equation of state and any required
!! mixture relations. Surface tension effects are accounted for and
!! are modeled by means of a volume force acting across the diffuse
!! material interface region. The implementation details of viscous
!! and capillary effects, into the Riemann solvers, may be found in
!! Perigaud and Saurel (2005). Note that both effects are available
!! only in the volume fraction model. At this time, the approximate
!! and exact Riemann solvers that are listed below are available:
!! 1) Harten-Lax-van Leer (HLL)
!! 2) Harten-Lax-van Leer-Contact (HLLC)
!! 3) Exact
module m_riemann_solvers
! Dependencies =============================================================
use m_derived_types !< Definitions of the derived types
use m_global_parameters !< Definitions of the global parameters
use m_mpi_proxy !< Message passing interface (MPI) module proxy
use m_variables_conversion !< State variables type conversion procedures
use m_bubbles !< To get the bubble wall pressure function
! ==========================================================================
implicit none
private; public :: s_initialize_riemann_solvers_module, &
s_riemann_solver, &
s_hll_riemann_solver, &
s_hllc_riemann_solver, &
s_convert_species_to_mixture_variables_riemann_acc, &
s_finalize_riemann_solvers_module
abstract interface ! =======================================================
!> Abstract interface to the subroutines that are utilized to compute the
!! Riemann problem solution. For additional information please reference:
!! 1) s_hll_riemann_solver
!! 2) s_hllc_riemann_solver
!! 3) s_exact_riemann_solver
!! @param qL_prim_vf The left WENO-reconstructed cell-boundary values of the
!! cell-average primitive variables
!! @param qR_prim_vf The right WENO-reconstructed cell-boundary values of the
!! cell-average primitive variables
!! @param dqL_prim_dx_vf The left WENO-reconstructed cell-boundary values of the
!! first-order x-dir spatial derivatives
!! @param dqL_prim_dy_vf The left WENO-reconstructed cell-boundary values of the
!! first-order y-dir spatial derivatives
!! @param dqL_prim_dz_vf The left WENO-reconstructed cell-boundary values of the
!! first-order z-dir spatial derivatives
!! @param dqR_prim_dx_vf The right WENO-reconstructed cell-boundary values of the
!! first-order x-dir spatial derivatives
!! @param dqR_prim_dy_vf The right WENO-reconstructed cell-boundary values of the
!! first-order y-dir spatial derivatives
!! @param dqR_prim_dz_vf The right WENO-reconstructed cell-boundary values of the
!! first-order z-dir spatial derivatives
!! @param gm_alphaL_vf Left averaged gradient magnitude
!! @param gm_alphaR_vf Right averaged gradient magnitude
!! @param flux_vf Intra-cell fluxes
!! @param flux_src_vf Intra-cell fluxes sources
!! @param flux_gsrc_vf Intra-cell geometric fluxes sources
!! @param norm_dir Dir. splitting direction
!! @param ix Index bounds in the x-dir
!! @param iy Index bounds in the y-dir
!! @param iz Index bounds in the z-dir
!! @param q_prim_vf Cell-averaged primitive variables
subroutine s_abstract_riemann_solver(qL_prim_rsx_vf_flat, qL_prim_rsy_vf_flat, qL_prim_rsz_vf_flat, dqL_prim_dx_vf, &
dqL_prim_dy_vf, &
dqL_prim_dz_vf, &
qL_prim_vf, &
qR_prim_rsx_vf_flat, qR_prim_rsy_vf_flat, qR_prim_rsz_vf_flat, dqR_prim_dx_vf, &
dqR_prim_dy_vf, &
dqR_prim_dz_vf, &
qR_prim_vf, &
q_prim_vf, &
flux_vf, flux_src_vf, &
flux_gsrc_vf, &
norm_dir, ix, iy, iz)
import :: scalar_field, int_bounds_info, sys_size, startx, starty, startz
real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf_flat, qL_prim_rsy_vf_flat, qL_prim_rsz_vf_flat, qR_prim_rsx_vf_flat, qR_prim_rsy_vf_flat, qR_prim_rsz_vf_flat
type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
type(scalar_field), &
allocatable, dimension(:), &
intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
dqL_prim_dy_vf, dqR_prim_dy_vf, &
dqL_prim_dz_vf, dqR_prim_dz_vf
type(scalar_field), &
dimension(sys_size), &
intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
integer, intent(IN) :: norm_dir
type(int_bounds_info), intent(IN) :: ix, iy, iz
end subroutine s_abstract_riemann_solver
!> The abstract interface to the subroutines that are used to calculate
!! the Roe and arithmetic average states. For more information refer to:
!! 1) s_compute_roe_average_state
!! 2) s_compute_arithmetic_average_state
!! @param i First coordinate location index
!! @param j Second coordinate location index
!! @param k Third coordinate location index
subroutine s_compute_abstract_average_state(qL_prim_rs_vf, qR_prim_rs_vf, i, j, k)
import :: scalar_field, int_bounds_info, sys_size
integer, intent(IN) :: i, j, k
type(scalar_field), dimension(sys_size), intent(IN) :: qL_prim_rs_vf, qR_prim_rs_vf
end subroutine s_compute_abstract_average_state
!> The abstract interface to the subroutines that are utilized to compute
!! the wave speeds of the Riemann problem either directly or by the means
!! of pressure-velocity estimates. For more information please refer to:
!! 1) s_compute_direct_wave_speeds
!! 2) s_compute_pressure_velocity_wave_speeds
!! @param i First coordinate location index
!! @param j Second coordinate location index
!! @param k Third coordinate location index
subroutine s_compute_abstract_wave_speeds(i, j, k)
integer, intent(IN) :: i, j, k
end subroutine s_compute_abstract_wave_speeds
!> The abstract interface to the subroutines that are utilized to compute
!! the viscous source fluxes for either Cartesian or cylindrical geometries.
!! For more information please refer to:
!! 1) s_compute_cartesian_viscous_source_flux
!! 2) s_compute_cylindrical_viscous_source_flux
subroutine s_compute_abstract_viscous_source_flux(velL_vf, & ! -------------
dvelL_dx_vf, &
dvelL_dy_vf, &
dvelL_dz_vf, &
velR_vf, &
dvelR_dx_vf, &
dvelR_dy_vf, &
dvelR_dz_vf, &
flux_src_vf, &
norm_dir, &
ix, iy, iz)
import :: scalar_field, int_bounds_info, num_dims, sys_size
type(scalar_field), &
dimension(num_dims), &
intent(IN) :: velL_vf, velR_vf, &
dvelL_dx_vf, dvelR_dx_vf, &
dvelL_dy_vf, dvelR_dy_vf, &
dvelL_dz_vf, dvelR_dz_vf
type(scalar_field), &
dimension(sys_size), &
intent(INOUT) :: flux_src_vf
integer, intent(IN) :: norm_dir
type(int_bounds_info), intent(IN) :: ix, iy, iz
end subroutine s_compute_abstract_viscous_source_flux
end interface ! ============================================================
type(scalar_field), allocatable, dimension(:) :: qL_prim_rs_vf
type(scalar_field), allocatable, dimension(:) :: qR_prim_rs_vf
type(scalar_field), allocatable, dimension(:) :: flux_rs_vf, flux_src_rs_vf
type(scalar_field), allocatable, dimension(:) :: flux_gsrc_rs_vf !<
type(scalar_field), allocatable, dimension(:) :: vel_src_rs_vf
!> The left (L) and right (R) WENO-reconstructed cell-boundary values of the
!! cell-average primitive variables that define the left and right states of
!! the Riemann problem. Variables qK_prim_rs_vf, K = L or R, are obtained by
!! reshaping (RS) qK_prim_vf in a coordinate direction that is normal to the
!! cell-boundaries along which the fluxes are to be determined.
!> @{
type(scalar_field), allocatable, dimension(:) :: qL_prim_rsx_vf
type(scalar_field), allocatable, dimension(:) :: qR_prim_rsx_vf
type(scalar_field), allocatable, dimension(:) :: qL_prim_rsy_vf
type(scalar_field), allocatable, dimension(:) :: qR_prim_rsy_vf
type(scalar_field), allocatable, dimension(:) :: qL_prim_rsz_vf
type(scalar_field), allocatable, dimension(:) :: qR_prim_rsz_vf
!> @}
!> @}
type(scalar_field), allocatable, dimension(:) :: flux_gsrc_rsx_vf !<
type(scalar_field), allocatable, dimension(:) :: flux_gsrc_rsy_vf !<
type(scalar_field), allocatable, dimension(:) :: flux_gsrc_rsz_vf !<
!! The cell-boundary values of the geometrical source flux that are computed
!! through the chosen Riemann problem solver by using the left and right
!! states given in qK_prim_rs_vf. Currently 2D axisymmetric for inviscid only.
! The cell-boundary values of the velocity. vel_src_rs_vf is determined as
! part of Riemann problem solution and is used to evaluate the source flux.
type(scalar_field), allocatable, dimension(:) :: vel_src_rsx_vf
type(scalar_field), allocatable, dimension(:) :: vel_src_rsy_vf
type(scalar_field), allocatable, dimension(:) :: vel_src_rsz_vf
!> @}
!> The cell-boundary values of the fluxes (src - source) that are computed
!! through the chosen Riemann problem solver, and the direct evaluation of
!! source terms, by using the left and right states given in qK_prim_rs_vf,
!! dqK_prim_ds_vf and kappaK_rs_vf, where ds = dx, dy or dz.
!> @{
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsx_vf_flat, flux_src_rsx_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsy_vf_flat, flux_src_rsy_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsz_vf_flat, flux_src_rsz_vf_flat
!> @}
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_gsrc_rsx_vf_flat !<
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_gsrc_rsy_vf_flat !<
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_gsrc_rsz_vf_flat !<
!! The cell-boundary values of the geometrical source flux that are computed
!! through the chosen Riemann problem solver by using the left and right
!! states given in qK_prim_rs_vf. Currently 2D axisymmetric for inviscid only.
! The cell-boundary values of the velocity. vel_src_rs_vf is determined as
! part of Riemann problem solution and is used to evaluate the source flux.
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: vel_src_rsx_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: vel_src_rsy_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: vel_src_rsz_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: mom_sp_rsx_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: mom_sp_rsy_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: mom_sp_rsz_vf_flat
!> @name Left and right, WENO-reconstructed, cell-boundary values of cell-average
!! partial densities, density, velocity, pressure, internal energy, energy, enthalpy, volume
!! fractions, mass fractions, the specific heat ratio and liquid stiffness functions, speed
!! of sound, shear and volume Reynolds numbers and the Weber numbers. These
!! variables are left and right states of the Riemann problem obtained from
!! qK_prim_rs_vf and kappaK_rs_vf.
!> @{
real(kind(0d0)), allocatable, dimension(:) :: alpha_rho_L, alpha_rho_R
real(kind(0d0)) :: rho_L, rho_R
real(kind(0d0)), allocatable, dimension(:) :: vel_L, vel_R
real(kind(0d0)) :: pres_L, pres_R
real(kind(0d0)) :: E_L, E_R
real(kind(0d0)) :: H_L, H_R
real(kind(0d0)), allocatable, dimension(:) :: alpha_L, alpha_R
real(kind(0d0)) :: Y_L, Y_R
real(kind(0d0)) :: gamma_L, gamma_R
real(kind(0d0)) :: pi_inf_L, pi_inf_R
real(kind(0d0)) :: c_L, c_R
real(kind(0d0)), dimension(2) :: Re_L, Re_R
real(kind(0d0)), allocatable, dimension(:) :: tau_e_L, tau_e_R
real(kind(0d0)), allocatable, dimension(:) :: G_L, G_R
!$acc declare create(alpha_rho_L, alpha_rho_R,rho_L, rho_R,vel_L, vel_R,pres_L, pres_R, &
!$acc E_L, E_R, H_L, H_R, alpha_L, alpha_R, Y_L, Y_R, gamma_L, gamma_R,pi_inf_L, pi_inf_R, &
!$acc c_L, c_R,Re_L, Re_R,tau_e_L, tau_e_R, G_L, G_R)
!> @}
!> @name Left and right, WENO-reconstructed, cell-boundary values of cell-average
!! bubble density, radius, radial velocity, pressure, wall pressure, and modified
!! pressure. These variables are left and right states of the Riemann problem obtained from
!! qK_prim_rs_vf and kappaK_rs_vf.
!> @{
real(kind(0d0)) :: nbub_L, nbub_R
real(kind(0d0)), allocatable, dimension(:) :: R0_L, R0_R
real(kind(0d0)), allocatable, dimension(:) :: V0_L, V0_R
real(kind(0d0)), allocatable, dimension(:) :: P0_L, P0_R
real(kind(0d0)), allocatable, dimension(:) :: pbw_L, pbw_R
real(kind(0d0)), allocatable, dimension(:, :) :: moms_L, moms_R
real(kind(0d0)) :: ptilde_L, ptilde_R
!> @}
!$acc declare create(nbub_L, nbub_R, R0_L, R0_R, V0_L, V0_R, P0_L, P0_R, pbw_L, pbw_R, moms_L, moms_R, ptilde_L, ptilde_R )
!> @name Gamma-related constants for use in exact Riemann solver (following Toro (1999) pp.153)
!> @{
real(kind(0d0)) :: G1_L, G1_R
real(kind(0d0)) :: G2_L, G2_R
real(kind(0d0)) :: G3_L, G3_R
real(kind(0d0)) :: G4_L, G4_R
real(kind(0d0)) :: G5_L, G5_R
real(kind(0d0)) :: G6_L, G6_R
real(kind(0d0)) :: G7_L, G7_R
real(kind(0d0)) :: G8_L, G8_R
!> @}
!> @name Star region pressure and velocity
!> @{
real(kind(0d0)) :: pres_S
real(kind(0d0)) :: vel_S
!> @}
!> @name Intercell solution used to calculated intercell flux
!> @{
real(kind(0d0)), allocatable, dimension(:) :: alpha_rho_IC
real(kind(0d0)) :: rho_IC
real(kind(0d0)), allocatable, dimension(:) :: vel_IC
real(kind(0d0)) :: pres_IC
real(kind(0d0)) :: E_IC
real(kind(0d0)), allocatable, dimension(:) :: alpha_IC
real(kind(0d0)), allocatable, dimension(:) :: tau_e_IC
!> @}
!> @name Surface tension pressure contribution
!> @{
real(kind(0d0)) :: dpres_L, dpres_R
!> @}
!$acc declare create(pres_S, vel_S, alpha_IC, alpha_rho_IC, vel_IC, pres_IC, E_IC, rho_IC, tau_e_IC, dpres_L, dpres_R)
!> @name Roe or arithmetic average density, velocity, enthalpy, volume fractions,
!! specific heat ratio function, speed of sound, shear and volume Reynolds
!! numbers, Weber numbers and curvatures, at the cell-boundaries, computed
!! from the left and the right states of the Riemann problem
!> @{
real(kind(0d0)) :: rho_avg
real(kind(0d0)), allocatable, dimension(:) :: vel_avg
real(kind(0d0)) :: H_avg
type(scalar_field), allocatable, dimension(:) :: alpha_avg_rs_vf
real(kind(0d0)) :: gamma_avg
real(kind(0d0)) :: c_avg
type(scalar_field), allocatable, dimension(:) :: Re_avg_rs_vf
type(scalar_field), allocatable, dimension(:) :: Re_avg_rsx_vf
type(scalar_field), allocatable, dimension(:) :: Re_avg_rsy_vf
type(scalar_field), allocatable, dimension(:) :: Re_avg_rsz_vf
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: Re_avg_rsx_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: Re_avg_rsy_vf_flat
real(kind(0d0)), allocatable, dimension(:, :, :, :) :: Re_avg_rsz_vf_flat
!$acc declare create(rho_avg, vel_avg, H_avg, alpha_avg_rs_vf, gamma_avg, c_avg, Re_avg_rs_vf, Re_avg_rsx_vf, Re_avg_rsy_vf, Re_avg_rsz_vf, Re_avg_rsx_vf_flat, Re_avg_rsy_vf_flat, Re_avg_rsz_vf_flat)
!> @}
!> @name Left, right and star (S) region wave speeds
!> @{
real(kind(0d0)) :: s_L, s_R, s_S
!> @}
!> @name Star region variables (HLLC)
!> @{
real(kind(0d0)) :: rho_Star, E_Star, p_Star, p_K_Star
!> @}
!> Minus (M) and plus (P) wave speeds
!> @{
real(kind(0d0)) :: s_M, s_P
!> @}
!> Minus and plus wave speeds functions
!> @{
real(kind(0d0)) :: xi_M, xi_P
!> @}
real(kind(0d0)) :: xi_L, xi_R
!$acc declare create(s_L, s_R, s_S, rho_Star, E_Star, p_Star, p_K_Star, s_M, s_P, xi_M, xi_P, xi_L, xi_R)
procedure(s_abstract_riemann_solver), &
pointer :: s_riemann_solver => null() !<
!! Pointer to the procedure that is utilized to calculate either the HLL,
!! HLLC or exact intercell fluxes, based on the choice of Riemann solver
procedure(s_compute_abstract_average_state), &
pointer :: s_compute_average_state => null() !<
!! Pointer to the subroutine utilized to calculate either the Roe or the
!! arithmetic average state variables, based on the chosen average state
procedure(s_compute_abstract_wave_speeds), &
pointer :: s_compute_wave_speeds => null() !<
!! Pointer to the subroutine that is utilized to compute the wave speeds of
!! the Riemann problem either directly or by the means of pressure-velocity
!! estimates, based on the selected method of estimation of the wave speeds
procedure(s_compute_abstract_viscous_source_flux), &
pointer :: s_compute_viscous_source_flux => null() !<
!! Pointer to the subroutine that is utilized to compute the viscous source
!! flux for either Cartesian or cylindrical geometries.
!> @name Indical bounds in the s1-, s2- and s3-directions
!> @{
type(int_bounds_info) :: is1, is2, is3
type(int_bounds_info) :: isx, isy, isz
!> @}
!$acc declare create(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, &
!$acc is1, is2, is3, isx, isy, isz, vel_src_rsx_vf, vel_src_rsy_vf, vel_src_rsz_vf, &
!$acc flux_gsrc_rsx_vf, flux_gsrc_rsy_vf, flux_gsrc_rsz_vf)
!$acc declare create(&
!$acc flux_rsx_vf_flat, flux_src_rsx_vf_flat, flux_rsy_vf_flat, flux_src_rsy_vf_flat, flux_rsz_vf_flat, flux_src_rsz_vf_flat, vel_src_rsx_vf_flat, vel_src_rsy_vf_flat, vel_src_rsz_vf_flat, &
!$acc flux_gsrc_rsx_vf_flat, flux_gsrc_rsy_vf_flat, flux_gsrc_rsz_vf_flat, mom_sp_rsx_vf_flat, mom_sp_rsy_vf_flat, mom_sp_rsz_vf_flat)
integer :: momxb, momxe
integer :: contxb, contxe
integer :: advxb, advxe
integer :: bubxb, bubxe
integer :: intxb, intxe
integer :: strxb, strxe
!$acc declare create(momxb, momxe, contxb, contxe, advxb, advxe, bubxb, bubxe, intxb, intxe, strxb, strxe)
real(kind(0d0)), allocatable, dimension(:) :: gammas, pi_infs, Gs
!$acc declare create(gammas, pi_infs, Gs)
integer, allocatable, dimension(:) :: rs, vs, ps, ms
!$acc declare create(rs, vs, ps, ms)
real(kind(0d0)), allocatable, dimension(:, :) :: Res
!$acc declare create(Res)
contains
subroutine s_hll_riemann_solver(qL_prim_rsx_vf_flat, qL_prim_rsy_vf_flat, qL_prim_rsz_vf_flat, dqL_prim_dx_vf, & ! -------
dqL_prim_dy_vf, &
dqL_prim_dz_vf, &
qL_prim_vf, &
qR_prim_rsx_vf_flat, qR_prim_rsy_vf_flat, qR_prim_rsz_vf_flat, dqR_prim_dx_vf, &
dqR_prim_dy_vf, &
dqR_prim_dz_vf, &
qR_prim_vf, &
q_prim_vf, &
flux_vf, flux_src_vf, &
flux_gsrc_vf, &
norm_dir, ix, iy, iz)
real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf_flat, qL_prim_rsy_vf_flat, qL_prim_rsz_vf_flat, qR_prim_rsx_vf_flat, qR_prim_rsy_vf_flat, qR_prim_rsz_vf_flat
type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
type(scalar_field), &
allocatable, dimension(:), &
intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
dqL_prim_dy_vf, dqR_prim_dy_vf, &
dqL_prim_dz_vf, dqR_prim_dz_vf
! Intercell fluxes
type(scalar_field), &
dimension(sys_size), &
intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
integer, intent(IN) :: norm_dir
type(int_bounds_info), intent(IN) :: ix, iy, iz
real(kind(0d0)), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R
real(kind(0d0)) :: rho_L, rho_R
real(kind(0d0)), dimension(num_dims) :: vel_L, vel_R
real(kind(0d0)) :: pres_L, pres_R
real(kind(0d0)) :: E_L, E_R
real(kind(0d0)) :: H_L, H_R
real(kind(0d0)), dimension(num_fluids) :: alpha_L, alpha_R
real(kind(0d0)) :: Y_L, Y_R
real(kind(0d0)) :: gamma_L, gamma_R
real(kind(0d0)) :: pi_inf_L, pi_inf_R
real(kind(0d0)) :: c_L, c_R
real(kind(0d0)), dimension(6) :: tau_e_L, tau_e_R
real(kind(0d0)) :: G_L, G_R
real(kind(0d0)), dimension(2) :: Re_L, Re_R
real(kind(0d0)) :: rho_avg
real(kind(0d0)), dimension(num_dims) :: vel_avg
real(kind(0d0)) :: H_avg
real(kind(0d0)) :: gamma_avg
real(kind(0d0)) :: c_avg
real(kind(0d0)) :: s_L, s_R, s_M, s_P, s_S
real(kind(0d0)) :: xi_L, xi_R !< Left and right wave speeds functions
real(kind(0d0)) :: xi_M, xi_P
real(kind(0d0)) :: nbub_L, nbub_R
real(kind(0d0)) :: ptilde_L, ptilde_R
real(kind(0d0)) :: vel_L_rms, vel_R_rms, vel_avg_rms
real(kind(0d0)) :: blkmod1, blkmod2
real(kind(0d0)) :: rho_Star, E_Star, p_Star, p_K_Star
real(kind(0d0)) :: Ms_L, Ms_R, pres_SL, pres_SR
real(kind(0d0)) :: alpha_L_sum, alpha_R_sum
integer :: i, j, k, l, q !< Generic loop iterators
! Populating the buffers of the left and right Riemann problem
! states variables, based on the choice of boundary conditions
call s_populate_riemann_states_variables_buffers( &
qL_prim_rsx_vf_flat, qL_prim_rsy_vf_flat, qL_prim_rsz_vf_flat, dqL_prim_dx_vf, &
dqL_prim_dy_vf, &
dqL_prim_dz_vf, &
qL_prim_vf, &
qR_prim_rsx_vf_flat, qR_prim_rsy_vf_flat, qR_prim_rsz_vf_flat, dqR_prim_dx_vf, &
dqR_prim_dy_vf, &
dqR_prim_dz_vf, &
qR_prim_vf, &
norm_dir, ix, iy, iz)
! Reshaping inputted data based on dimensional splitting direction
call s_initialize_riemann_solver( &
q_prim_vf, &
flux_vf, flux_src_vf, &
flux_gsrc_vf, &
norm_dir, ix, iy, iz)
#:for NORM_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
if (norm_dir == ${NORM_DIR}$) then
!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, vel_avg, tau_e_L, tau_e_R, G_L, G_R, Re_L, Re_R, &
!$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S)
do l = is3%beg, is3%end
do k = is2%beg, is2%end
do j = is1%beg, is1%end
!$acc loop seq
do i = 1, contxe
alpha_rho_L(i) = qL_prim_rs${XYZ}$_vf_flat(j, k, l, i)
alpha_rho_R(i) = qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, i)
end do
!$acc loop seq
do i = 1, num_dims
vel_L(i) = qL_prim_rs${XYZ}$_vf_flat(j, k, l, contxe + i)
vel_R(i) = qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, contxe + i)
end do
vel_L_rms = 0d0; vel_R_rms = 0d0
!$acc loop seq
do i = 1, num_dims
vel_L_rms = vel_L_rms + vel_L(i)**2d0
vel_R_rms = vel_R_rms + vel_R(i)**2d0
end do
!$acc loop seq
do i = 1, num_fluids
alpha_L(i) = qL_prim_rs${XYZ}$_vf_flat(j, k, l, E_idx + i)
alpha_R(i) = qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, E_idx + i)
end do
pres_L = qL_prim_rs${XYZ}$_vf_flat(j, k, l, E_idx)
pres_R = qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, E_idx)
rho_L = 0d0
gamma_L = 0d0
pi_inf_L = 0d0
rho_R = 0d0
gamma_R = 0d0
pi_inf_R = 0d0
alpha_L_sum = 0d0
alpha_R_sum = 0d0
if (mpp_lim) then
!$acc loop seq
do i = 1, num_fluids
alpha_rho_L(i) = max(0d0, alpha_rho_L(i))
alpha_L(i) = min(max(0d0, alpha_L(i)), 1d0)
alpha_L_sum = alpha_L_sum + alpha_L(i)
end do
alpha_L = alpha_L/max(alpha_L_sum, sgm_eps)
!$acc loop seq
do i = 1, num_fluids
alpha_rho_R(i) = max(0d0, alpha_rho_R(i))
alpha_R(i) = min(max(0d0, alpha_R(i)), 1d0)
alpha_R_sum = alpha_R_sum + alpha_R(i)
end do
alpha_R = alpha_R/max(alpha_R_sum, sgm_eps)
end if
!$acc loop seq
do i = 1, num_fluids
rho_L = rho_L + alpha_rho_L(i)
gamma_L = gamma_L + alpha_L(i)*gammas(i)
pi_inf_L = pi_inf_L + alpha_L(i)*pi_infs(i)
rho_R = rho_R + alpha_rho_R(i)
gamma_R = gamma_R + alpha_R(i)*gammas(i)
pi_inf_R = pi_inf_R + alpha_R(i)*pi_infs(i)
end do
if (any(Re_size > 0)) then
!$acc loop seq
do i = 1, 2
Re_L(i) = dflt_real
if (Re_size(i) > 0) Re_L(i) = 0d0
!$acc loop seq
do q = 1, Re_size(i)
Re_L(i) = alpha_L(Re_idx(i, q))/Res(i, q) &
+ Re_L(i)
end do
Re_L(i) = 1d0/max(Re_L(i), sgm_eps)
end do
!$acc loop seq
do i = 1, 2
Re_R(i) = dflt_real
if (Re_size(i) > 0) Re_R(i) = 0d0
!$acc loop seq
do q = 1, Re_size(i)
Re_R(i) = alpha_R(Re_idx(i, q))/Res(i, q) &
+ Re_R(i)
end do
Re_R(i) = 1d0/max(Re_R(i), sgm_eps)
end do
end if
E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms
E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms
H_L = (E_L + pres_L)/rho_L
H_R = (E_R + pres_R)/rho_R
if (hypoelasticity) then
!$acc loop seq
do i = 1, strxe - strxb + 1
tau_e_L(i) = qL_prim_rs${XYZ}$_vf_flat(j, k, l, strxb - 1 + i)
tau_e_R(i) = qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, strxb - 1 + i)
end do
G_L = 0d0
G_R = 0d0
!$acc loop seq
do i = 1, num_fluids
G_L = G_L + alpha_L(i)*Gs(i)
G_R = G_R + alpha_R(i)*Gs(i)
end do
do i = 1, strxe - strxb + 1
! Elastic contribution to energy if G large enough
!TODO take out if statement if stable without
if ((G_L > 1000) .and. (G_R > 1000)) then
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4d0*G_L)
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4d0*G_R)
! Additional terms in 2D and 3D
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4d0*G_L)
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4d0*G_R)
end if
end if
end do
end if
if (avg_state == 2) then
rho_avg = 5d-1*(rho_L + rho_R)
!$acc loop seq
do i = 1, num_dims
vel_avg(i) = 5d-1*(vel_L(i) + vel_R(i))
end do
H_avg = 5d-1*(H_L + H_R)
gamma_avg = 5d-1*(gamma_L + gamma_R)
elseif (avg_state == 1) then
rho_avg = sqrt(rho_L*rho_R)
!$acc loop seq
do i = 1, num_dims
vel_avg(i) = (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))/ &
(sqrt(rho_L) + sqrt(rho_R))
end do
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
(sqrt(rho_L) + sqrt(rho_R))
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
(sqrt(rho_L) + sqrt(rho_R))
end if
vel_avg_rms = 0d0
!$acc loop seq
do i = 1, num_dims
vel_avg_rms = vel_avg_rms + vel_avg(i)**2d0
end do
if (mixture_err) then
if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then
c_avg = sgm_eps
else
c_avg = sqrt((H_avg - 5d-1*vel_avg_rms)/gamma_avg)
end if
else
c_avg = sqrt((H_avg - 5d-1*vel_avg_rms)/gamma_avg)
end if
if (alt_soundspeed) then
blkmod1 = ((gammas(1) + 1d0)*pres_L + &
pi_infs(1))/gammas(1)
blkmod2 = ((gammas(2) + 1d0)*pres_L + &
pi_infs(2))/gammas(2)
c_L = 1d0/(rho_L*(alpha_L(1)/blkmod1 + alpha_L(2)/blkmod2))
blkmod1 = ((gammas(1) + 1d0)*pres_R + &
pi_infs(1))/gammas(1)
blkmod2 = ((gammas(2) + 1d0)*pres_R + &
pi_infs(2))/gammas(2)
c_R = 1d0/(rho_R*(alpha_R(1)/blkmod1 + alpha_R(2)/blkmod2))
elseif (model_eqns == 3) then
c_L = 0d0
c_R = 0d0
!$acc loop seq
do i = 1, num_fluids
c_L = c_L + qL_prim_rs${XYZ}$_vf_flat(j, k, l, i + advxb - 1)*(1d0/gammas(i) + 1d0)* &
(qL_prim_rs${XYZ}$_vf_flat(j, k, l, E_idx) + pi_infs(i)/(gammas(i) + 1d0))
c_R = c_R + qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, i + advxb - 1)*(1d0/gammas(i) + 1d0)* &
(qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, E_idx) + pi_infs(i)/(gammas(i) + 1d0))
end do
c_L = c_L/rho_L
c_R = c_R/rho_R
elseif ((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles)) then
! Sound speed for bubble mmixture to order O(\alpha)
if (mpp_lim .and. (num_fluids > 1)) then
c_L = (1d0/gamma_L + 1d0)* &
(pres_L + pi_inf_L)/rho_L
c_R = (1d0/gamma_R + 1d0)* &
(pres_R + pi_inf_R)/rho_R
else
c_L = &
(1d0/gamma_L + 1d0)* &
(pres_L + pi_inf_L)/ &
(rho_L*(1d0 - alpha_L(num_fluids)))
c_R = &
(1d0/gamma_R + 1d0)* &
(pres_R + pi_inf_R)/ &
(rho_R*(1d0 - alpha_R(num_fluids)))
end if
else
c_L = ((H_L - 5d-1*vel_L_rms)/gamma_L)
c_R = ((H_R - 5d-1*vel_R_rms)/gamma_R)
end if
if (mixture_err .and. c_L < 0d0) then
c_L = 100.d0*sgm_eps
else
c_L = sqrt(c_L)
end if
if (mixture_err .and. c_R < 0d0) then
c_R = 100.d0*sgm_eps
else
c_R = sqrt(c_R)
end if
if (any(Re_size > 0)) then
!$acc loop seq
do i = 1, 2
Re_avg_rs${XYZ}$_vf_flat(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i))
end do
end if
if (wave_speeds == 1) then
if (hypoelasticity) then
s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + &
(((4d0*G_L)/3d0) + &
tau_e_L(dir_idx_tau(1)))/rho_L) &
, vel_R(dir_idx(1)) - sqrt(c_R*c_R + &
(((4d0*G_R)/3d0) + &
tau_e_R(dir_idx_tau(1)))/rho_R))
s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + &
(((4d0*G_R)/3d0) + &
tau_e_R(dir_idx_tau(1)))/rho_R) &
, vel_L(dir_idx(1)) + sqrt(c_L*c_L + &
(((4d0*G_L)/3d0) + &
tau_e_L(dir_idx_tau(1)))/rho_L))
else
s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L)
end if
s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* &
(s_L - vel_L(dir_idx(1))) - &
rho_R*vel_R(dir_idx(1))* &
(s_R - vel_R(dir_idx(1)))) &
/(rho_L*(s_L - vel_L(dir_idx(1))) - &
rho_R*(s_R - vel_R(dir_idx(1))))
elseif (wave_speeds == 2) then
pres_SL = 5d-1*(pres_L + pres_R + rho_avg*c_avg* &
(vel_L(dir_idx(1)) - &
vel_R(dir_idx(1))))
pres_SR = pres_SL
Ms_L = max(1d0, sqrt(1d0 + ((5d-1 + gamma_L)/(1d0 + gamma_L))* &
(pres_SL/pres_L - 1d0)*pres_L/ &
((pres_L + pi_inf_L/(1d0 + gamma_L)))))
Ms_R = max(1d0, sqrt(1d0 + ((5d-1 + gamma_R)/(1d0 + gamma_R))* &
(pres_SR/pres_R - 1d0)*pres_R/ &
((pres_R + pi_inf_R/(1d0 + gamma_R)))))
s_L = vel_L(dir_idx(1)) - c_L*Ms_L
s_R = vel_R(dir_idx(1)) + c_R*Ms_R
s_S = 5d-1*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + &
(pres_L - pres_R)/ &
(rho_avg*c_avg))
end if
s_M = min(0d0, s_L); s_P = max(0d0, s_R)
xi_M = (5d-1 + sign(5d-1, s_L)) &
+ (5d-1 - sign(5d-1, s_L)) &
*(5d-1 + sign(5d-1, s_R))
xi_P = (5d-1 - sign(5d-1, s_R)) &
+ (5d-1 - sign(5d-1, s_L)) &
*(5d-1 + sign(5d-1, s_R))
! Mass
!$acc loop seq
do i = 1, contxe
flux_rs${XYZ}$_vf_flat(j, k, l, i) = &
(s_M*alpha_rho_R(i)*vel_R(dir_idx(1)) &
- s_P*alpha_rho_L(i)*vel_L(dir_idx(1)) &
+ s_M*s_P*(alpha_rho_L(i) &
- alpha_rho_R(i))) &
/(s_M - s_P)
end do
! Momentum
if (bubbles) then
!$acc loop seq
do i = 1, num_dims
flux_rs${XYZ}$_vf_flat(j, k, l, contxe + dir_idx(i)) = &
(s_M*(rho_R*vel_R(dir_idx(1)) &
*vel_R(dir_idx(i)) &
+ dir_flg(dir_idx(i))*(pres_R - ptilde_R)) &
- s_P*(rho_L*vel_L(dir_idx(1)) &
*vel_L(dir_idx(i)) &
+ dir_flg(dir_idx(i))*(pres_L - ptilde_L)) &
+ s_M*s_P*(rho_L*vel_L(dir_idx(i)) &
- rho_R*vel_R(dir_idx(i)))) &
/(s_M - s_P)
end do
else if (hypoelasticity) then
!$acc loop seq
do i = 1, num_dims
flux_rs${XYZ}$_vf_flat(j, k, l, contxe + dir_idx(i)) = &
(s_M*(rho_R*vel_R(dir_idx(1)) &
*vel_R(dir_idx(i)) &
+ dir_flg(dir_idx(i))*pres_R &
- tau_e_R(dir_idx_tau(i))) &
- s_P*(rho_L*vel_L(dir_idx(1)) &
*vel_L(dir_idx(i)) &
+ dir_flg(dir_idx(i))*pres_L &
- tau_e_L(dir_idx_tau(i))) &
+ s_M*s_P*(rho_L*vel_L(dir_idx(i)) &
- rho_R*vel_R(dir_idx(i)))) &
/(s_M - s_P)
end do
else
!$acc loop seq
do i = 1, num_dims
flux_rs${XYZ}$_vf_flat(j, k, l, contxe + dir_idx(i)) = &
(s_M*(rho_R*vel_R(dir_idx(1)) &
*vel_R(dir_idx(i)) &
+ dir_flg(dir_idx(i))*pres_R) &
- s_P*(rho_L*vel_L(dir_idx(1)) &
*vel_L(dir_idx(i)) &
+ dir_flg(dir_idx(i))*pres_L) &
+ s_M*s_P*(rho_L*vel_L(dir_idx(i)) &
- rho_R*vel_R(dir_idx(i)))) &
/(s_M - s_P)
end do
end if
! Energy
if (bubbles) then
flux_rs${XYZ}$_vf_flat(j, k, l, E_idx) = &
(s_M*vel_R(dir_idx(1))*(E_R + pres_R - ptilde_R) &
- s_P*vel_L(dir_idx(1))*(E_L + pres_L - ptilde_L) &
+ s_M*s_P*(E_L - E_R)) &
/(s_M - s_P)
else if (hypoelasticity) then
!TODO: simplify this so it's not split into 3
if (num_dims == 1) then
flux_rs${XYZ}$_vf_flat(j, k, l, E_idx) = &
(s_M*(vel_R(dir_idx(1))*(E_R + pres_R) &
- (tau_e_R(dir_idx_tau(1))*vel_R(dir_idx(1)))) &
- s_P*(vel_L(dir_idx(1))*(E_L + pres_L) &
- (tau_e_L(dir_idx_tau(1))*vel_L(dir_idx(1)))) &
+ s_M*s_P*(E_L - E_R)) &
/(s_M - s_P)
else if (num_dims == 2) then
flux_rs${XYZ}$_vf_flat(j, k, l, E_idx) = &
(s_M*(vel_R(dir_idx(1))*(E_R + pres_R) &
- (tau_e_R(dir_idx_tau(1))*vel_R(dir_idx(1))) &
- (tau_e_R(dir_idx_tau(2))*vel_R(dir_idx(2)))) &
- s_P*(vel_L(dir_idx(1))*(E_L + pres_L) &
- (tau_e_L(dir_idx_tau(1))*vel_L(dir_idx(1))) &
- (tau_e_L(dir_idx_tau(2))*vel_L(dir_idx(2)))) &
+ s_M*s_P*(E_L - E_R)) &
/(s_M - s_P)
else if (num_dims == 3) then
flux_rs${XYZ}$_vf_flat(j, k, l, E_idx) = &
(s_M*(vel_R(dir_idx(1))*(E_R + pres_R) &
- (tau_e_R(dir_idx_tau(1))*vel_R(dir_idx(1))) &
- (tau_e_R(dir_idx_tau(2))*vel_R(dir_idx(2))) &
- (tau_e_R(dir_idx_tau(3))*vel_R(dir_idx(3)))) &
- s_P*(vel_L(dir_idx(1))*(E_L + pres_L) &
- (tau_e_L(dir_idx_tau(1))*vel_L(dir_idx(1))) &
- (tau_e_L(dir_idx_tau(2))*vel_L(dir_idx(2))) &
- (tau_e_L(dir_idx_tau(3))*vel_L(dir_idx(3)))) &
+ s_M*s_P*(E_L - E_R)) &
/(s_M - s_P)
end if
else
flux_rs${XYZ}$_vf_flat(j, k, l, E_idx) = &
(s_M*vel_R(dir_idx(1))*(E_R + pres_R) &
- s_P*vel_L(dir_idx(1))*(E_L + pres_L) &
+ s_M*s_P*(E_L - E_R)) &
/(s_M - s_P)
end if
! Elastic Stresses
if (hypoelasticity) then
do i = 1, strxe - strxb + 1 !TODO: this indexing may be slow
flux_rs${XYZ}$_vf_flat(j, k, l, strxb - 1 + i) = &
(s_M*(rho_R*vel_R(dir_idx(1)) &
*tau_e_R(i)) &
- s_P*(rho_L*vel_L(dir_idx(1)) &
*tau_e_L(i)) &
+ s_M*s_P*(rho_L*tau_e_L(i) &
- rho_R*tau_e_R(i))) &
/(s_M - s_P)
end do
end if
! Advection
!$acc loop seq
do i = advxb, advxe
flux_rs${XYZ}$_vf_flat(j, k, l, i) = &
(qL_prim_rs${XYZ}$_vf_flat(j, k, l, i) &
- qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, i)) &
*s_M*s_P/(s_M - s_P)
flux_src_rs${XYZ}$_vf_flat(j, k, l, i) = &
(s_M*qR_prim_rs${XYZ}$_vf_flat(j + 1, k, l, i) &
- s_P*qL_prim_rs${XYZ}$_vf_flat(j, k, l, i)) &
/(s_M - s_P)
end do
! Div(U)?
!$acc loop seq
do i = 1, num_dims
vel_src_rs${XYZ}$_vf_flat(j, k, l, dir_idx(i)) = &
(xi_M*(rho_L*vel_L(dir_idx(i))* &
(s_L - vel_L(dir_idx(1))) - &
pres_L*dir_flg(dir_idx(i))) - &
xi_P*(rho_R*vel_R(dir_idx(i))* &
(s_R - vel_R(dir_idx(1))) - &
pres_R*dir_flg(dir_idx(i)))) &
/(xi_M*rho_L*(s_L - vel_L(dir_idx(1))) - &
xi_P*rho_R*(s_R - vel_R(dir_idx(1))))
end do
if (bubbles) then
! From HLLC: Kills mass transport @ bubble gas density
if (num_fluids > 1) then
flux_rs${XYZ}$_vf_flat(j, k, l, contxe) = 0d0
end if
end if
end do
end do
end do
end if
#:endfor
if (any(Re_size > 0)) then
if (weno_Re_flux) then
call s_compute_viscous_source_flux( &
qL_prim_vf(momxb:momxe), &
dqL_prim_dx_vf(momxb:momxe), &
dqL_prim_dy_vf(momxb:momxe), &