diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 5b3f54b202..c74f421a96 100644 --- a/docs/documentation/case.md +++ b/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 | diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index 59e12cd216..8e07c857a3 100644 --- a/src/common/m_checker_common.fpp +++ b/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. diff --git a/src/post_process/m_checker.fpp b/src/post_process/m_checker.fpp index cfe8d8ad27..8d2f23b664 100644 --- a/src/post_process/m_checker.fpp +++ b/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.") diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index bdacecf1c4..33f226d0b1 100644 --- a/src/post_process/m_data_output.fpp +++ b/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 diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 4bdca155af..a20950e208 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/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 diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 796d0026f4..d53352ca93 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/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', & diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 9de428eb64..de7691c220 100644 --- a/src/post_process/m_start_up.f90 +++ b/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 diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index d35bcce4cc..a3b9a6ebd6 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/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 diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index a8be987681..ec2aa0e3a7 100644 --- a/src/simulation/m_checker.fpp +++ b/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") diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 3eeb6c0b9b..3a86ee4efd 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -1006,7 +1006,7 @@ 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 @@ -1014,16 +1014,17 @@ contains 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 diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 86c4fc5529..d223cbcbda 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/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, diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py index a8f21f0e07..252732457a 100644 --- a/toolchain/mfc/test/case.py +++ b/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',