@@ -633,7 +633,7 @@ contains
633633 0.5_wp * mu_R* vflux_R_arr(3 )* vel_R(1 )* (1._wp / dx(j))
634634 end if
635635
636- if (igr_order == 3 ) then
636+ if (igr_order == 3 ) then
637637 E_L = (1._wp / 6._wp )* (2._wp * q_cons_vf(E_idx)%sf(j, k, l) + &
638638 5._wp * q_cons_vf(E_idx)%sf(j + 1 , k, l) - &
639639 1._wp * q_cons_vf(E_idx)%sf(j + 2 , k, l))
@@ -816,124 +816,124 @@ contains
816816 do k = 0 , n
817817 do j = - 1 , m
818818
819- dvel = 0._wp
820- vflux_L_arr = 0._wp
821- vflux_R_arr = 0._wp
819+ dvel = 0._wp
820+ vflux_L_arr = 0._wp
821+ vflux_R_arr = 0._wp
822822
823- #:if MFC_CASE_OPTIMIZATION
824- #:if igr_order == 5
825- !DIR$ unroll 6
826- #:elif igr_irder == 3
827- !DIR$ unroll 4
828- #:endif
823+ #:if MFC_CASE_OPTIMIZATION
824+ #:if igr_order == 5
825+ !DIR$ unroll 6
826+ #:elif igr_irder == 3
827+ !DIR$ unroll 4
829828 #:endif
829+ #:endif
830+ !$acc loop seq
831+ do q = vidxb, vidxe
832+ dvel_small = 0._wp
833+ !x- direction contributions
830834 !$acc loop seq
831- do q = vidxb, vidxe
832- dvel_small = 0._wp
833- !x- direction contributions
835+ do i = - 1 , 1
836+ rho_L = 0._wp
834837 !$acc loop seq
835- do i = - 1 , 1
836- rho_L = 0._wp
837- !$acc loop seq
838- do r = 1 , num_fluids
839- rho_L = rho_L + q_cons_vf(r)%sf(j + i + q, k, l)
840- end do
841- rho_sf_small(i) = rho_L
838+ do r = 1 , num_fluids
839+ rho_L = rho_L + q_cons_vf(r)%sf(j + i + q, k, l)
842840 end do
841+ rho_sf_small(i) = rho_L
842+ end do
843843
844- dvel_small(1 ) = (1 / (2._wp * dx(j)))* ( &
845- q_cons_vf(momxb)%sf(j + 1 + q, k, l)/ rho_sf_small(1 ) - &
846- q_cons_vf(momxb)%sf(j - 1 + q, k, l)/ rho_sf_small(- 1 ))
847- dvel_small(2 ) = (1 / (2._wp * dx(j)))* ( &
848- q_cons_vf(momxb + 1 )%sf(j + 1 + q, k, l)/ rho_sf_small(1 ) - &
849- q_cons_vf(momxb + 1 )%sf(j - 1 + q, k, l)/ rho_sf_small(- 1 ))
850- dvel_small(3 ) = (1 / (2._wp * dx(j)))* ( &
851- q_cons_vf(momxb + 2 )%sf(j + 1 + q, k, l)/ rho_sf_small(1 ) - &
852- q_cons_vf(momxb + 2 )%sf(j - 1 + q, k, l)/ rho_sf_small(- 1 ))
844+ dvel_small(1 ) = (1 / (2._wp * dx(j)))* ( &
845+ q_cons_vf(momxb)%sf(j + 1 + q, k, l)/ rho_sf_small(1 ) - &
846+ q_cons_vf(momxb)%sf(j - 1 + q, k, l)/ rho_sf_small(- 1 ))
847+ dvel_small(2 ) = (1 / (2._wp * dx(j)))* ( &
848+ q_cons_vf(momxb + 1 )%sf(j + 1 + q, k, l)/ rho_sf_small(1 ) - &
849+ q_cons_vf(momxb + 1 )%sf(j - 1 + q, k, l)/ rho_sf_small(- 1 ))
850+ dvel_small(3 ) = (1 / (2._wp * dx(j)))* ( &
851+ q_cons_vf(momxb + 2 )%sf(j + 1 + q, k, l)/ rho_sf_small(1 ) - &
852+ q_cons_vf(momxb + 2 )%sf(j - 1 + q, k, l)/ rho_sf_small(- 1 ))
853853
854- if (q == 0 ) dvel(:, 1 ) = dvel_small
855- if (q > vidxb) then
856- vflux_L_arr(1 ) = vflux_L_arr(1 ) + coeff_L(q)* (dvel_small(2 ))
857- vflux_L_arr(2 ) = vflux_L_arr(2 ) + coeff_L(q)* (dvel_small(3 ))
858- vflux_L_arr(3 ) = vflux_L_arr(3 ) + coeff_L(q)* (4._wp * dvel_small(1 ))/ 3._wp
859- end if
860- if (q < vidxe) then
861- vflux_R_arr(1 ) = vflux_R_arr(1 ) + coeff_R(q)* (dvel_small(2 ))
862- vflux_R_arr(2 ) = vflux_R_arr(2 ) + coeff_R(q)* (dvel_small(3 ))
863- vflux_R_arr(3 ) = vflux_R_arr(3 ) + coeff_R(q)* (4._wp * dvel_small(1 ))/ 3._wp
864- end if
854+ if (q == 0 ) dvel(:, 1 ) = dvel_small
855+ if (q > vidxb) then
856+ vflux_L_arr(1 ) = vflux_L_arr(1 ) + coeff_L(q)* (dvel_small(2 ))
857+ vflux_L_arr(2 ) = vflux_L_arr(2 ) + coeff_L(q)* (dvel_small(3 ))
858+ vflux_L_arr(3 ) = vflux_L_arr(3 ) + coeff_L(q)* (4._wp * dvel_small(1 ))/ 3._wp
859+ end if
860+ if (q < vidxe) then
861+ vflux_R_arr(1 ) = vflux_R_arr(1 ) + coeff_R(q)* (dvel_small(2 ))
862+ vflux_R_arr(2 ) = vflux_R_arr(2 ) + coeff_R(q)* (dvel_small(3 ))
863+ vflux_R_arr(3 ) = vflux_R_arr(3 ) + coeff_R(q)* (4._wp * dvel_small(1 ))/ 3._wp
864+ end if
865865
866- !y- direction contributions
866+ !y- direction contributions
867+ !$acc loop seq
868+ do i = - 1 , 1
869+ rho_L = 0._wp
867870 !$acc loop seq
868- do i = - 1 , 1
869- rho_L = 0._wp
870- !$acc loop seq
871- do r = 1 , num_fluids
872- rho_L = rho_L + q_cons_vf(r)%sf(j + q, k + i, l)
873- end do
874- rho_sf_small(i) = rho_L
871+ do r = 1 , num_fluids
872+ rho_L = rho_L + q_cons_vf(r)%sf(j + q, k + i, l)
875873 end do
874+ rho_sf_small(i) = rho_L
875+ end do
876876
877- dvel_small(1 ) = (1 / (2._wp * dy(k)))* ( &
878- q_cons_vf(momxb)%sf(j + q, k + 1 , l)/ rho_sf_small(1 ) - &
879- q_cons_vf(momxb)%sf(j + q, k - 1 , l)/ rho_sf_small(- 1 ))
880- dvel_small(2 ) = (1 / (2._wp * dy(k)))* ( &
881- q_cons_vf(momxb + 1 )%sf(j + q, k + 1 , l)/ rho_sf_small(1 ) - &
882- q_cons_vf(momxb + 1 )%sf(j + q, k - 1 , l)/ rho_sf_small(- 1 ))
883- if (q == 0 ) dvel_small(3 ) = (1 / (2._wp * dy(k)))* ( &
884- q_cons_vf(momxb + 2 )%sf(j + q, k + 1 , l)/ rho_sf_small(1 ) - &
885- q_cons_vf(momxb + 2 )%sf(j + q, k - 1 , l)/ rho_sf_small(- 1 ))
886- if (q == 0 ) dvel(:, 2 ) = dvel_small
877+ dvel_small(1 ) = (1 / (2._wp * dy(k)))* ( &
878+ q_cons_vf(momxb)%sf(j + q, k + 1 , l)/ rho_sf_small(1 ) - &
879+ q_cons_vf(momxb)%sf(j + q, k - 1 , l)/ rho_sf_small(- 1 ))
880+ dvel_small(2 ) = (1 / (2._wp * dy(k)))* ( &
881+ q_cons_vf(momxb + 1 )%sf(j + q, k + 1 , l)/ rho_sf_small(1 ) - &
882+ q_cons_vf(momxb + 1 )%sf(j + q, k - 1 , l)/ rho_sf_small(- 1 ))
883+ if (q == 0 ) dvel_small(3 ) = (1 / (2._wp * dy(k)))* ( &
884+ q_cons_vf(momxb + 2 )%sf(j + q, k + 1 , l)/ rho_sf_small(1 ) - &
885+ q_cons_vf(momxb + 2 )%sf(j + q, k - 1 , l)/ rho_sf_small(- 1 ))
886+ if (q == 0 ) dvel(:, 2 ) = dvel_small
887887
888- if (q > vidxb) then
889- vflux_L_arr(1 ) = vflux_L_arr(1 ) + coeff_L(q)* (dvel_small(1 ))
890- vflux_L_arr(3 ) = vflux_L_arr(3 ) + coeff_L(q)* (- 2._wp * dvel_small(2 ))/ 3._wp
891- end if
892- if (q < vidxe) then
893- vflux_R_arr(1 ) = vflux_R_arr(1 ) + coeff_R(q)* (dvel_small(1 ))
894- vflux_R_arr(3 ) = vflux_R_arr(3 ) + coeff_R(q)* (- 2._wp * dvel_small(2 ))/ 3._wp
895- end if
888+ if (q > vidxb) then
889+ vflux_L_arr(1 ) = vflux_L_arr(1 ) + coeff_L(q)* (dvel_small(1 ))
890+ vflux_L_arr(3 ) = vflux_L_arr(3 ) + coeff_L(q)* (- 2._wp * dvel_small(2 ))/ 3._wp
891+ end if
892+ if (q < vidxe) then
893+ vflux_R_arr(1 ) = vflux_R_arr(1 ) + coeff_R(q)* (dvel_small(1 ))
894+ vflux_R_arr(3 ) = vflux_R_arr(3 ) + coeff_R(q)* (- 2._wp * dvel_small(2 ))/ 3._wp
895+ end if
896896
897- !z- direction contributions
897+ !z- direction contributions
898+ !$acc loop seq
899+ do i = - 1 , 1
900+ rho_L = 0._wp
898901 !$acc loop seq
899- do i = - 1 , 1
900- rho_L = 0._wp
901- !$acc loop seq
902- do r = 1 , num_fluids
903- rho_L = rho_L + q_cons_vf(r)%sf(j + q, k, l + i)
904- end do
905- rho_sf_small(i) = rho_L
902+ do r = 1 , num_fluids
903+ rho_L = rho_L + q_cons_vf(r)%sf(j + q, k, l + i)
906904 end do
905+ rho_sf_small(i) = rho_L
906+ end do
907907
908- dvel_small(1 ) = (1 / (2._wp * dz(l)))* ( &
909- q_cons_vf(momxb)%sf(j + q, k, l + 1 )/ rho_sf_small(1 ) - &
910- q_cons_vf(momxb)%sf(j + q, k, l - 1 )/ rho_sf_small(- 1 ))
911- if (q == 0 ) dvel_small(2 ) = (1 / (2._wp * dz(l)))* ( &
912- q_cons_vf(momxb + 1 )%sf(j + q, k, l + 1 )/ rho_sf_small(1 ) - &
913- q_cons_vf(momxb + 1 )%sf(j + q, k, l - 1 )/ rho_sf_small(- 1 ))
914- dvel_small(3 ) = (1 / (2._wp * dz(l)))* ( &
915- q_cons_vf(momxb + 2 )%sf(j + q, k, l + 1 )/ rho_sf_small(1 ) - &
916- q_cons_vf(momxb + 2 )%sf(j + q, k, l - 1 )/ rho_sf_small(- 1 ))
917- if (q == 0 ) dvel(:, 3 ) = dvel_small
908+ dvel_small(1 ) = (1 / (2._wp * dz(l)))* ( &
909+ q_cons_vf(momxb)%sf(j + q, k, l + 1 )/ rho_sf_small(1 ) - &
910+ q_cons_vf(momxb)%sf(j + q, k, l - 1 )/ rho_sf_small(- 1 ))
911+ if (q == 0 ) dvel_small(2 ) = (1 / (2._wp * dz(l)))* ( &
912+ q_cons_vf(momxb + 1 )%sf(j + q, k, l + 1 )/ rho_sf_small(1 ) - &
913+ q_cons_vf(momxb + 1 )%sf(j + q, k, l - 1 )/ rho_sf_small(- 1 ))
914+ dvel_small(3 ) = (1 / (2._wp * dz(l)))* ( &
915+ q_cons_vf(momxb + 2 )%sf(j + q, k, l + 1 )/ rho_sf_small(1 ) - &
916+ q_cons_vf(momxb + 2 )%sf(j + q, k, l - 1 )/ rho_sf_small(- 1 ))
917+ if (q == 0 ) dvel(:, 3 ) = dvel_small
918918
919- if (q > vidxb) then
920- vflux_L_arr(2 ) = vflux_L_arr(2 ) + coeff_L(q)* (dvel_small(1 ))
921- vflux_L_arr(3 ) = vflux_L_arr(3 ) + coeff_L(q)* (- 2._wp * dvel_small(3 ))/ 3._wp
922- end if
923- if (q < vidxe) then
924- vflux_R_arr(2 ) = vflux_R_arr(2 ) + coeff_R(q)* (dvel_small(1 ))
925- vflux_R_arr(3 ) = vflux_R_arr(3 ) + coeff_R(q)* (- 2._wp * dvel_small(3 ))/ 3._wp
926- end if
919+ if (q > vidxb) then
920+ vflux_L_arr(2 ) = vflux_L_arr(2 ) + coeff_L(q)* (dvel_small(1 ))
921+ vflux_L_arr(3 ) = vflux_L_arr(3 ) + coeff_L(q)* (- 2._wp * dvel_small(3 ))/ 3._wp
922+ end if
923+ if (q < vidxe) then
924+ vflux_R_arr(2 ) = vflux_R_arr(2 ) + coeff_R(q)* (dvel_small(1 ))
925+ vflux_R_arr(3 ) = vflux_R_arr(3 ) + coeff_R(q)* (- 2._wp * dvel_small(3 ))/ 3._wp
926+ end if
927927
928- if (q == 0 ) then
929- jac_rhs(j, k, l) = alf_igr* (2._wp * (dvel(1 , 2 )* dvel(2 , 1 ) &
930- + dvel(1 , 3 )* dvel(3 , 1 ) &
931- + dvel(2 , 3 )* dvel(3 , 2 )) &
932- + dvel(1 , 1 )** 2._wp + dvel(2 , 2 )** 2._wp &
933- + dvel(3 , 3 )** 2._wp &
934- + (dvel(1 , 1 ) + dvel(2 , 2 ) + dvel(3 , 3 ))** 2._wp )
935- end if
936- end do
928+ if (q == 0 ) then
929+ jac_rhs(j, k, l) = alf_igr* (2._wp * (dvel(1 , 2 )* dvel(2 , 1 ) &
930+ + dvel(1 , 3 )* dvel(3 , 1 ) &
931+ + dvel(2 , 3 )* dvel(3 , 2 )) &
932+ + dvel(1 , 1 )** 2._wp + dvel(2 , 2 )** 2._wp &
933+ + dvel(3 , 3 )** 2._wp &
934+ + (dvel(1 , 1 ) + dvel(2 , 2 ) + dvel(3 , 3 ))** 2._wp )
935+ end if
936+ end do
937937
938938 if (igr_order == 3 ) then
939939 !$acc loop seq
0 commit comments