Merge pull request #151 from lee-hyeoksu/master

Author: sbryngelson

examples/3D_turb_mixing/case.py

@@ -22,9 +22,9 @@
 Lz = 59.0
 
 # Number of grid cells
-Nx = 255
-Ny = 255
-Nz = 255
+Nx = 191
+Ny = 191
+Nz = 191
 
 # Grid spacing
 dx      = Lx/float(Nx)
@@ -75,7 +75,6 @@
     'num_fluids'                    : 1,
     'adv_alphan'                    : 'T',
     'time_stepper'                  : 3,
-    'weno_vars'                     : 2,
     'weno_order'                    : 5,
     'weno_eps'                      : 1.E-16,
     'weno_Re_flux'                  : 'T',

src/common/include/inline_conversions.fpp

@@ -32,11 +32,11 @@
 
             if (mpp_lim .and. (num_fluids > 1)) then
                 c = (1d0/gamma + 1d0)* &
-                      (pres + pi_inf)/rho
+                      (pres + pi_inf/(gamma+1d0))/rho
             else
                 c = &
                     (1d0/gamma + 1d0)* &    
-                    (pres + pi_inf)/ &
+                    (pres + pi_inf/(gamma+1d0))/ &
                     (rho*(1d0 - adv(num_fluids)))
             end if
         else 
@@ -48,6 +48,7 @@
         else
             c = sqrt(c)
         end if
+
     end subroutine s_compute_speed_of_sound
 #:enddef
 

src/common/m_variables_conversion.fpp

@@ -233,19 +233,33 @@ contains
         real(kind(0d0)), intent(OUT), target :: gamma
         real(kind(0d0)), intent(OUT), target :: pi_inf
 
+        real(kind(0d0)), dimension(num_fluids) :: alpha_rho_K, alpha_K
         real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K
 
         real(kind(0d0)), optional, intent(OUT) :: G_K
         real(kind(0d0)), optional, dimension(num_fluids), intent(IN) :: G
 
-        integer :: i
+        integer :: i, q
 
         ! Constraining the partial densities and the volume fractions within
         ! their physical bounds to make sure that any mixture variables that
         ! are derived from them result within the limits that are set by the
         ! fluids physical parameters that make up the mixture
-        ! alpha_rho_K(1) = qK_vf(i)%sf(i,j,k)
-        ! alpha_K(1)     = qK_vf(E_idx+i)%sf(i,j,k)
+        do i = 1, num_fluids
+            alpha_rho_K(i) = q_vf(i)%sf(j, k, l)
+            alpha_K(i) = q_vf(advxb + i - 1)%sf(j, k, l)
+        end do
+
+        if (mpp_lim) then
+
+            do i = 1, num_fluids
+                alpha_rho_K(i) = max(0d0, alpha_rho_K(i))
+                alpha_K(i) = min(max(0d0, alpha_K(i)), 1d0)
+            end do
+
+            alpha_K = alpha_K/max(sum(alpha_K), 1d-16)
+
+        end if
 
         ! Performing the transfer of the density, the specific heat ratio
         ! function as well as the liquid stiffness function, respectively
@@ -284,6 +298,26 @@ contains
             end if
         end if
 
+#ifdef MFC_SIMULATION
+        ! Computing the shear and bulk Reynolds numbers from species analogs
+        if (any(Re_size > 0)) then 
+            if (num_fluids == 1) then ! need to consider case with num_fluids >= 2
+                do i = 1, 2
+
+                    Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0d0
+
+                    do q = 1, Re_size(i)
+                        Re_K(i) = (1-alpha_K(Re_idx(i, q)))/fluid_pp(Re_idx(i, q))%Re(i) &
+                                + Re_K(i)
+                    end do
+
+                    Re_K(i) = 1d0/max(Re_K(i), sgm_eps)
+
+                end do
+            end if
+        end if 
+#endif
+
         ! Post process requires rho_sf/gamma_sf/pi_inf_sf to also be updated
 #ifdef MFC_POST_PROCESS
         rho_sf   (j, k, l) = rho
