Hyper edits

docs/documentation/case.md

@@ -508,6 +508,7 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running.
 | `vel_wrt(i)`         | Logical | Add the $i$-direction velocity to the database	  |
 | `E_wrt`              | Logical | Add the total energy to the database	 |
 | `pres_wrt`           | Logical | Add the pressure to the database	|
+| `tau_wrt`            | Logical | Add the elastic stresses to the database |
 | `alpha_wrt(i)`       | Logical | Add the volume fraction of fluid $i$ to the database	|
 | `gamma_wrt`          | Logical | Add the specific heat ratio function to the database	|
 | `heat_ratio_wrt`     | Logical | Add the specific heat ratio to the database	|

src/common/m_checker_common.fpp

@@ -161,23 +161,12 @@ contains
     !> Checks constraints on the elasticity parameters.
         !! Called by s_check_inputs_common for all three stages
     subroutine s_check_inputs_elasticity
-        @:PROHIBIT((hypoelasticity .or. hyperelasticity) .and. model_eqns == 1, &
-            "Elasticity does not work for model_eqns = 1")
-        @:PROHIBIT((hypoelasticity .or. hyperelasticity) .and. model_eqns > 3, &
-            "Elasticity works only for model_eqns 2 and 3")
+        @:PROHIBIT(elasticity .and. .not. (model_eqns == 2 .or. model_eqns == 3))
         #:for X in ['x', 'y', 'z']
             #:for BOUND in ['beg', 'end']
-                @:PROHIBIT(hyperelasticity .and. ((bc_${X}$%${BOUND}$ .lt. -3)), &
-                    "bc_${X}$%${BOUND}$ is not supported")
+                @:PROHIBIT(hyperelasticity .and. bc_${X}$%${BOUND}$ /= dflt_int .and. (bc_${X}$%${BOUND}$ < -3))
             #:endfor
         #:endfor
-
-#ifdef MFC_SIMULATION
-        @:PROHIBIT(elasticity .and. fd_order /= 4)
-        @:PROHIBIT(hyperelasticity .and. hyper_model .le. 0, &
-            "Set the hyper_model in the input file")
-#endif
-
     end subroutine s_check_inputs_elasticity
 
     !> Checks constraints on dimensionality and the number of cells for the grid.

src/post_process/m_checker.fpp

@@ -33,6 +33,7 @@ contains
         call s_check_inputs_volume_fraction
         call s_check_inputs_vorticity
         call s_check_inputs_schlieren
+        call s_check_inputs_elasticity
         call s_check_inputs_surface_tension
         call s_check_inputs_no_flow_variables
 
@@ -129,6 +130,14 @@ contains
         end do
     end subroutine s_check_inputs_schlieren
 
