Add viscous and surface_tension logicals

src/common/m_mpi_common.fpp

@@ -245,7 +245,7 @@ contains
                         MPI_DOUBLE_PRECISION, MPI_MAX, 0, &
                         MPI_COMM_WORLD, ierr)
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             call MPI_REDUCE(vcfl_max_loc, vcfl_max_glb, 1, &
                             MPI_DOUBLE_PRECISION, MPI_MAX, 0, &
                             MPI_COMM_WORLD, ierr)
@@ -258,7 +258,7 @@ contains
 
         icfl_max_glb = icfl_max_loc
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             vcfl_max_glb = vcfl_max_loc
             Rc_min_glb = Rc_min_loc
         end if

src/common/m_variables_conversion.fpp

@@ -340,7 +340,7 @@ contains
 
 #ifdef MFC_SIMULATION
         ! Computing the shear and bulk Reynolds numbers from species analogs
-        if (any(Re_size > 0)) then
+        if (viscous) then
             if (num_fluids == 1) then ! need to consider case with num_fluids >= 2
                 do i = 1, 2
 
@@ -530,7 +530,7 @@ contains
             G_K = max(0d0, G_K)
         end if
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
 
             do i = 1, 2
                 Re_K(i) = dflt_real
@@ -596,7 +596,7 @@ contains
             qv_K = qvs(1)
         end if
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             if (num_fluids == 1) then ! need to consider case with num_fluids >= 2
 
                 do i = 1, 2
@@ -661,7 +661,7 @@ contains
 
 #ifdef MFC_SIMULATION
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size)))
             do i = 1, 2
                 do j = 1, Re_size(i)

src/post_process/m_checker.fpp

@@ -115,7 +115,7 @@ contains
 
     !> Checks constraints on surface tension parameters (cf_wrt and sigma)
     subroutine s_check_inputs_surface_tension
-        @:PROHIBIT(cf_wrt .and. f_is_default(sigma), &
+        @:PROHIBIT(cf_wrt .and. .not. surface_tension, &
             "cf_wrt can only be enabled if the surface coefficient is set")
     end subroutine s_check_inputs_surface_tension
 

src/post_process/m_global_parameters.fpp

@@ -275,6 +275,7 @@ module m_global_parameters
     !> @name surface tension coefficient
     !> @{
     real(kind(0d0)) :: sigma
+    logical :: surface_tension
     !> #}
 
     !> @name Index variables used for m_variables_conversion
@@ -397,6 +398,7 @@ contains
         poly_sigma = dflt_real
         sigR = dflt_real
         sigma = dflt_real
+        surface_tension = .false.
         adv_n = .false.
 
     end subroutine s_assign_default_values_to_user_inputs
@@ -534,7 +536,7 @@ contains
                 sys_size = stress_idx%end
             end if
 
-            if (.not. f_is_default(sigma)) then
+            if (surface_tension) then
                 c_idx = sys_size + 1
                 sys_size = c_idx
             end if
@@ -559,7 +561,7 @@ contains
             sys_size = internalEnergies_idx%end
             alf_idx = 1 ! dummy, cannot actually have a void fraction
 
-            if (.not. f_is_default(sigma)) then
+            if (surface_tension) then
                 c_idx = sys_size + 1
                 sys_size = c_idx
             end if

src/post_process/m_start_up.f90

@@ -83,7 +83,7 @@ subroutine s_read_input_file
             polydisperse, poly_sigma, file_per_process, relax, &
             relax_model, cf_wrt, sigma, adv_n, ib, &
             cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
-            cfl_target
+            cfl_target, surface_tension
 
         ! Inquiring the status of the post_process.inp file
         file_loc = 'post_process.inp'

src/pre_process/m_assign_variables.fpp

@@ -663,7 +663,7 @@ contains
             end do
         end if
 
-        if (.not. f_is_default(sigma)) then
+        if (surface_tension) then
             q_prim_vf(c_idx)%sf(j, k, l) = eta*patch_icpp(patch_id)%cf_val + &
                                            (1d0 - eta)*patch_icpp(smooth_patch_id)%cf_val
         end if

src/pre_process/m_global_parameters.fpp

@@ -228,6 +228,7 @@ module m_global_parameters
     !> @name Surface Tension Modeling
     !> @{
     real(kind(0d0)) :: sigma
+    logical :: surface_tension
     !> @}
 
     !> @name Index variables used for m_variables_conversion
@@ -400,6 +401,7 @@ contains
         Re_inv = dflt_real
         Web = dflt_real
         poly_sigma = dflt_real
+        surface_tension = .false.
 
         adv_n = .false.
 
@@ -615,7 +617,7 @@ contains
                 sys_size = stress_idx%end
             end if
 
-            if (.not. f_is_default(sigma)) then
+            if (surface_tension) then
                 c_idx = sys_size + 1
                 sys_size = c_idx
             end if
@@ -639,7 +641,7 @@ contains
             internalEnergies_idx%end = adv_idx%end + num_fluids
             sys_size = internalEnergies_idx%end
 
-            if (.not. f_is_default(sigma)) then
+            if (surface_tension) then
                 c_idx = sys_size + 1
                 sys_size = c_idx
             end if

src/pre_process/m_mpi_proxy.fpp

@@ -56,7 +56,7 @@ contains
             & 'perturb_flow', 'perturb_sph', 'mixlayer_vel_profile',           &
             & 'mixlayer_perturb', 'bubbles', 'polytropic', 'polydisperse',     &
             & 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt',       &
-            & 'cfl_const_dt', 'cfl_dt']
+            & 'cfl_const_dt', 'cfl_dt', 'surface_tension']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
         call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)