@@ -467,13 +501,16 @@ contains
 
     subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, &
                                                                   gamma_K, pi_inf_K, &
-                                                                  alpha_K, alpha_rho_K, k, l, r)
+                                                                  alpha_K, alpha_rho_K, Re_K, k, l, r)
 !$acc routine seq
 
         real(kind(0d0)), intent(INOUT) :: rho_K, gamma_K, pi_inf_K
 
         real(kind(0d0)), dimension(num_fluids), intent(IN) :: alpha_rho_K, alpha_K !<
             !! Partial densities and volume fractions
+
+        real(kind(0d0)), dimension(2), intent(OUT) :: Re_K
+
         integer, intent(IN) :: k, l, r
         integer :: i, j !< Generic loop iterators
 
@@ -499,6 +536,25 @@ contains
             gamma_K = gammas(1)
             pi_inf_K = pi_infs(1)
         end if
+
+        if (any(Re_size > 0)) then
+            if (num_fluids == 1) then ! need to consider case with num_fluids >= 2
+
+                do i = 1, 2
+                    Re_K(i) = dflt_real
+
+                    if (Re_size(i) > 0) Re_K(i) = 0d0
+
+                    do j = 1, Re_size(i)
+                        Re_K(i) = (1d0-alpha_K(Re_idx(i, j)))/Res(i, j) &
+                                + Re_K(i)
+                    end do
+
+                    Re_K(i) = 1d0/max(Re_K(i), sgm_eps)
+
+                end do
+            end if
+        end if
 #endif
 
     end subroutine s_convert_species_to_mixture_variables_bubbles_acc
@@ -706,7 +762,7 @@ contains
                                                                             alpha_rho_K, Re_K, j, k, l, G_K, Gs)
                         else if (bubbles) then
                             call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, &
-                                                                                alpha_K, alpha_rho_K, j, k, l)
+                                                                                alpha_K, alpha_rho_K, Re_K, j, k, l)
                         else 
                             call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, &
                                                                                 alpha_K, alpha_rho_K, Re_K, j, k, l)
@@ -1004,9 +1060,9 @@ contains
                         call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, &
                                                                         alpha_K, alpha_rho_K, Re_K, &
                                                                         j, k, l, G_K, Gs)
-!                    else if (bubbles) then
-!                        call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, &
-!                                                                pi_inf_K, alpha_K, alpha_rho_K, j, k, l)
+                   else if (bubbles) then
+                       call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, &
+                                                               pi_inf_K, alpha_K, alpha_rho_K, Re_K, j, k, l)
                     else
                         call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, &
                                                                         alpha_K, alpha_rho_K, Re_K, j, k, l)

src/pre_process/m_assign_variables.f90

@@ -271,7 +271,8 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
         end if
 
         ! Computing Mixture Variables from Original Primitive Variables
