diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 2732ba69eb..ec82b5bbd8 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -536,7 +536,6 @@ See @ref equations "Equations" for the mathematical models these parameters cont | `surface_tension` | Logical | Activate surface tension | | `viscous` | Logical | Activate viscosity | | `hypoelasticity` | Logical | Activate hypoelasticity* | -| `pre_stress` | Logical | Enable pre-stress initialization for hypoelasticity | | `igr` | Logical | Enable solution via information geometric regularization (IGR) \cite Cao24 | | `igr_order` | Integer | Order of reconstruction for IGR [3,5] | | `alf_factor` | Real | Alpha factor for IGR entropic pressure (default 10) | @@ -1124,7 +1123,6 @@ Note: For relativistic flow, the conservative and primitive densities are differ | Parameter | Type | Description | | ---: | :---: | :--- | | `hypoelasticity` | Logical | Enable hypoelasticity simulation | -| `hyperelasticity` | Logical | Enable hyperelasticity simulation | | `cont_damage` | Logical | Enable continuum damage model | | `tau_star` | Real | Threshold stress for continuum damage model | | `cont_damage_s` | Real | Power `s` for continuum damage model | diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index 27d4962c4e..0c120da36f 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -573,34 +573,6 @@ where \f$\mathbf{l} = \nabla \mathbf{u}\f$ is the velocity gradient and \f$\math This adds 6 additional transport equations in 3D (symmetric stress tensor: \f$\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e\f$). -### 7.2 Hyperelastic Model (`hyperelasticity = .true.`) (\cite Kamrin12; \cite Wilfong26 Sec. 4.1.6) - -**Source:** `src/simulation/m_hyperelastic.fpp` - -**Reference map evolution:** - -\f[\frac{\partial (\rho\,\boldsymbol{\xi})}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\xi} \otimes \mathbf{u}) = 0\f] - -**Deformation gradient from reference map:** - -\f[\mathbf{F} = (\nabla \boldsymbol{\xi})^{-1}\f] - -**Left Cauchy-Green tensor:** - -\f[\mathbf{b} = \mathbf{F}\,\mathbf{F}^T\f] - -**Neo-Hookean Cauchy stress:** - -\f[\boldsymbol{\tau}^e = \frac{G}{J}\left(\mathbf{b} - \frac{\text{tr}(\mathbf{b})}{3}\,\mathbf{I}\right)\f] - -where \f$J = \det(\mathbf{F})\f$. - -**Hyperelastic energy:** - -\f[e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text{tr}(\mathbf{b})\f] - ---- - ## 8. Phase Change (`relax = .true.`) (\cite Wilfong26 Sec. 4.1.3) **Source:** `src/common/m_phase_change.fpp` diff --git a/docs/module_categories.json b/docs/module_categories.json index a48f8b023a..9100364765 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -28,7 +28,6 @@ "m_bubbles_EL", "m_bubbles_EL_kernels", "m_qbmm", - "m_hyperelastic", "m_hypoelastic", "m_phase_change", "m_chemistry", diff --git a/src/common/m_boundary_primitives.fpp b/src/common/m_boundary_primitives.fpp index 2fce625617..9732aa3907 100644 --- a/src/common/m_boundary_primitives.fpp +++ b/src/common/m_boundary_primitives.fpp @@ -135,16 +135,12 @@ contains q_T_sf%sf(-j, k, l) = q_T_sf%sf(j - 1, k, l) end if - if (elasticity) then + if (hypoelasticity) then do i = 1, shear_BC_flip_num q_prim_vf(shear_BC_flip_indices(1, i))%sf(-j, k, l) = -q_prim_vf(shear_BC_flip_indices(1, & & i))%sf(j - 1, k, l) end do end if - - if (hyperelasticity) then - q_prim_vf(eqn_idx%xi%beg)%sf(-j, k, l) = -q_prim_vf(eqn_idx%xi%beg)%sf(j - 1, k, l) - end if end do if (qbmm .and. .not. polytropic .and. present(pb_in) .and. present(mv_in)) then @@ -173,16 +169,12 @@ contains q_T_sf%sf(m + j, k, l) = q_T_sf%sf(m - (j - 1), k, l) end if - if (elasticity) then + if (hypoelasticity) then do i = 1, shear_BC_flip_num q_prim_vf(shear_BC_flip_indices(1, i))%sf(m + j, k, l) = -q_prim_vf(shear_BC_flip_indices(1, & & i))%sf(m - (j - 1), k, l) end do end if - - if (hyperelasticity) then - q_prim_vf(eqn_idx%xi%beg)%sf(m + j, k, l) = -q_prim_vf(eqn_idx%xi%beg)%sf(m - (j - 1), k, l) - end if end do if (qbmm .and. .not. polytropic .and. present(pb_in) .and. present(mv_in)) then do i = 1, nb @@ -212,16 +204,12 @@ contains q_T_sf%sf(k, -j, l) = q_T_sf%sf(k, j - 1, l) end if - if (elasticity) then + if (hypoelasticity) then do i = 1, shear_BC_flip_num q_prim_vf(shear_BC_flip_indices(2, i))%sf(k, -j, l) = -q_prim_vf(shear_BC_flip_indices(2, i))%sf(k, & & j - 1, l) end do end if - - if (hyperelasticity) then - q_prim_vf(eqn_idx%xi%beg + 1)%sf(k, -j, l) = -q_prim_vf(eqn_idx%xi%beg + 1)%sf(k, j - 1, l) - end if end do if (qbmm .and. .not. polytropic .and. present(pb_in) .and. present(mv_in)) then @@ -250,16 +238,12 @@ contains q_T_sf%sf(k, n + j, l) = q_T_sf%sf(k, n - (j - 1), l) end if - if (elasticity) then + if (hypoelasticity) then do i = 1, shear_BC_flip_num q_prim_vf(shear_BC_flip_indices(2, i))%sf(k, n + j, l) = -q_prim_vf(shear_BC_flip_indices(2, & & i))%sf(k, n - (j - 1), l) end do end if - - if (hyperelasticity) then - q_prim_vf(eqn_idx%xi%beg + 1)%sf(k, n + j, l) = -q_prim_vf(eqn_idx%xi%beg + 1)%sf(k, n - (j - 1), l) - end if end do if (qbmm .and. .not. polytropic .and. present(pb_in) .and. present(mv_in)) then @@ -290,16 +274,12 @@ contains q_T_sf%sf(k, l, -j) = q_T_sf%sf(k, l, j - 1) end if - if (elasticity) then + if (hypoelasticity) then do i = 1, shear_BC_flip_num q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, -j) = -q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, & & l, j - 1) end do end if - - if (hyperelasticity) then - q_prim_vf(eqn_idx%xi%end)%sf(k, l, -j) = -q_prim_vf(eqn_idx%xi%end)%sf(k, l, j - 1) - end if end do if (qbmm .and. .not. polytropic .and. present(pb_in) .and. present(mv_in)) then @@ -328,16 +308,12 @@ contains q_T_sf%sf(k, l, p + j) = q_T_sf%sf(k, l, p - (j - 1)) end if - if (elasticity) then + if (hypoelasticity) then do i = 1, shear_BC_flip_num q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, p + j) = -q_prim_vf(shear_BC_flip_indices(3, & & i))%sf(k, l, p - (j - 1)) end do end if - - if (hyperelasticity) then - q_prim_vf(eqn_idx%xi%end)%sf(k, l, p + j) = -q_prim_vf(eqn_idx%xi%end)%sf(k, l, p - (j - 1)) - end if end do if (qbmm .and. .not. polytropic .and. present(pb_in) .and. present(mv_in)) then diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index e52054befa..872c04b390 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -140,7 +140,6 @@ module m_derived_types type(idx_bounds_info) :: adv !< Volume fractions (advection equations) type(idx_bounds_info) :: bub !< Bubble equation range (beg/end only) type(idx_bounds_info) :: stress !< Stress tensor components - type(idx_bounds_info) :: xi !< Reference map equations type(idx_bounds_info) :: B !< Magnetic field components type(idx_bounds_info) :: int_en !< Internal energy equations type(idx_bounds_info) :: species !< Chemistry species equations diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 4e65b887e0..991d54e030 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -44,10 +44,8 @@ module m_global_parameters_common !> @name Annotations of the structure of the state and flux vectors in terms of the size and configuration of the system of !! equations !> @{ - integer :: sys_size !< Number of unknowns in system of equations - type(eqn_idx_info) :: eqn_idx !< All conserved-variable equation index ranges and scalars - integer :: b_size !< Number of elements in the symmetric b tensor, plus one - integer :: tensor_size !< Number of elements in the full tensor plus one + integer :: sys_size !< Number of unknowns in system of equations + type(eqn_idx_info) :: eqn_idx !< All conserved-variable equation index ranges and scalars !> @} !> @name Chemistry modeling (Fypp compile-time constant; same value in all targets) @@ -55,9 +53,8 @@ module m_global_parameters_common logical, parameter :: chemistry = .${chemistry}$. !> @} - !> @name Elasticity and shear stress state (identical across all three executables) + !> @name Hypoelastic shear stress state (identical across all three executables) !> @{ - logical :: elasticity !< elasticity modeling, true for hyper or hypo integer :: shear_num !< Number of shear stress components integer, dimension(3) :: shear_indices !< Indices of the stress components that represent shear stress integer :: shear_BC_flip_num !< Number of shear stress components to reflect for boundary conditions @@ -65,7 +62,7 @@ module m_global_parameters_common !> @} #ifdef MFC_SIMULATION - $:GPU_DECLARE(create='[sys_size, eqn_idx, b_size, tensor_size]') + $:GPU_DECLARE(create='[sys_size, eqn_idx]') $:GPU_DECLARE(create='[shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices]') ! Device residency for namelist/case-opt state declared above via the generated ! includes: declare directives must live in the declaring module (Cray ftn rejects @@ -77,7 +74,7 @@ module m_global_parameters_common $:GPU_DECLARE(create='[muscl_eps]') $:GPU_DECLARE(create='[mpp_lim, model_eqns, mixture_err, alt_soundspeed]') $:GPU_DECLARE(create='[avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity]') - $:GPU_DECLARE(create='[hyperelasticity, elasticity, low_Mach]') + $:GPU_DECLARE(create='[low_Mach]') $:GPU_DECLARE(create='[cont_damage, hyper_cleaning]') $:GPU_DECLARE(create='[relax, relax_model, palpha_eps, ptgalpha_eps]') $:GPU_DECLARE(create='[down_sample]') @@ -125,9 +122,8 @@ module m_global_parameters_common contains - !> Initialize equation-index state (eqn_idx, sys_size, b_size, tensor_size) from the namelist parameters. This is the shared - !! skeleton: it covers the model_eqns dispatch, all eqn_idx field assignments, and the elasticity/surface-tension/chemistry - !! extensions. + !> Initialize equation-index state (eqn_idx and sys_size) from the namelist parameters. This is the shared skeleton: it covers + !! the model_eqns dispatch, all eqn_idx field assignments, and the hypoelastic/surface-tension/chemistry extensions. !! !! @param nmom_in Number of carried moments per R0 location (per-target: pre/post pass an !! integer variable; sim passes its integer parameter nmom = 6). Used only in the 5eq @@ -255,8 +251,7 @@ contains end if if (model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) then - if (hypoelasticity .or. hyperelasticity) then - elasticity = .true. + if (hypoelasticity) then eqn_idx%stress%beg = sys_size + 1 eqn_idx%stress%end = sys_size + (num_dims*(num_dims + 1))/2 if (cyl_coord) eqn_idx%stress%end = eqn_idx%stress%end + 1 @@ -283,16 +278,6 @@ contains end if end if - if (hyperelasticity) then - ! number of entries in the symmetric b tensor plus the jacobian - b_size = (num_dims*(num_dims + 1))/2 + 1 - tensor_size = num_dims**2 + 1 - eqn_idx%xi%beg = sys_size + 1 - eqn_idx%xi%end = sys_size + num_dims - ! adding equations for the xi field and the elastic energy - sys_size = eqn_idx%xi%end + 1 - end if - if (surface_tension) then eqn_idx%c = sys_size + 1 sys_size = eqn_idx%c @@ -408,10 +393,6 @@ contains relax = .false. relax_model = dflt_int hypoelasticity = .false. - hyperelasticity = .false. - elasticity = .false. - b_size = dflt_int - tensor_size = dflt_int cont_damage = .false. hyper_cleaning = .false. diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index fe17b89c27..a230416b5e 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -515,7 +515,7 @@ contains if (model_eqns /= model_eqns_4eq) then #ifdef MFC_SIMULATION ! If in simulation, use acc mixture subroutines - if (elasticity) then + if (hypoelasticity) then call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & & Re_K, G_K, Gs_vc) else @@ -524,7 +524,7 @@ contains end if #else ! If pre-processing, use non acc mixture subroutines - if (elasticity) then + if (hypoelasticity) then call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, rho_K, gamma_K, pi_inf_K, qv_K, Re_K, G_K, & & fluid_pp(:)%G) else @@ -725,7 +725,7 @@ contains end do end if - if (elasticity) then + if (hypoelasticity) then $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%stress%beg, eqn_idx%stress%end qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/rho_K @@ -749,13 +749,6 @@ contains end do end if - if (hyperelasticity) then - $:GPU_LOOP(parallelism='[seq]') - do i = eqn_idx%xi%beg, eqn_idx%xi%end - qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/rho_K - end do - end if - if (.not. igr .or. num_fluids > 1) then $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%adv%beg, eqn_idx%adv%end @@ -978,7 +971,7 @@ contains end do end if - if (elasticity) then + if (hypoelasticity) then ! adding the elastic contribution Multiply \tau to \rho \tau do i = eqn_idx%stress%beg, eqn_idx%stress%end q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) @@ -1001,14 +994,6 @@ contains end do end if - ! using \rho xi as the conservative formulation stated in Kamrin et al. JFM 2022 - if (hyperelasticity) then - ! Multiply \xi to \rho \xi - do i = eqn_idx%xi%beg, eqn_idx%xi%end - q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) - end do - end if - if (surface_tension) then q_cons_vf(eqn_idx%c)%sf(j, k, l) = q_prim_vf(eqn_idx%c)%sf(j, k, l) end if @@ -1098,7 +1083,7 @@ contains end do pres_K = qK_prim_vf(j, k, l, eqn_idx%E) - if (elasticity) then + if (hypoelasticity) then call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & & Re_K, G_K, Gs_vc) else diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index b24941c46c..28c65fa3fd 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -407,7 +407,7 @@ impure subroutine s_read_parallel_conservative_data(t_step, m_MOK, n_MOK, p_MOK, str_MOK = int(name_len, MPI_OFFSET_KIND) NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) - if (bubbles_euler .or. elasticity .or. mhd) then + if (bubbles_euler .or. hypoelasticity .or. mhd) then do i = 1, sys_size var_MOK = int(i, MPI_OFFSET_KIND) call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 7e1c804700..c9729081a5 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -15,7 +15,7 @@ module m_global_parameters use m_helper_basic use m_thermochem, only: species_names use m_constants, only: format_silo, precision_single - ! Shared state: generated_decls, num_dims, num_vels, sys_size, eqn_idx, b_size, tensor_size, chemistry, elasticity, shear_* + ! Shared state: generated_decls, num_dims, num_vels, sys_size, eqn_idx, chemistry, shear_* use m_global_parameters_common implicit none @@ -75,7 +75,7 @@ module m_global_parameters !> @name Simulation Algorithm Parameters !> @{ - ! sys_size, elasticity, b_size, tensor_size, chemistry, eqn_idx: in m_global_parameters_common + ! sys_size, chemistry, eqn_idx: in m_global_parameters_common !> @} !> @name Annotations of the structure, i.e. the organization, of the state vectors @@ -152,7 +152,7 @@ contains integer :: i !< Generic loop iterator - ! Shared defaults (case_dir, m/n/p, cyl_coord, cfl flags, model_eqns, elasticity, BC blocks, + ! Shared defaults (case_dir, m/n/p, cyl_coord, cfl flags, model_eqns, BC blocks, ! recon/weno/muscl/num_fluids/igr/mhd/relativity under case-opt guard, Tait EOS, bubble flags, ! IB flags, parallel I/O flags, fft_wrt) @@ -331,7 +331,7 @@ contains ! (guard matches the original site: inside the 5-equation branch) if (model_eqns == model_eqns_5eq .and. qbmm) nmom = 6 - ! Populate eqn_idx, sys_size, b_size, tensor_size, elasticity, shear_* (shared logic) + ! Populate eqn_idx, sys_size, shear_* (shared logic) call s_initialize_eqn_idx(nmom, nb) ! post-only: 6eq alf is a dummy (no void fraction in 6eq) diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 3b294ff075..2c5c9182f4 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -442,7 +442,7 @@ contains end do end if - if (elasticity) then + if (hypoelasticity) then do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 if (prim_vars_wrt) then write (varname, '(A,I0)') 'tau', i @@ -452,15 +452,6 @@ contains end do end if - if (hyperelasticity) then - do i = 1, eqn_idx%xi%end - eqn_idx%xi%beg + 1 - if (prim_vars_wrt) then - write (varname, '(A,I0)') 'xi', i - call s_write_field(varname, t_step, q_prim_vf(i - 1 + eqn_idx%xi%beg), x_beg, x_end, y_beg, y_end, z_beg, z_end) - end if - end do - end if - if (cont_damage) then write (varname, '(A)') 'damage_state' call s_write_field(varname, t_step, q_cons_vf(eqn_idx%damage), x_beg, x_end, y_beg, y_end, z_beg, z_end) diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index 1325eaaf63..20340576a6 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -231,8 +231,6 @@ contains real(wp) :: orig_qv real(wp) :: muR, muV real(wp) :: R3bar - real(wp) :: rcoord, theta, phi, xi_sph - real(wp), dimension(3) :: xi_cart real(wp) :: Ys(1:num_species) real(stp), dimension(sys_size) :: orig_prim_vf !< Vector to hold original values of cell for smoothing purposes integer :: i @@ -376,35 +374,13 @@ contains end if end if - if (elasticity) then + if (hypoelasticity) then do i = 1, (eqn_idx%stress%end - eqn_idx%stress%beg) + 1 q_prim_vf(i + eqn_idx%stress%beg - 1)%sf(j, k, & & l) = (eta*patch_icpp(patch_id)%tau_e(i) + (1._wp - eta)*orig_prim_vf(i + eqn_idx%stress%beg - 1)) end do end if - if (hyperelasticity) then - if (pre_stress) then ! pre stressed initial condition in spatial domain - rcoord = sqrt((x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)) - theta = atan2(y_cc(k), x_cc(j)) - phi = atan2(sqrt(x_cc(j)**2 + y_cc(k)**2), z_cc(l)) - ! spherical coord, assuming Rmax=1 - xi_sph = (rcoord**3 - R0ref**3 + 1._wp)**(1._wp/3._wp) - xi_cart(1) = xi_sph*sin(phi)*cos(theta) - xi_cart(2) = xi_sph*sin(phi)*sin(theta) - xi_cart(3) = xi_sph*cos(phi) - else - xi_cart(1) = x_cc(j) - xi_cart(2) = y_cc(k) - xi_cart(3) = z_cc(l) - end if - - ! assigning the reference map to the q_prim vector field - do i = 1, num_dims - q_prim_vf(i + eqn_idx%xi%beg - 1)%sf(j, k, l) = eta*xi_cart(i) + (1._wp - eta)*orig_prim_vf(i + eqn_idx%xi%beg - 1) - end do - end if - if (mpp_lim .and. bubbles_euler) then ! adjust volume fractions, according to modeled gas void fraction alf_sum%sf = 0._wp diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 0376529957..54c6e84b1a 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -710,7 +710,6 @@ contains call write_range(eqn_idx%bub%beg, eqn_idx%bub%end, " Bubbles") call write_range(eqn_idx%stress%beg, eqn_idx%stress%end, " Stress") call write_range(eqn_idx%int_en%beg, eqn_idx%int_en%end, " Internal Energies") - call write_range(eqn_idx%xi%beg, eqn_idx%xi%end, " Reference Map") call write_range(eqn_idx%B%beg, eqn_idx%B%end, " Magnetic Field") call write_range(eqn_idx%c, eqn_idx%c, " Color Function") call write_range(eqn_idx%species%beg, eqn_idx%species%end, " Chemistry") diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index bbf2444803..7ef66f2ac5 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -13,7 +13,7 @@ module m_global_parameters use m_derived_types ! Definitions of the derived types use m_helper_basic ! Functions to compare floating point numbers - ! Shared state: generated_decls, sys_size, eqn_idx, b_size, tensor_size, chemistry, elasticity, shear_* + ! Shared state: generated_decls, sys_size, eqn_idx, chemistry, shear_* use m_global_parameters_common implicit none @@ -43,7 +43,7 @@ module m_global_parameters type(bounds_info) :: x_domain_glb, y_domain_glb, z_domain_glb ! Simulation Algorithm Parameters - ! sys_size, eqn_idx, b_size, tensor_size, chemistry, elasticity, shear_*: in m_global_parameters_common + ! sys_size, eqn_idx, chemistry, shear_*: in m_global_parameters_common ! weno_polyn, muscl_polyn, num_dims, num_vels: in m_global_parameters_common ! Annotations of the structure, i.e. the organization, of the state vectors type(qbmm_idx_info) :: qbmm_idx !< QBMM moment index mappings. @@ -117,7 +117,7 @@ contains integer :: i !< Generic loop operator - ! Shared defaults (case_dir, m/n/p, cyl_coord, cfl flags, model_eqns, elasticity, BC blocks, + ! Shared defaults (case_dir, m/n/p, cyl_coord, cfl flags, model_eqns, BC blocks, ! recon/weno/muscl/num_fluids/igr/mhd/relativity under case-opt guard, Tait EOS, bubble flags, ! IB flags, parallel I/O flags, fft_wrt) @@ -181,8 +181,6 @@ contains palpha_eps = dflt_real ptgalpha_eps = dflt_real igr_order = dflt_int - pre_stress = .false. - precision = 2 viscous = .false. mixlayer_vel_profile = .false. @@ -449,7 +447,7 @@ contains ! (guards match the original site: 5-equation bubbles with 4-node qbmm) if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. qbmm .and. nnode == 4) nmom = 6 - ! Populate eqn_idx, sys_size, b_size, tensor_size, elasticity, shear_* (shared logic) + ! Populate eqn_idx, sys_size, shear_* (shared logic) call s_initialize_eqn_idx(nmom, nb) ! Per-target (pre_process): qbmm_idx allocations and fills diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index ad621c796f..f6b58f4e39 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1213,7 +1213,7 @@ contains end if ! Computing/Sharing necessary state variables - if (elasticity) then + if (hypoelasticity) then call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv, Re, G_local, & & fluid_pp(:)%G) else @@ -1225,7 +1225,7 @@ contains dyn_p = 0.5_wp*rho*dot_product(vel, vel) - if (elasticity) then + if (hypoelasticity) then if (cont_damage) then damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k, l) G_local = G_local*max((1._wp - damage_state), 0._wp) @@ -1242,7 +1242,7 @@ contains if (model_eqns == model_eqns_4eq) then lit_gamma = gammas(1) - else if (elasticity) then + else if (hypoelasticity) then tau_e(1) = q_cons_vf(eqn_idx%stress%end)%sf(j - 2, k, l)/rho end if @@ -1329,7 +1329,7 @@ contains dyn_p = 0.5_wp*rho*dot_product(vel, vel) - if (elasticity) then + if (hypoelasticity) then if (cont_damage) then damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k - 2, l) G_local = G_local*max((1._wp - damage_state), 0._wp) @@ -1346,7 +1346,7 @@ contains if (model_eqns == model_eqns_4eq) then lit_gamma = gs_min(1) - else if (elasticity) then + else if (hypoelasticity) then do s = 1, 3 tau_e(s) = q_cons_vf(s)%sf(j - 2, k - 2, l)/rho end do @@ -1416,7 +1416,7 @@ contains end do end if - if (elasticity) then + if (hypoelasticity) then if (cont_damage) then damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k - 2, l - 2) G_local = G_local*max((1._wp - damage_state), 0._wp) @@ -1466,7 +1466,7 @@ contains end if end if - if (elasticity) then + if (hypoelasticity) then do s = 1, (num_dims*(num_dims + 1))/2 tmp = tau_e(s) call s_mpi_allreduce_sum(tmp, tau_e(s)) @@ -1508,7 +1508,7 @@ contains write (i + 30, '(6X,10F24.8)') nondim_time, rho, vel(1), vel(2), pres, alf, nR(1), nRdot(1), R(1), & & Rdot(1) #:endif - else if (elasticity) then + else if (hypoelasticity) then #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8)') nondim_time, rho, & & vel(1), vel(2), pres, tau_e(1), tau_e(2), tau_e(3) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 5db4b0f72e..6ec9c3107b 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -14,8 +14,7 @@ module m_global_parameters use m_derived_types use m_helper_basic - ! Shared state: generated_decls, generated_case_opt_decls, sys_size, eqn_idx, b_size, tensor_size, chemistry, elasticity, - ! shear_* + ! Shared state: generated_decls, generated_case_opt_decls, sys_size, eqn_idx, chemistry, shear_* use m_global_parameters_common ! $:USE_GPU_MODULE() @@ -67,8 +66,7 @@ module m_global_parameters logical :: cfl_dt ! Simulation Algorithm Parameters generated_case_opt_decls.fpp: now in m_global_parameters_common - integer :: hyper_model !< hyperelasticity solver algorithm - ! elasticity, chemistry: in m_global_parameters_common + ! chemistry: in m_global_parameters_common logical :: shear_stress !< Shear stresses logical :: bulk_stress !< Bulk stresses logical :: bodyForces @@ -92,7 +90,6 @@ module m_global_parameters integer :: cpu_start, cpu_end, cpu_rate - $:GPU_DECLARE(create='[hyper_model]') $:GPU_DECLARE(create='[shear_stress, bulk_stress]') logical :: bc_io @@ -138,7 +135,7 @@ module m_global_parameters type(mpi_io_levelset_norm_var), public :: MPI_IO_levelsetnorm_DATA real(wp), allocatable, dimension(:,:), public :: MPI_IO_DATA_lag_bubbles - ! sys_size, eqn_idx, b_size, tensor_size: in m_global_parameters_common (GPU_DECLARE there too) + ! sys_size and eqn_idx: in m_global_parameters_common (GPU_DECLARE there too) type(qbmm_idx_info) :: qbmm_idx !< QBMM moment index mappings (allocatable; GPU-managed separately). ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). Stands for "InDices With INTerior". @@ -321,7 +318,7 @@ contains integer :: i, j !< Generic loop iterator - ! Shared defaults (case_dir, m/n/p, cyl_coord, cfl flags, model_eqns, elasticity, BC blocks, + ! Shared defaults (case_dir, m/n/p, cyl_coord, cfl flags, model_eqns, BC blocks, ! recon/weno/muscl/num_fluids/igr/mhd/relativity under case-opt guard, Tait EOS, bubble flags, ! IB flags, parallel I/O flags, fft_wrt) @@ -391,7 +388,6 @@ contains int_comp = 0 ic_eps = dflt_ic_eps ic_beta = dflt_ic_beta - hyper_model = dflt_int rdma_mpi = .false. shear_stress = .false. bulk_stress = .false. @@ -736,12 +732,12 @@ contains Re_size = 0 Re_size_max = 0 - ! Populate eqn_idx, sys_size, b_size, tensor_size, elasticity, shear_* (shared logic) + ! Populate eqn_idx, sys_size, shear_* (shared logic) call s_initialize_eqn_idx(nmom, nb) ! sim-only: GPU update for shear state after s_initialize_eqn_idx populated it if (model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) then - if (hypoelasticity .or. hyperelasticity) then + if (hypoelasticity) then $:GPU_UPDATE(device='[shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices]') end if end if @@ -906,7 +902,7 @@ contains if (ib) allocate (MPI_IO_IB_DATA%var%sf(0:m,0:n,0:p)) - if (elasticity .or. mhd .or. probe_wrt .or. ib .or. bubbles_lagrange) then + if (hypoelasticity .or. mhd .or. probe_wrt .or. ib .or. bubbles_lagrange) then fd_number = max(1, fd_order/2) end if @@ -930,14 +926,11 @@ contains end if $:GPU_UPDATE(device='[sys_size, buff_size, eqn_idx, adv_n, adap_dt, pi_fac, adap_dt_tol, adap_dt_max_iters]') - $:GPU_UPDATE(device='[b_size, tensor_size]') - $:GPU_UPDATE(device='[cfl_target, m, n, p]') $:GPU_UPDATE(device='[alt_soundspeed, acoustic_source, num_source]') $:GPU_UPDATE(device='[dt, sys_size, buff_size, pref, rhoref, eqn_idx, mpp_lim, bubbles_euler, hypoelasticity, & - & alt_soundspeed, avg_state, model_eqns, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, & - & hyperelasticity, hyper_model, elasticity, low_Mach]') + & alt_soundspeed, avg_state, model_eqns, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, low_Mach]') $:GPU_UPDATE(device='[Bx0]') diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp deleted file mode 100644 index e0017a3685..0000000000 --- a/src/simulation/m_hyperelastic.fpp +++ /dev/null @@ -1,270 +0,0 @@ -!> -!! @file -!! @brief Contains module m_hyperelastic - -#:include 'macros.fpp' - -!> @brief Computes the left Cauchy--Green deformation tensor and hyperelastic stress source terms - -module m_hyperelastic - - use m_derived_types - use m_global_parameters - use m_variables_conversion - use m_finite_differences - - implicit none - - private; public :: s_hyperelastic_rmt_stress_update, s_initialize_hyperelastic_module, s_finalize_hyperelastic_module - - ! The btensor at the cell-interior Gaussian quadrature points. These tensor is needed to be calculated once and make the code - ! DRY. - type(vector_field) :: btensor - $:GPU_DECLARE(create='[btensor]') - - real(wp), allocatable, dimension(:,:) :: fd_coeff_x_hyper - real(wp), allocatable, dimension(:,:) :: fd_coeff_y_hyper - real(wp), allocatable, dimension(:,:) :: fd_coeff_z_hyper - $:GPU_DECLARE(create='[fd_coeff_x_hyper, fd_coeff_y_hyper, fd_coeff_z_hyper]') - real(wp), allocatable, dimension(:) :: Gs_hyper - $:GPU_DECLARE(create='[Gs_hyper]') - -contains - - !> Initialize the hyperelastic module - impure subroutine s_initialize_hyperelastic_module - - integer :: i !< generic iterator - - @:ALLOCATE(btensor%vf(1:b_size)) - do i = 1, b_size - @:ALLOCATE(btensor%vf(i)%sf(0:m, 0:n, 0:p)) - end do - @:ACC_SETUP_VFs(btensor) - - @:ALLOCATE(Gs_hyper(1:num_fluids)) - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - Gs_hyper(i) = fluid_pp(i)%G - end do - $:GPU_UPDATE(device='[Gs_hyper]') - - @:ALLOCATE(fd_coeff_x_hyper(-fd_number:fd_number, 0:m)) - if (n > 0) then - @:ALLOCATE(fd_coeff_y_hyper(-fd_number:fd_number, 0:n)) - end if - if (p > 0) then - @:ALLOCATE(fd_coeff_z_hyper(-fd_number:fd_number, 0:p)) - end if - - ! Computing centered finite difference coefficients - call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_hyper, buff_size, fd_number, fd_order) - $:GPU_UPDATE(device='[fd_coeff_x_hyper]') - if (n > 0) then - call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_hyper, buff_size, fd_number, fd_order) - $:GPU_UPDATE(device='[fd_coeff_y_hyper]') - end if - if (p > 0) then - call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_hyper, buff_size, fd_number, fd_order) - $:GPU_UPDATE(device='[fd_coeff_z_hyper]') - end if - - end subroutine s_initialize_hyperelastic_module - - !> Compute the left Cauchy-Green deformation tensor and update the hyperelastic stress - subroutine s_hyperelastic_rmt_stress_update(q_cons_vf, q_prim_vf) - - type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - - #:if USING_AMD - real(wp), dimension(10) :: tensora, tensorb - #:else - real(wp), dimension(tensor_size) :: tensora, tensorb - #:endif - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_k, alpha_rho_k - #:else - real(wp), dimension(num_fluids) :: alpha_k, alpha_rho_k - #:endif - real(wp), dimension(2) :: Re - real(wp) :: rho, gamma, pi_inf, qv - real(wp) :: G_local - integer :: j, k, l, i, r - - $:GPU_PARALLEL_LOOP(collapse=3, & - & private='[i, j, k, l, alpha_K, alpha_rho_K, rho, gamma, pi_inf, qv, G_local, Re, tensora, tensorb]') - do l = 0, p - do k = 0, n - do j = 0, m - call s_compute_species_fraction(q_cons_vf, j, k, l, alpha_rho_k, alpha_k) - - ! If in simulation, use acc mixture subroutines - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, alpha_rho_k, Re, G_local, & - & Gs_hyper) - rho = max(rho, sgm_eps) - G_local = max(G_local, sgm_eps) - - if (G_local > verysmall) then - $:GPU_LOOP(parallelism='[seq]') - do i = 1, tensor_size - tensora(i) = 0._wp - end do - ! STEP 1: computing the grad_xi tensor using finite differences grad_xi definition / organization number for - ! the tensor 1-3: dxix_dx, dxiy_dx, dxiz_dx 4-6 : dxix_dy, dxiy_dy, dxiz_dy 7-9 : dxix_dz, dxiy_dz, dxiz_dz - $:GPU_LOOP(parallelism='[seq]') - do r = -fd_number, fd_number - ! derivatives in the x-direction - tensora(1) = tensora(1) + q_prim_vf(eqn_idx%xi%beg)%sf(j + r, k, l)*fd_coeff_x_hyper(r, j) - tensora(2) = tensora(2) + q_prim_vf(eqn_idx%xi%beg + 1)%sf(j + r, k, l)*fd_coeff_x_hyper(r, j) - tensora(3) = tensora(3) + q_prim_vf(eqn_idx%xi%end)%sf(j + r, k, l)*fd_coeff_x_hyper(r, j) - ! derivatives in the y-direction - tensora(4) = tensora(4) + q_prim_vf(eqn_idx%xi%beg)%sf(j, k + r, l)*fd_coeff_y_hyper(r, k) - tensora(5) = tensora(5) + q_prim_vf(eqn_idx%xi%beg + 1)%sf(j, k + r, l)*fd_coeff_y_hyper(r, k) - tensora(6) = tensora(6) + q_prim_vf(eqn_idx%xi%end)%sf(j, k + r, l)*fd_coeff_y_hyper(r, k) - ! derivatives in the z-direction - tensora(7) = tensora(7) + q_prim_vf(eqn_idx%xi%beg)%sf(j, k, l + r)*fd_coeff_z_hyper(r, l) - tensora(8) = tensora(8) + q_prim_vf(eqn_idx%xi%beg + 1)%sf(j, k, l + r)*fd_coeff_z_hyper(r, l) - tensora(9) = tensora(9) + q_prim_vf(eqn_idx%xi%end)%sf(j, k, l + r)*fd_coeff_z_hyper(r, l) - end do - ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse - tensorb(1) = tensora(5)*tensora(9) - tensora(6)*tensora(8) - tensorb(2) = -(tensora(2)*tensora(9) - tensora(3)*tensora(8)) - tensorb(3) = tensora(2)*tensora(6) - tensora(3)*tensora(5) - tensorb(4) = -(tensora(4)*tensora(9) - tensora(6)*tensora(7)) - tensorb(5) = tensora(1)*tensora(9) - tensora(3)*tensora(7) - tensorb(6) = -(tensora(1)*tensora(6) - tensora(4)*tensora(3)) - tensorb(7) = tensora(4)*tensora(8) - tensora(5)*tensora(7) - tensorb(8) = -(tensora(1)*tensora(8) - tensora(2)*tensora(7)) - tensorb(9) = tensora(1)*tensora(5) - tensora(2)*tensora(4) - - ! STEP 2b: computing the determinant of the grad_xi tensor - tensorb(tensor_size) = tensora(1)*(tensora(5)*tensora(9) - tensora(6)*tensora(8)) - tensora(2)*(tensora(4) & - & *tensora(9) - tensora(6)*tensora(7)) + tensora(3)*(tensora(4)*tensora(8) - tensora(5)*tensora(7)) - - if (tensorb(tensor_size) > verysmall) then - ! STEP 2c: computing the inverse of grad_xi tensor = F tensorb is the adjoint, tensora becomes F - $:GPU_LOOP(parallelism='[seq]') - do i = 1, tensor_size - 1 - tensora(i) = tensorb(i)/tensorb(tensor_size) - end do - - ! STEP 2d: computing the J = det(F) = 1/det(\grad{\xi}) - tensorb(tensor_size) = 1._wp/tensorb(tensor_size) - - ! STEP 3: computing F transpose F - tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensora(3)**2 - tensorb(5) = tensora(4)**2 + tensora(5)**2 + tensora(6)**2 - tensorb(9) = tensora(7)**2 + tensora(8)**2 + tensora(9)**2 - tensorb(2) = tensora(1)*tensora(4) + tensora(2)*tensora(5) + tensora(3)*tensora(6) - tensorb(3) = tensora(1)*tensora(7) + tensora(2)*tensora(8) + tensora(3)*tensora(9) - tensorb(6) = tensora(4)*tensora(7) + tensora(5)*tensora(8) + tensora(6)*tensora(9) - ! STEP 4: update the btensor, this is consistent with Riemann solvers - #:for BIJ, TXY in [(1,1),(2,2),(3,5),(4,3),(5,6),(6,9)] - btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) - #:endfor - ! store the determinant at the last entry of the btensor - btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) - ! STEP 5a: updating the Cauchy stress primitive scalar field - if (hyper_model == 1) then - call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G_local, j, k, l) - else if (hyper_model == 2) then - call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G_local, j, k, l) - end if - ! STEP 5b: updating the pressure field - q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, & - & l) - G_local*q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l)/gamma - ! STEP 5c: updating the Cauchy stress conservative scalar field - $:GPU_LOOP(parallelism='[seq]') - do i = 1, b_size - 1 - q_cons_vf(eqn_idx%stress%beg + i - 1)%sf(j, k, & - & l) = rho*q_prim_vf(eqn_idx%stress%beg + i - 1)%sf(j, k, l) - end do - end if - end if - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_hyperelastic_rmt_stress_update - - !> Compute the neo-Hookean Cauchy stress from the left Cauchy-Green tensor - subroutine s_neoHookean_cauchy_solver(btensor_in, q_prim_vf, G_param, j, k, l) - - $:GPU_ROUTINE(parallelism='[seq]') - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(b_size), intent(inout) :: btensor_in - real(wp), intent(in) :: G_param - integer, intent(in) :: j, k, l - real(wp) :: trace - real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i !< Generic loop iterators - ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor_in(1)%sf(j, k, l) + btensor_in(3)%sf(j, k, l) + btensor_in(6)%sf(j, k, l) - - ! Deviatoric left Cauchy-Green tensor: dev(b) = b - (tr(b)/3)*I - #:for IJ in [1,3,6] - btensor_in(${IJ}$)%sf(j, k, l) = btensor_in(${IJ}$)%sf(j, k, l) - f13*trace - #:endfor - ! dividing by the jacobian for neo-Hookean model setting the tensor to the stresses for riemann solver - $:GPU_LOOP(parallelism='[seq]') - do i = 1, b_size - 1 - q_prim_vf(eqn_idx%stress%beg + i - 1)%sf(j, k, l) = G_param*btensor_in(i)%sf(j, k, l)/btensor_in(b_size)%sf(j, k, l) - end do - ! First invariant strain energy: W = G/2 * (I1 - 3), neo-Hookean model - q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l) = 0.5_wp*(trace - 3.0_wp)/btensor_in(b_size)%sf(j, k, l) - - end subroutine s_neoHookean_cauchy_solver - - !> Compute the Mooney-Rivlin Cauchy stress from the left Cauchy-Green tensor - subroutine s_Mooney_Rivlin_cauchy_solver(btensor_in, q_prim_vf, G_param, j, k, l) - - $:GPU_ROUTINE(parallelism='[seq]') - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(b_size), intent(inout) :: btensor_in - real(wp), intent(in) :: G_param - integer, intent(in) :: j, k, l - real(wp) :: trace - real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i !< Generic loop iterators - ! TODO: Make 1D and 2D capable - trace = btensor_in(1)%sf(j, k, l) + btensor_in(3)%sf(j, k, l) + btensor_in(6)%sf(j, k, l) - - ! Deviatoric left Cauchy-Green tensor: dev(b) = b - (tr(b)/3)*I - btensor_in(1)%sf(j, k, l) = btensor_in(1)%sf(j, k, l) - f13*trace - btensor_in(3)%sf(j, k, l) = btensor_in(3)%sf(j, k, l) - f13*trace - btensor_in(6)%sf(j, k, l) = btensor_in(6)%sf(j, k, l) - f13*trace - - ! dividing by the jacobian for neo-Hookean model setting the tensor to the stresses for riemann solver - $:GPU_LOOP(parallelism='[seq]') - do i = 1, b_size - 1 - q_prim_vf(eqn_idx%stress%beg + i - 1)%sf(j, k, l) = G_param*btensor_in(i)%sf(j, k, l)/btensor_in(b_size)%sf(j, k, l) - end do - ! First invariant strain energy: W = G/2 * (I1 - 3), neo-Hookean model - q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l) = 0.5_wp*(trace - 3.0_wp)/btensor_in(b_size)%sf(j, k, l) - - end subroutine s_Mooney_Rivlin_cauchy_solver - - !> Finalize the hyperelastic module - impure subroutine s_finalize_hyperelastic_module() - - integer :: i !< iterator - ! Deallocating memory - - do i = 1, b_size - @:DEALLOCATE(btensor%vf(i)%sf) - end do - @:DEALLOCATE(Gs_hyper) - @:DEALLOCATE(fd_coeff_x_hyper) - if (n > 0) then - @:DEALLOCATE(fd_coeff_y_hyper) - if (p > 0) then - @:DEALLOCATE(fd_coeff_z_hyper) - end if - end if - - end subroutine s_finalize_hyperelastic_module - -end module m_hyperelastic diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 624784ee38..d2e598ed22 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -295,7 +295,7 @@ contains if (model_eqns /= model_eqns_4eq) then ! If in simulation, use acc mixture subroutines - if (elasticity) then + if (hypoelasticity) then call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K, & & G_K, Gs) else diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 337aa8a3ed..5a443fc124 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -23,7 +23,6 @@ module m_rhs use m_bubbles_EL use m_qbmm use m_hypoelastic - use m_hyperelastic use m_acoustic_src use m_viscous use m_ibm @@ -564,10 +563,6 @@ contains call nvtxEndRange end if - call nvtxStartRange("RHS-ELASTIC") - if (hyperelasticity) call s_hyperelastic_rmt_stress_update(q_cons_qp%vf, q_prim_qp%vf) - call nvtxEndRange - if (cfl_dt) then if (mytime >= t_stop) return else diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index ddfc38cfd9..e222642c43 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -76,7 +76,6 @@ contains real(wp) :: G_L, G_R real(wp) :: damage_L, damage_R real(wp), dimension(2) :: Re_L, Re_R - real(wp), dimension(3) :: xi_field_L, xi_field_R real(wp) :: rho_avg real(wp) :: H_avg real(wp) :: qv_avg @@ -114,15 +113,15 @@ contains #:set SF = lambda offs: COORDS.format(STENCIL_IDX = SV + offs) if (norm_dir == ${NORM_DIR}$) then $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, & - & tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, Ys_L, Ys_R, xi_field_L, & - & xi_field_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, & - & h_avg_2, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, zcoef, vel_L_tmp, vel_R_tmp, & - & rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, & - & T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, & - & gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, & - & damage_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, & - & vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, & - & flux_tau_R]', copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]') + & tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, Ys_L, Ys_R, Cp_iL, & + & Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, c_fast, & + & pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, zcoef, vel_L_tmp, vel_R_tmp, rho_L, rho_R, & + & pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, & + & Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, & + & gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, & + & rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, & + & Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R]', & + & copyin='[norm_dir]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -362,11 +361,6 @@ contains & /rho_R), & & vel_L(dir_idx(1)) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1))) & & /rho_L)) - else if (hyperelasticity) then - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L), & - & vel_R(dir_idx(1)) - sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R), & - & vel_L(dir_idx(1)) + sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L)) else s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index 7ee6eca2b6..2f1ab5a542 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -102,23 +102,15 @@ contains real(wp), dimension(nb) :: pbw_L, pbw_R #:endif - real(wp) :: alpha_L_sum, alpha_R_sum, nbub_L, nbub_R - real(wp) :: ptilde_L, ptilde_R - real(wp) :: PbwR3Lbar, PbwR3Rbar - real(wp) :: R3Lbar, R3Rbar - real(wp) :: R3V2Lbar, R3V2Rbar - real(wp), dimension(6) :: tau_e_L, tau_e_R - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: xi_field_L, xi_field_R - #:else - real(wp), dimension(num_dims) :: xi_field_L, xi_field_R - #:endif - real(wp) :: G_L, G_R + real(wp) :: alpha_L_sum, alpha_R_sum, nbub_L, nbub_R + real(wp) :: ptilde_L, ptilde_R + real(wp) :: PbwR3Lbar, PbwR3Rbar + real(wp) :: R3Lbar, R3Rbar + real(wp) :: R3V2Lbar, R3V2Rbar real(wp) :: vel_L_rms, vel_R_rms, vel_avg_rms real(wp) :: vel_L_tmp, vel_R_tmp real(wp) :: rho_Star, E_Star, p_Star, p_K_Star, vel_K_star real(wp) :: pres_SL, pres_SR, Ms_L, Ms_R - real(wp) :: flux_ene_e real(wp) :: zcoef, pcorr !< low Mach number correction integer :: i, j, k, l, q !< Generic loop iterators integer :: Re_size_loc1, Re_size_loc2 !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU @@ -145,15 +137,14 @@ contains ! 6-equation model (model_eqns=3): separate phasic internal energies $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, & & alpha_rho_L, alpha_rho_R, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, & - & Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, tau_e_L, tau_e_R, flux_ene_e, xi_field_L, & - & xi_field_R, pcorr, zcoef, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, & - & Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, & - & Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, & - & qv_R, qv_avg, c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, & - & vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, & - & alpha_L_sum, alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, vel_K_star, s_L, s_R, & - & s_M, s_P, s_S, xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, xi_PP]', & - & firstprivate='[Re_size_loc1, Re_size_loc2]') + & Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, pcorr, zcoef, rho_L, rho_R, pres_L, pres_R, E_L, & + & E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, & + & MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, & + & pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, rho_avg, H_avg, c_avg, gamma_avg, & + & ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, & + & Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, & + & vel_K_star, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, & + & xi_PP]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -235,32 +226,6 @@ contains E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R - ! Hyperelastic stress contribution: strain energy added to total energy - if (hyperelasticity) then - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - xi_field_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%xi%beg - 1 + i) - xi_field_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%xi%beg - 1 + i) - end do - G_L = 0._wp; G_R = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - ! Mixture left and right shear modulus - G_L = G_L + alpha_L(i)*Gs_rs(i) - G_R = G_R + alpha_R(i)*Gs_rs(i) - end do - ! Elastic contribution to energy if G large enough - if (G_L > verysmall .and. G_R > verysmall) then - E_L = E_L + G_L*qL_prim_rsx_vf(${SF('')}$, eqn_idx%xi%end + 1) - E_R = E_R + G_R*qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%xi%end + 1) - end if - $:GPU_LOOP(parallelism='[seq]') - do i = 1, b_size - 1 - tau_e_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%stress%beg - 1 + i) - tau_e_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%stress%beg - 1 + i) - end do - end if - H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -291,27 +256,11 @@ contains ! COMPUTING THE DIRECT WAVE SPEEDS if (wave_speeds == wave_speeds_direct) then - if (elasticity) then - ! Elastic wave speed, Rodriguez et al. JCP (2019) - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1) & - & ))/rho_L), & - & vel_R(dir_idx(1)) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) & - & + tau_e_R(dir_idx_tau(1)))/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1) & - & ))/rho_R), & - & vel_L(dir_idx(1)) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) & - & + tau_e_L(dir_idx_tau(1)))/rho_L)) - s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) & - & + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) & - & *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R & - & - vel_R(dir_idx(1)))) - else - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) & - & - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L & - & - vel_L(dir_idx(1))) - rho_R*(s_R - vel_R(dir_idx(1)))) - end if + s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) + s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) + s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) & + & - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) & + & - rho_R*(s_R - vel_R(dir_idx(1)))) else if (wave_speeds == wave_speeds_pressure) then pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(dir_idx(1)) - vel_R(dir_idx(1)))) @@ -387,25 +336,6 @@ contains ! ENERGY FLUX. f = u*(E-\sigma), q = E, q_star = \xi*E+(s-u)(\rho s_star - \sigma/(s-u)) flux_rsx_vf(${SF('')}$, eqn_idx%E) = (E_star + p_Star)*vel_K_Star + (s_M/s_L)*(s_P/s_R)*pcorr*s_S - ! ELASTICITY. Elastic shear stress additions for the momentum and energy flux - if (elasticity) then - flux_ene_e = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - ! MOMENTUM ELASTIC FLUX. - flux_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(i)) = flux_rsx_vf(${SF('')}$, & - & eqn_idx%cont%end + dir_idx(i)) - xi_M*tau_e_L(dir_idx_tau(i)) & - & - xi_P*tau_e_R(dir_idx_tau(i)) - ! ENERGY ELASTIC FLUX. - flux_ene_e = flux_ene_e - xi_M*(vel_L(dir_idx(i))*tau_e_L(dir_idx_tau(i)) & - & + s_M*(xi_L*((s_S - vel_L(i))*(tau_e_L(dir_idx_tau(i)) & - & /(s_L - vel_L(i)))))) - xi_P*(vel_R(dir_idx(i)) & - & *tau_e_R(dir_idx_tau(i)) + s_P*(xi_R*((s_S - vel_R(i)) & - & *(tau_e_R(dir_idx_tau(i))/(s_R - vel_R(i)))))) - end do - flux_rsx_vf(${SF('')}$, eqn_idx%E) = flux_rsx_vf(${SF('')}$, eqn_idx%E) + flux_ene_e - end if - ! VOLUME FRACTION FLUX. $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%adv%beg, eqn_idx%adv%end @@ -445,17 +375,6 @@ contains flux_src_rsx_vf(${SF('')}$, eqn_idx%adv%beg) = vel_src_rsx_vf(${SF('')}$, dir_idx(1)) - ! Hyperelastic reference map flux for material deformation tracking - if (hyperelasticity) then - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - flux_rsx_vf(${SF('')}$, & - & eqn_idx%xi%beg - 1 + i) = xi_M*(s_S/(s_L - s_S))*(s_L*rho_L*xi_field_L(i) & - & - rho_L*vel_L(dir_idx(1))*xi_field_L(i)) + xi_P*(s_S/(s_R - s_S)) & - & *(s_R*rho_R*xi_field_R(i) - rho_R*vel_R(dir_idx(1))*xi_field_R(i)) - end do - end if - ! COLOR FUNCTION FLUX if (surface_tension) then flux_rsx_vf(${SF('')}$, eqn_idx%c) = (xi_M*qL_prim_rsx_vf(${SF('')}$, & @@ -508,7 +427,7 @@ contains & nbub_L, nbub_R, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, & & T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, & & Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, & - & G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, & + & rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, & & vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, & & alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, vel_K_star, s_L, s_R, s_M, s_P, s_S, & & xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, xi_PP, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, & @@ -1100,8 +1019,8 @@ contains & xi_R, xi_L_m1, xi_R_m1, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, & & alpha_L, alpha_R, alpha_rho_L, alpha_rho_R, alpha_lim_L, alpha_lim_R, s_L, s_R, s_S, & & vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, & - & Gamma_iR, Cp_iL, Cp_iR, tau_e_L, tau_e_R, xi_field_L, xi_field_R, Yi_avg, Phi_avg, & - & h_iL, h_iR, h_avg_2, G_L, G_R]', copyin='[is1, is2, is3]', firstprivate='[Re_size_loc1, Re_size_loc2]') + & Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]', copyin='[is1, is2, & + & is3]', firstprivate='[Re_size_loc1, Re_size_loc2]') do l = ${Z_BND}$%beg, ${Z_BND}$%end do k = ${Y_BND}$%beg, ${Y_BND}$%end do j = ${X_BND}$%beg, ${X_BND}$%end @@ -1228,33 +1147,6 @@ contains H_R = (E_R + pres_R)/rho_R end if - ! Hyperelastic stress contribution: strain energy added to total energy - if (hyperelasticity) then - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - xi_field_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%xi%beg - 1 + i) - xi_field_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%xi%beg - 1 + i) - end do - G_L = 0._wp - G_R = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - ! Mixture left and right shear modulus - G_L = G_L + alpha_L(i)*Gs_rs(i) - G_R = G_R + alpha_R(i)*Gs_rs(i) - end do - ! Elastic contribution to energy if G large enough - if (G_L > verysmall .and. G_R > verysmall) then - E_L = E_L + G_L*qL_prim_rsx_vf(${SF('')}$, eqn_idx%xi%end + 1) - E_R = E_R + G_R*qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%xi%end + 1) - end if - $:GPU_LOOP(parallelism='[seq]') - do i = 1, b_size - 1 - tau_e_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%stress%beg - 1 + i) - tau_e_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%stress%beg - 1 + i) - end do - end if - H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1287,27 +1179,11 @@ contains end if if (wave_speeds == wave_speeds_direct) then - if (elasticity) then - ! Elastic wave speed, Rodriguez et al. JCP (2019) - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1) & - & ))/rho_L), & - & vel_R(dir_idx(1)) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) & - & + tau_e_R(dir_idx_tau(1)))/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1) & - & ))/rho_R), & - & vel_L(dir_idx(1)) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) & - & + tau_e_L(dir_idx_tau(1)))/rho_L)) - s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) & - & + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) & - & *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R & - & - vel_R(dir_idx(1)))) - else - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) & - & - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L & - & - vel_L(dir_idx(1))) - rho_R*(s_R - vel_R(dir_idx(1)))) - end if + s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) + s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) + s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) & + & - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) & + & - rho_R*(s_R - vel_R(dir_idx(1)))) else if (wave_speeds == wave_speeds_pressure) then pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(dir_idx(1)) - vel_R(dir_idx(1)))) @@ -1378,25 +1254,6 @@ contains & - vel_R(dir_idx(1)))*(rho_R*s_S + pres_R/(s_R - vel_R(dir_idx(1)))))) + (s_M/s_L) & & *(s_P/s_R)*pcorr*s_S - ! ELASTICITY. Elastic shear stress additions for the momentum and energy flux - if (elasticity) then - flux_ene_e = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - ! MOMENTUM ELASTIC FLUX. - flux_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(i)) = flux_rsx_vf(${SF('')}$, & - & eqn_idx%cont%end + dir_idx(i)) - xi_M*tau_e_L(dir_idx_tau(i)) & - & - xi_P*tau_e_R(dir_idx_tau(i)) - ! ENERGY ELASTIC FLUX. - flux_ene_e = flux_ene_e - xi_M*(vel_L(dir_idx(i))*tau_e_L(dir_idx_tau(i)) & - & + s_M*(xi_L*((s_S - vel_L(i))*(tau_e_L(dir_idx_tau(i)) & - & /(s_L - vel_L(i)))))) - xi_P*(vel_R(dir_idx(i)) & - & *tau_e_R(dir_idx_tau(i)) + s_P*(xi_R*((s_S - vel_R(i)) & - & *(tau_e_R(dir_idx_tau(i))/(s_R - vel_R(i)))))) - end do - flux_rsx_vf(${SF('')}$, eqn_idx%E) = flux_rsx_vf(${SF('')}$, eqn_idx%E) + flux_ene_e - end if - ! VOLUME FRACTION FLUX. $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%adv%beg, eqn_idx%adv%end @@ -1420,17 +1277,6 @@ contains & + xi_P*qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%c)*(vel_R(dir_idx(1)) + s_P*xi_R_m1) end if - ! Hyperelastic reference map flux for material deformation tracking - if (hyperelasticity) then - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - flux_rsx_vf(${SF('')}$, & - & eqn_idx%xi%beg - 1 + i) = xi_M*(s_S/(s_L - s_S))*(s_L*rho_L*xi_field_L(i) & - & - rho_L*vel_L(dir_idx(1))*xi_field_L(i)) + xi_P*(s_S/(s_R - s_S)) & - & *(s_R*rho_R*xi_field_R(i) - rho_R*vel_R(dir_idx(1))*xi_field_R(i)) - end do - end if - flux_src_rsx_vf(${SF('')}$, eqn_idx%adv%beg) = vel_src_rsx_vf(${SF('')}$, dir_idx(1)) if (chemistry) then diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index 087bfc0aaf..5fe5499c04 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -107,7 +107,7 @@ contains $:GPU_UPDATE(device='[is1, is2, is3]') - if (elasticity) then + if (hypoelasticity) then if (norm_dir == 1) then dir_idx_tau = (/1, 2, 4/) else if (norm_dir == 2) then diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 3e380a8081..6ba61d86cc 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -96,7 +96,7 @@ contains call s_compute_species_fraction(q_prim_vf, j, k, l, alpha_rho, alpha) - if (elasticity) then + if (hypoelasticity) then call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, G_local, Gs) else call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) @@ -128,11 +128,6 @@ contains E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_sum + qv end if - ! Adjust energy for hyperelasticity - if (hyperelasticity) then - E = E + G_local*q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l) - end if - H = (E + pres)/rho end subroutine s_compute_enthalpy diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 97a1c76afd..488c858535 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -28,7 +28,6 @@ module m_start_up use m_qbmm use m_derived_variables use m_hypoelastic - use m_hyperelastic use m_phase_change use m_viscous use m_bubbles_EE @@ -234,7 +233,7 @@ contains end if end do - if (bubbles_euler .or. elasticity) then + if (bubbles_euler .or. hypoelasticity) then ! Read pb and mv for non-polytropic qbmm if (qbmm .and. .not. polytropic) then do i = 1, nb @@ -394,7 +393,7 @@ contains WP_MOK = int(storage_size(0._stp)/8, MPI_OFFSET_KIND) MOK = int(1._wp, MPI_OFFSET_KIND) - if (bubbles_euler .or. elasticity) then + if (bubbles_euler .or. hypoelasticity) then do i = 1, sys_size var_MOK = int(i, MPI_OFFSET_KIND) @@ -456,7 +455,7 @@ contains WP_MOK = int(storage_size(0._stp)/8, MPI_OFFSET_KIND) MOK = int(1._wp, MPI_OFFSET_KIND) - if (bubbles_euler .or. elasticity) then + if (bubbles_euler .or. hypoelasticity) then do i = 1, sys_size var_MOK = int(i, MPI_OFFSET_KIND) disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) @@ -922,7 +921,6 @@ contains if (bubbles_lagrange) call s_initialize_bubbles_EL_module(q_cons_ts(1)%vf, bc_type) if (hypoelasticity) call s_initialize_hypoelastic_module() - if (hyperelasticity) call s_initialize_hyperelastic_module() end subroutine s_initialize_modules @@ -1080,7 +1078,6 @@ contains call s_finalize_time_steppers_module() if (hypoelasticity) call s_finalize_hypoelastic_module() - if (hyperelasticity) call s_finalize_hyperelastic_module() call s_finalize_derived_variables_module() call s_finalize_data_output_module() call s_finalize_rhs_module() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index cc0c71a66c..daf1c2a734 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -17,7 +17,6 @@ module m_time_steppers use m_bubbles_EE use m_bubbles_EL use m_ibm - use m_hyperelastic use m_mpi_proxy use m_boundary_common use m_helper @@ -250,7 +249,7 @@ contains end do end if - if (elasticity) then + if (hypoelasticity) then do i = eqn_idx%stress%beg, eqn_idx%stress%end @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) @@ -258,14 +257,6 @@ contains end do end if - if (hyperelasticity) then - do i = eqn_idx%xi%beg, eqn_idx%xi%end + 1 - @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - @:ACC_SETUP_SFs(q_prim_vf(i)) - end do - end if - if (cont_damage) then @:ALLOCATE(q_prim_vf(eqn_idx%damage)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) @@ -986,18 +977,12 @@ contains end do end if - if (elasticity) then + if (hypoelasticity) then do i = eqn_idx%stress%beg, eqn_idx%stress%end @:DEALLOCATE(q_prim_vf(i)%sf) end do end if - if (hyperelasticity) then - do i = eqn_idx%xi%beg, eqn_idx%xi%end + 1 - @:DEALLOCATE(q_prim_vf(i)%sf) - end do - end if - if (cont_damage) then @:DEALLOCATE(q_prim_vf(eqn_idx%damage)%sf) end if diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index b90c5306ae..94d7212683 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1091,7 +1091,6 @@ def check_igr_simulation(self): acoustic_source = self.get("acoustic_source", "F") == "T" relax = self.get("relax", "F") == "T" mhd = self.get("mhd", "F") == "T" - hyperelasticity = self.get("hyperelasticity", "F") == "T" cyl_coord = self.get("cyl_coord", "F") == "T" probe_wrt = self.get("probe_wrt", "F") == "T" int_comp = self.get("int_comp", 0) @@ -1110,7 +1109,6 @@ def check_igr_simulation(self): self.prohibit(acoustic_source, "IGR does not support acoustic sources") self.prohibit(relax, "IGR does not support phase change") self.prohibit(mhd, "IGR does not support magnetohydrodynamics") - self.prohibit(hyperelasticity, "IGR does not support hyperelasticity") self.prohibit(cyl_coord, "IGR does not support cylindrical or axisymmetric coordinates") self.prohibit(probe_wrt, "IGR does not support probe writes") self.prohibit(int_comp > 0, "IGR does not support int_comp > 0") @@ -1389,21 +1387,6 @@ def check_probe_integral_output(self): self.prohibit(integral_wrt and fd_order is None, "fd_order must be specified for integral_wrt") self.prohibit(integral_wrt and not bubbles_euler, "integral_wrt requires bubbles_euler to be enabled") - def check_hyperelasticity(self): - """Checks hyperelasticity constraints""" - hyperelasticity = self.get("hyperelasticity", "F") == "T" - pre_stress = self.get("pre_stress", "F") == "T" - - self.prohibit(pre_stress and not hyperelasticity, "pre_stress requires hyperelasticity to be enabled") - - if not hyperelasticity: - return - - model_eqns = self.get("model_eqns") - - self.prohibit(model_eqns == 1, "hyperelasticity is not supported for model_eqns = 1") - self.prohibit(model_eqns is not None and model_eqns > 3, "hyperelasticity is not supported for model_eqns > 3") - # Pre-Process Specific Checks def check_restart(self): @@ -2324,7 +2307,6 @@ def validate_common(self): self.check_qbmm_and_polydisperse() self.check_adv_n() self.check_hypoelasticity() - self.check_hyperelasticity() self.check_phase_change() self.check_ibm() self.check_stiffened_eos() diff --git a/toolchain/mfc/gen_case_constraints_docs.py b/toolchain/mfc/gen_case_constraints_docs.py index 621e7f2ddf..2c73038783 100644 --- a/toolchain/mfc/gen_case_constraints_docs.py +++ b/toolchain/mfc/gen_case_constraints_docs.py @@ -310,7 +310,7 @@ def render_markdown(rules: Iterable[Rule]) -> str: # Define major feature groups (excluding IGR) major_features = { - "Physics Models": ["mhd", "surface_tension", "hypoelasticity", "hyperelasticity", "relax", "viscous", "acoustic_source"], + "Physics Models": ["mhd", "surface_tension", "hypoelasticity", "relax", "viscous", "acoustic_source"], "Bubble Models": ["bubbles_euler", "bubbles_lagrange", "qbmm", "polydisperse", "adv_n"], "Numerics": ["riemann_solver", "weno_order", "muscl_order"], "Geometry": ["ib", "cyl_coord"], diff --git a/toolchain/mfc/lint_docs.py b/toolchain/mfc/lint_docs.py index e852284821..64d9ea4d61 100644 --- a/toolchain/mfc/lint_docs.py +++ b/toolchain/mfc/lint_docs.py @@ -451,7 +451,6 @@ def check_physics_docs_coverage(repo_root: Path) -> list[str]: "check_body_forces", "check_continuum_damage", "check_grcbc", - "check_hyperelasticity", "check_interface_compression", "check_ibm", "check_igr_simulation", diff --git a/toolchain/mfc/params/ast_analyzer.py b/toolchain/mfc/params/ast_analyzer.py index ad7d94fd42..5517979f12 100644 --- a/toolchain/mfc/params/ast_analyzer.py +++ b/toolchain/mfc/params/ast_analyzer.py @@ -659,7 +659,6 @@ def classify_message(msg: str) -> str: "alt_soundspeed": {"title": "Alternative Sound Speed", "category": "physics"}, "surface_tension": {"title": "Surface Tension Model", "category": "physics"}, "hypoelasticity": {"title": "Hypoelasticity", "category": "physics"}, - "hyperelasticity": {"title": "Hyperelasticity", "category": "physics"}, "relax": {"title": "Phase Change (Relaxation)", "category": "physics"}, "viscous": {"title": "Viscosity", "category": "physics"}, "acoustic_source": {"title": "Acoustic Sources", "category": "physics"}, diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 8eb3002f44..886d4252db 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -125,7 +125,7 @@ def _fc(name: str, default: int) -> int: "grid": "Grid", "weno": "WENO", "viscosity": "Viscosity", - "elasticity": "Elasticity", + "hypoelasticity": "Hypoelasticity", "surface_tension": "Surface tension", "acoustic": "Acoustic", "ib": "Immersed boundary", @@ -618,9 +618,8 @@ def _load(): # Viscosity _r("viscous", LOG, {"viscosity"}) - # Elasticity - for n in ["hypoelasticity", "hyperelasticity"]: - _r(n, LOG, {"elasticity"}) + # Hypoelasticity + _r("hypoelasticity", LOG, {"hypoelasticity"}) # Surface tension _r("sigma", REAL, {"surface_tension"}, math=r"\f$\sigma\f$") @@ -772,7 +771,6 @@ def _load(): "mixlayer_perturb", "perturb_flow", "perturb_sph", - "pre_stress", "elliptic_smoothing", "simplex_perturb", "alt_soundspeed", @@ -864,9 +862,9 @@ def _load(): for f in range(1, NF + 1): _r(f"{px}alpha({f})", A_REAL, math=r"\f$\alpha_" + str(f) + r"\f$") _r(f"{px}alpha_rho({f})", A_REAL, math=r"\f$\alpha \rho\f$") - # Elasticity stress tensor + # Hypoelastic stress tensor for j in range(1, 7): - _r(f"{px}tau_e({j})", A_REAL, {"elasticity"}, math=r"\f$\tau_e\f$") + _r(f"{px}tau_e({j})", A_REAL, {"hypoelasticity"}, math=r"\f$\tau_e\f$") if i >= 2: for j in range(1, i): _r(f"{px}alter_patch({j})", LOG) @@ -890,7 +888,7 @@ def _load(): px = f"fluid_pp({f})%" for a, sym in [("gamma", r"\f$\gamma_k\f$"), ("pi_inf", r"\f$\pi_{\infty,k}\f$"), ("cv", r"\f$c_{v,k}\f$"), ("qv", r"\f$q_{v,k}\f$"), ("qvp", r"\f$q'_{v,k}\f$")]: _r(f"{px}{a}", REAL, math=sym) - _r(f"{px}G", REAL, {"elasticity"}, math=r"\f$G_k\f$") + _r(f"{px}G", REAL, {"hypoelasticity"}, math=r"\f$G_k\f$") _r(f"{px}Re(1)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (shear)") _r(f"{px}Re(2)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (bulk)") _r(f"{px}non_newtonian", LOG, {"viscosity"}, math=r"\mathrm{non\text{-}Newtonian}_k") @@ -1250,7 +1248,6 @@ def _nv(targets: set, *names: str) -> None: "sigma", "adv_n", "hypoelasticity", - "hyperelasticity", "surface_tension", "relativity", "ib", @@ -1416,7 +1413,6 @@ def _nv(targets: set, *names: str) -> None: "mixlayer_perturb", "mixlayer_perturb_nk", "mixlayer_perturb_k0", - "pre_stress", "elliptic_smoothing", "elliptic_smoothing_iters", "simplex_perturb", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 989e351948..bd09bb4f65 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -101,7 +101,6 @@ "Re_inv": "Inverse Reynolds number", "viscous": "Enable viscous effects", "hypoelasticity": "Enable hypoelastic model", - "hyperelasticity": "Enable hyperelastic model", "surface_tension": "Enable surface tension effects", "chemistry": "Enable chemical reactions", "reactive_burn": "Enable condensed-phase reactive burn (programmed pressure burn on the multi-fluid model)", @@ -186,7 +185,6 @@ "perturb_flow": "Enable flow perturbation", "perturb_sph": "Enable spherical perturbation", "cfl_dt": "Enable CFL-based time stepping", - "pre_stress": "Enable pre-stress initialization", "elliptic_smoothing": "Enable elliptic smoothing", "simplex_perturb": "Enable simplex noise perturbation", "n_start_old": "Starting index from previous simulation", @@ -700,7 +698,7 @@ def get_math_symbol(param_name: str) -> str: "time": "Time stepping and integration", "output": "Output and visualization", "chemistry": "Chemical reactions and species transport", - "elasticity": "Elastic and hyperelastic materials", + "hypoelasticity": "Hypoelastic materials", "acoustic": "Acoustic sources and wave generation", "ib": "Immersed boundary method", "grid": "Computational grid and domain",