src/pre_process/m_start_up.fpp

@@ -141,7 +141,8 @@ contains
             sigR, sigV, dist_type, rhoRV, R0_type, &
             file_per_process, relax, relax_model, &
             palpha_eps, ptgalpha_eps, ib, num_ibs, patch_ib, &
-            sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, n_start_old
+            sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, &
+            n_start_old, surface_tension
 
         ! Inquiring the status of the pre_process.inp file
         file_loc = 'pre_process.inp'

src/simulation/m_data_output.fpp

@@ -148,15 +148,15 @@ contains
 
         ! Generating table header for the stability criteria to be outputted
         if (cfl_dt) then
-            if (any(Re_size > 0)) then
+            if (viscous) then
                 write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// &
                     'Max ==== VCFL Max ====== Rc Min ======='
             else
                 write (1, '(A)') '=========== Time-steps ============== dt ===== Time '// &
                     '============== ICFL Max ============='
             end if
         else
-            if (any(Re_size > 0)) then
+            if (viscous) then
                 write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// &
                     'Max ==== VCFL Max ====== Rc Min ======='
             else
@@ -266,7 +266,7 @@ contains
                     ! Compute mixture sound speed
                     call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0d0, c)
 
-                    if (any(Re_size > 0)) then
+                    if (viscous) then
                         call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf)
                     else
                         call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf)
@@ -282,13 +282,13 @@ contains
 #ifdef CRAY_ACC_WAR
         !$acc update host(icfl_sf)
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             !$acc update host(vcfl_sf, Rc_sf)
         end if
 
         icfl_max_loc = maxval(icfl_sf)
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             vcfl_max_loc = maxval(vcfl_sf)
             Rc_min_loc = minval(Rc_sf)
         end if
@@ -297,7 +297,7 @@ contains
         icfl_max_loc = maxval(icfl_sf)
         !$acc end kernels
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             !$acc kernels
             vcfl_max_loc = maxval(vcfl_sf)
             Rc_min_loc = minval(Rc_sf)
@@ -317,21 +317,21 @@ contains
                                                          Rc_min_glb)
         else
             icfl_max_glb = icfl_max_loc
-            if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc
-            if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc
+            if (viscous) vcfl_max_glb = vcfl_max_loc
+            if (viscous) Rc_min_glb = Rc_min_loc
         end if
 
         ! Determining the stability criteria extrema over all the time-steps
         if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb
             if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb
         end if
 
         ! Outputting global stability criteria extrema at current time-step
         if (proc_rank == 0) then
-            if (any(Re_size > 0)) then
+            if (viscous) then
                 write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') &
                     t_step, dt, t_step*dt, icfl_max_glb, &
                     vcfl_max_glb, &
@@ -348,7 +348,7 @@ contains
                 call s_mpi_abort('ICFL is greater than 1.0. Exiting ...')
             end if
 
-            if (any(Re_size > 0)) then
+            if (viscous) then
                 if (vcfl_max_glb /= vcfl_max_glb) then
                     call s_mpi_abort('VCFL is NaN. Exiting ...')
                 elseif (vcfl_max_glb > 1d0) then
@@ -1578,8 +1578,8 @@ contains
         write (3, '(A)') ''
 
         write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max
