Remove non-operational hyperelasticity capability

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                |

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`

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",

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

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

src/common/m_global_parameters_common.fpp

@@ -44,28 +44,25 @@ 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)
     !> @{
     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
     integer, dimension(3, 2) :: shear_BC_flip_indices  !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num)
     !> @}
 
 #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.
 

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

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)

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)