+    !> Checks constraints on elasticity parameters
+    subroutine s_check_inputs_elasticity
+        @:PROHIBIT(.not. (hypoelasticity .or. hyperelasticity) .and. tau_wrt)
+        ! Note: 'elasticity' variable isn't initialized yet; use (hypoelasticity .or. hyperelasticity) instead
+        @:PROHIBIT(.not. hyperelasticity .and. kymograph)
+        @:PROHIBIT(format == 2 .and. kymograph, 'Binary output format does not support kymograph')
+    end subroutine
+
     !> Checks constraints on surface tension parameters (cf_wrt and sigma)
     subroutine s_check_inputs_surface_tension
         @:PROHIBIT(cf_wrt .and. .not. surface_tension, &
@@ -138,7 +147,7 @@ contains
     !> Checks constraints on the absence of flow variables
     subroutine s_check_inputs_no_flow_variables
         @:PROHIBIT(.not. any([ &
-            (/rho_wrt, E_wrt, pres_wrt, gamma_wrt, heat_ratio_wrt, pi_inf_wrt, &
+            (/rho_wrt, E_wrt, pres_wrt, tau_wrt, gamma_wrt, heat_ratio_wrt, pi_inf_wrt, &
               pres_inf_wrt, cons_vars_wrt, prim_vars_wrt, c_wrt, schlieren_wrt/), &
             alpha_rho_wrt, mom_wrt, vel_wrt, flux_wrt, alpha_wrt, omega_wrt]), &
             "None of the flow variables have been selected for post-process. Exiting.")

src/post_process/m_data_output.fpp

@@ -345,7 +345,7 @@ contains
             if (pres_wrt .or. prim_vars_wrt) dbvars = dbvars + 1
 
             ! Elastic stresses
-            if (hypoelasticity) dbvars = dbvars + (num_dims*(num_dims + 1))/2
+            if (tau_wrt .or. prim_vars_wrt) dbvars = dbvars + (num_dims*(num_dims + 1))/2
 
             ! Volume fraction(s)
             if ((model_eqns == 2) .or. (model_eqns == 3)) then

src/post_process/m_global_parameters.fpp

@@ -219,6 +219,7 @@ module m_global_parameters
     logical, dimension(3) :: flux_wrt
     logical :: E_wrt
     logical :: pres_wrt
+    logical :: tau_wrt
     logical, dimension(num_fluids_max) :: alpha_wrt
     logical :: gamma_wrt
     logical :: heat_ratio_wrt
@@ -385,6 +386,7 @@ contains
         file_per_process = .false.
         E_wrt = .false.
         pres_wrt = .false.
+        tau_wrt = .false.
         alpha_wrt = .false.
         gamma_wrt = .false.
         heat_ratio_wrt = .false.
@@ -649,22 +651,23 @@ contains
             end if
         end if
 
-        if (hypoelasticity .or. hyperelasticity) elasticity = .true.
+        elasticity = hypoelasticity .or. hyperelasticity
 
         if (elasticity) then
             stress_idx%beg = sys_size + 1
             stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
             ! number of distinct stresses is 1 in 1D, 3 in 2D, 6 in 3D
             sys_size = stress_idx%end
-            if (hyperelasticity) then
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                tensor_size = num_dims**2 + 1
-            end if
+        end if
+
+        if (hyperelasticity) then
+            xi_idx%beg = sys_size + 1
+            xi_idx%end = sys_size + num_dims
+            ! adding three more equations for the \xi field and the elastic energy
+            sys_size = xi_idx%end + 1
+            ! number of entries in the symmetric btensor plus the jacobian
+            b_size = (num_dims*(num_dims + 1))/2 + 1
+            tensor_size = num_dims**2 + 1
         end if
 
         if (chemistry) then

src/post_process/m_mpi_proxy.fpp

@@ -163,7 +163,7 @@ contains
 
         #:for VAR in [ 'cyl_coord', 'mpp_lim', 'mixture_err',                  &
             & 'alt_soundspeed', 'hypoelasticity', 'parallel_io', 'rho_wrt',    &
-            & 'E_wrt', 'pres_wrt', 'gamma_wrt', 'sim_data', 'kymograph',        &
+            & 'E_wrt', 'pres_wrt', 'tau_wrt', 'gamma_wrt', 'sim_data', 'kymograph',        &
             & 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', &
             & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles_euler', 'qbmm',   &
             & 'polytropic', 'polydisperse', 'file_per_process', 'relax', 'cf_wrt',     &

src/post_process/m_start_up.f90

@@ -72,7 +72,7 @@ subroutine s_read_input_file
             hypoelasticity, G, &
             chem_wrt_Y, chem_wrt_T, avg_state, &
             alpha_rho_wrt, rho_wrt, mom_wrt, vel_wrt, &
-            E_wrt, pres_wrt, alpha_wrt, gamma_wrt, &
+            E_wrt, pres_wrt, tau_wrt, alpha_wrt, gamma_wrt, &
             heat_ratio_wrt, pi_inf_wrt, pres_inf_wrt, &
             cons_vars_wrt, prim_vars_wrt, c_wrt, &
             omega_wrt, qm_wrt, schlieren_wrt, schlieren_alpha, &
@@ -361,35 +361,30 @@ subroutine s_save_data(t_step, varname, pres, c, H)
         end if
 
         ! Adding the elastic shear stresses to the formatted database file
-        if (elasticity) then
-            if (prim_vars_wrt) then
-                do i = 1, stress_idx%end - stress_idx%beg + 1
-                    q_sf = q_prim_vf(i - 1 + stress_idx%beg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
-                    write (varname, '(A,I0)') 'tau', i
-                    call s_write_variable_to_formatted_database_file(varname, t_step)
+        if (elasticity .and. (tau_wrt .or. prim_vars_wrt)) then
+            do i = 1, stress_idx%end - stress_idx%beg + 1
+                q_sf = q_prim_vf(i - 1 + stress_idx%beg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                write (varname, '(A,I0)') 'tau', i
+                call s_write_variable_to_formatted_database_file(varname, t_step)
 
-                    varname(:) = ' '
-                end do
-            end if
+                varname(:) = ' '
+            end do
         end if
 
-        if (hyperelasticity) then
-            if (prim_vars_wrt) then
-                do i = 1, xiend - xibeg + 1
-                    q_sf = q_prim_vf(i - 1 + xibeg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
-                    write (varname, '(A,I0)') 'xi', i
-                    call s_write_variable_to_formatted_database_file(varname, t_step)
-
-                    varname(:) = ' '
-                end do
-
-                q_sf = q_prim_vf(xiend + 1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
-
-                write (varname, '(A,I0)') 'vonMises'
+        if (hyperelasticity .and. (tau_wrt .or. prim_vars_wrt)) then
+            do i = 1, xiend - xibeg + 1
+                q_sf = q_prim_vf(i - 1 + xibeg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                write (varname, '(A,I0)') 'xi', i
                 call s_write_variable_to_formatted_database_file(varname, t_step)
+
                 varname(:) = ' '
+            end do
 
-            end if
+            q_sf = q_prim_vf(xiend + 1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+
+            write (varname, '(A,I0)') 'vonMises'
+            call s_write_variable_to_formatted_database_file(varname, t_step)
+            varname(:) = ' '
         end if
 
         ! Adding the pressure to the formatted database file

src/pre_process/m_global_parameters.fpp

@@ -720,22 +720,24 @@ contains
             end if
         end if
 
-        if (hypoelasticity .or. hyperelasticity) elasticity = .true.
+        elasticity = hypoelasticity .or. hyperelasticity
+
         if (elasticity) then
             ! creates stress indices for both hypo and hyperelasticity
             stress_idx%beg = sys_size + 1
             stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
             ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
             sys_size = stress_idx%end
-            if (hyperelasticity) then
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                tensor_size = num_dims**2 + 1
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-            end if
+        end if
+
+        if (hyperelasticity) then
+            ! number of entries in the symmetric btensor plus the jacobian
+            b_size = (num_dims*(num_dims + 1))/2 + 1
+            tensor_size = num_dims**2 + 1
+            xi_idx%beg = sys_size + 1
+            xi_idx%end = sys_size + num_dims
+            ! adding three more equations for the \xi field and the elastic energy
+            sys_size = xi_idx%end + 1
         end if
 
         if (chemistry) then

src/simulation/m_checker.fpp

@@ -30,6 +30,7 @@ contains
         call s_check_inputs_time_stepping
         call s_check_inputs_model_eqns
         call s_check_inputs_acoustic_src
+        call s_check_inputs_elasticity
         call s_check_inputs_bubbles_euler
         call s_check_inputs_bubbles_lagrange
         call s_check_inputs_adapt_dt
@@ -256,6 +257,12 @@ contains
 
     end subroutine s_check_inputs_acoustic_src
 
+    !> Checks constraints on elasticity parameters
+    subroutine s_check_inputs_elasticity
+        @:PROHIBIT(hyperelasticity .and. hyper_model == dflt_int)
+        @:PROHIBIT(elasticity .and. fd_order /= 4)
+    end subroutine
+
     !> Checks constraints on bubble parameters
     subroutine s_check_inputs_bubbles_euler
         @:PROHIBIT(bubbles_euler .and. bubbles_lagrange, "Activate only one of the bubble subgrid models")

src/simulation/m_global_parameters.fpp

@@ -1006,24 +1006,25 @@ contains
         end if
         ! END: Volume Fraction Model
 
-        if (hypoelasticity .or. hyperelasticity) elasticity = .true.
+        elasticity = hypoelasticity .or. hyperelasticity
 
         if (elasticity) then
             ! creates stress indices for both hypo and hyperelasticity
             stress_idx%beg = sys_size + 1
             stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
             ! number of distinct stresses is 1 in 1D, 3 in 2D, 6 in 3D
             sys_size = stress_idx%end
-            if (hyperelasticity) then
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                ! storing the jacobian in the last entry
-                tensor_size = num_dims**2 + 1
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-            end if
+        end if
+
+        if (hyperelasticity) then
+            ! number of entries in the symmetric btensor plus the jacobian
+            b_size = (num_dims*(num_dims + 1))/2 + 1
+            ! storing the jacobian in the last entry
+            tensor_size = num_dims**2 + 1
+            xi_idx%beg = sys_size + 1
+            xi_idx%end = sys_size + num_dims
+            ! adding three more equations for the \xi field and the elastic energy
+            sys_size = xi_idx%end + 1
         end if
 
         if (chemistry) then

toolchain/mfc/run/case_dicts.py

@@ -363,6 +363,7 @@ def analytic(self):
     'flux_wrt': ParamType.LOG,
     'E_wrt': ParamType.LOG,
     'pres_wrt': ParamType.LOG,
+    'tau_wrt': ParamType.LOG,
     'alpha_wrt': ParamType.LOG,
     'kappa_wrt': ParamType.LOG,
     'gamma_wrt': ParamType.LOG,

toolchain/mfc/test/case.py

@@ -220,7 +220,8 @@ def create_directory(self):
         'prim_vars_wrt': 'T', 'alpha_rho_wrt(1)': 'T',
         'rho_wrt'      : 'T', 'mom_wrt(1)'      : 'T',
         'vel_wrt(1)'   : 'T', 'E_wrt'           : 'T',
-        'pres_wrt'     : 'T', 'alpha_wrt(1)'    : 'T',
+        'pres_wrt'     : 'T',
+        'tau_wrt'      : 'T', 'alpha_wrt(1)'    : 'T',
         'gamma_wrt'    : 'T', 'heat_ratio_wrt'  : 'T',
         'pi_inf_wrt'   : 'T', 'pres_inf_wrt'    : 'T',
         'c_wrt'        : 'T',