-        if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max
-        if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min
+        if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max
+        if (viscous) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min
 
         call cpu_time(run_time)
 
@@ -1615,7 +1615,7 @@ contains
         @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p))
         icfl_max = 0d0
 
-        if (any(Re_size > 0)) then
+        if (viscous) then
             @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p))
             @:ALLOCATE_GLOBAL(Rc_sf  (0:m, 0:n, 0:p))
 
@@ -1652,7 +1652,7 @@ contains
 
         ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria
         @:DEALLOCATE_GLOBAL(icfl_sf)
-        if (any(Re_size > 0)) then
+        if (viscous) then
             @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf)
         end if
 

src/simulation/m_global_parameters.fpp

@@ -163,6 +163,9 @@ module m_global_parameters
     logical :: hypoelasticity !< hypoelasticity modeling
     logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
     logical :: cu_tensor
+    logical :: viscous       !< Viscous effects
+    logical :: shear_stress  !< Shear stresses
+    logical :: bulk_stress   !< Bulk stresses
 
     !$acc declare create(chemistry)
 
@@ -183,7 +186,7 @@ module m_global_parameters
         !$acc declare create(num_dims, weno_polyn, weno_order, weno_num_stencils, num_fluids, wenojs, mapped_weno, wenoz, teno, wenoz_q)
     #:endif
 
-    !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach)
+    !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach, viscous, shear_stress, bulk_stress)
 
     logical :: relax          !< activate phase change
     integer :: relax_model    !< Relaxation model
@@ -459,7 +462,8 @@ module m_global_parameters
     !> @name Surface tension parameters
     !> @{
     real(kind(0d0)) :: sigma
-    !$acc declare create(sigma)
+    logical :: surface_tension
+    !$acc declare create(sigma, surface_tension)
     !> @}
 
     integer :: momxb, momxe
@@ -562,6 +566,9 @@ contains
         weno_flat = .true.
         riemann_flat = .true.
         rdma_mpi = .false.
+        viscous = .false.
+        shear_stress = .false.
+        bulk_stress = .false.
 
         #:if not MFC_CASE_OPTIMIZATION
             mapped_weno = .false.
@@ -650,6 +657,7 @@ contains
 
         ! Surface tension
         sigma = dflt_real
+        surface_tension = .false.
 
         ! Cuda aware MPI
         cu_tensor = .false.
@@ -888,7 +896,7 @@ contains
                     sys_size = stress_idx%end
                 end if
 
-                if (.not. f_is_default(sigma)) then
+                if (surface_tension) then
                     c_idx = sys_size + 1
                     sys_size = c_idx
                 end if
@@ -906,7 +914,7 @@ contains
                 internalEnergies_idx%end = adv_idx%end + num_fluids
                 sys_size = internalEnergies_idx%end
 
-                if (.not. f_is_default(sigma)) then
+                if (surface_tension) then
                     c_idx = sys_size + 1
                     sys_size = c_idx
                 end if
@@ -973,11 +981,14 @@ contains
                 if (fluid_pp(i)%Re(2) > 0) Re_size(2) = Re_size(2) + 1
             end do
 
-            !$acc update device(Re_size)
+            if (Re_size(1) > 0d0) shear_stress = .true.
+            if (Re_size(2) > 0d0) bulk_stress = .true.
+
+            !$acc update device(Re_size, viscous, shear_stress, bulk_stress)
 
             ! Bookkeeping the indexes of any viscous fluids and any pairs of
             ! fluids whose interface will support effects of surface tension
-            if (any(Re_size > 0)) then
+            if (viscous) then
 
                 @:ALLOCATE_GLOBAL(Re_idx(1:2, 1:maxval(Re_size)))
 
@@ -1049,7 +1060,7 @@ contains
         ! sufficient boundary conditions data as to iterate the solution in
         ! the physical computational domain from one time-step iteration to
         ! the next one
-        if (any(Re_size > 0)) then
+        if (viscous) then
             buff_size = 2*weno_polyn + 2
 !        else if (hypoelasticity) then !TODO: check if necessary
 !            buff_size = 2*weno_polyn + 2
@@ -1193,7 +1204,7 @@ contains
         ! Deallocating the variables bookkeeping the indexes of any viscous
         ! fluids and any pairs of fluids whose interfaces supported effects
         ! of surface tension
-        if (any(Re_size > 0)) then
+        if (viscous) then
             @:DEALLOCATE_GLOBAL(Re_idx)
         end if