-        call s_convert_species_to_mixture_variables( &
+        ! call s_convert_species_to_mixture_variables( &
+        call s_convert_to_mixture_variables( &
             q_prim_vf, j, k, l, &
             orig_rho, &
             orig_gamma, &
@@ -346,9 +347,10 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
                 end if
             end do
         end if
-
+  
         ! Density and the specific heat ratio and liquid stiffness functions
-        call s_convert_species_to_mixture_variables( &
+        ! call s_convert_species_to_mixture_variables( &
+        call s_convert_to_mixture_variables( &
             q_prim_vf, j, k, l, &
             patch_icpp(patch_id)%rho, &
             patch_icpp(patch_id)%gamma, &
@@ -417,7 +419,8 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
         end if
 
         ! Density and the specific heat ratio and liquid stiffness functions
-        call s_convert_species_to_mixture_variables( &
+        ! call s_convert_species_to_mixture_variables( &
+        call s_convert_to_mixture_variables( &
             q_prim_vf, j, k, l, &
             patch_icpp(smooth_patch_id)%rho, &
             patch_icpp(smooth_patch_id)%gamma, &
@@ -480,7 +483,8 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
         end if
 
         ! Density and the specific heat ratio and liquid stiffness functions
-        call s_convert_species_to_mixture_variables(q_prim_vf, j, k, l, &
+        ! call s_convert_species_to_mixture_variables(q_prim_vf, j, k, l, &
+        call s_convert_to_mixture_variables(q_prim_vf, j, k, l, &
                                                     rho, gamma, pi_inf)
 
         ! Velocity

src/pre_process/m_initial_condition.fpp

@@ -381,10 +381,14 @@ contains
         ! Set fluid flow properties
         gam = 1.+1./fluid_pp(1)%gamma
         pi_inf = fluid_pp(1)%pi_inf*(gam-1.)/gam
-        rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1)
+        if (bubbles .and. num_fluids == 1) then
+            rho1 = patch_icpp(1)%alpha_rho(1)/(1d0-patch_icpp(1)%alpha(1))
+        else
+            rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1)
+        end if
         c1 = sqrt((gam*(patch_icpp(1)%pres+pi_inf))/rho1)
         mach = 1./c1
-        
+
         ! Assign mean profiles
         do j=0,n
             u_mean(j)=tanh(y_cc(j))

src/pre_process/m_patches.f90

@@ -506,7 +506,7 @@ subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf) ! ---------------------
         real(kind(0d0)) :: pi_inf, gamma, lit_gamma !< Equation of state parameters
 
         integer :: i, j !< generic loop iterators
-
+        
         pi_inf = fluid_pp(1)%pi_inf
         gamma = fluid_pp(1)%gamma
         lit_gamma = (1d0 + gamma)/gamma

src/simulation/m_bubbles.fpp

@@ -195,8 +195,9 @@ contains
                             n_tait = gammas(1)
                             B_tait = pi_infs(1)
                         end if
-
+                        
                         n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
+                        B_tait = B_tait*(n_tait-1)/n_tait ! make this the usual pi_inf
 
                         myRho = q_prim_vf(1)%sf(j, k, l)
                         myP = q_prim_vf(E_idx)%sf(j, k, l)

src/simulation/m_cbc.fpp

@@ -752,7 +752,7 @@ contains
                     end do
 
                     if (bubbles) then
-                        call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, adv, alpha_rho, 0, k, r)
+                        call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r)
 
                     else
                         call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r)
@@ -769,7 +769,6 @@ contains
 
                     ! Compute mixture sound speed
                     call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_K_sum, c)
-
                     ! ============================================================
 
                     ! First-Order Spatial Derivatives of Primitive Variables =====

src/simulation/m_data_output.fpp

@@ -255,7 +255,7 @@ contains
                     end do
 
                     if (bubbles) then
-                        call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, alpha, alpha_rho, j, k, l)
+                        call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, alpha, alpha_rho, Re, j, k, l)
                     else
                         call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, alpha, alpha_rho, Re, j, k, l)
                     end if
@@ -270,7 +270,7 @@ contains
                     end do
 
                     pres = q_prim_vf(E_idx)%sf(j, k, l)
-
+                    
                     E = gamma*pres + pi_inf + 5d-1*rho*vel_sum
 
                     H = (E + pres)/rho
@@ -302,7 +302,6 @@ contains
                         end if
 
                         if (any(Re_size > 0)) then
-
                             if (grid_geometry == 3) then
                                 vcfl_sf(j, k, l) = maxval(dt/Re) &
                                                    /min(dx(j), dy(k), fltr_dtheta)**2d0

src/simulation/m_qbmm.fpp

@@ -247,6 +247,7 @@ contains
                         n_tait = gammas(1)
                         n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
                         B_tait = pi_infs(1)
+                        B_tait = B_tait*(n_tait-1)/n_tait ! make this the usual pi_inf
                         c = n_tait*(pres + B_tait)/(rho*(1.d0 - alf))
                         if (c > 0.d0) then
                             c = DSQRT(c)