Skip to content

Restore constrained neo-Hookean reference-map hyperelasticity #1677

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
sbryngelson wants to merge 7 commits into
base: master
Choose a base branch
from
Draft

Restore constrained neo-Hookean reference-map hyperelasticity #1677

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e3068ab
remove non-operational hyperelasticity
sbryngelson Jul 25, 2026
3b146a9
Revert "remove non-operational hyperelasticity"
sbryngelson Jul 25, 2026
4aa4488
repair constrained hyperelasticity implementation
sbryngelson Jul 25, 2026
b042bb2
Merge branch 'master' into agent/remove-hyperelasticity
sbryngelson Jul 25, 2026
f9f859e
Merge branch 'agent/remove-hyperelasticity' into agent/restore-hypere…
sbryngelson Jul 25, 2026
8a2f2ef
remove residual elasticity paths
sbryngelson Jul 26, 2026
f93fba3
Merge updated hyperelasticity removal baseline
sbryngelson Jul 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/documentation/case.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) |
Expand Down Expand Up @@ -1125,6 +1124,7 @@ Note: For relativistic flow, the conservative and primitive densities are differ
| ---: | :---: | :--- |
| `hypoelasticity` | Logical | Enable hypoelasticity simulation |
| `hyperelasticity` | Logical | Enable hyperelasticity simulation |
| `hyper_model` | Integer | Hyperelastic constitutive model [1]: neo-Hookean |
| `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 |
Expand Down
4 changes: 2 additions & 2 deletions src/common/m_variables_conversion.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -751,7 +751,7 @@ contains

if (hyperelasticity) then
$:GPU_LOOP(parallelism='[seq]')
do i = eqn_idx%xi%beg, eqn_idx%xi%end
do i = eqn_idx%xi%beg, eqn_idx%xi%end + 1
qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/rho_K
end do
end if
Expand Down Expand Up @@ -1004,7 +1004,7 @@ contains
! 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
do i = eqn_idx%xi%beg, eqn_idx%xi%end + 1
q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l)
end do
end if
Expand Down
22 changes: 7 additions & 15 deletions src/pre_process/m_assign_variables.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +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
Expand Down Expand Up @@ -384,25 +383,18 @@ contains
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
xi_cart(1) = x_cc(j)
xi_cart(2) = y_cc(k)
xi_cart(3) = z_cc(l)

! 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
do i = eqn_idx%stress%beg, eqn_idx%stress%end
q_prim_vf(i)%sf(j, k, l) = 0._wp
end do
q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l) = 0._wp
end if

if (mpp_lim .and. bubbles_euler) then
Expand Down
2 changes: 0 additions & 2 deletions src/pre_process/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
3 changes: 0 additions & 3 deletions src/simulation/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,6 @@ 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
logical :: shear_stress !< Shear stresses
logical :: bulk_stress !< Bulk stresses
Expand All @@ -92,7 +91,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
Expand Down Expand Up @@ -391,7 +389,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.
Expand Down
58 changes: 17 additions & 41 deletions src/simulation/m_hyperelastic.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,14 @@ contains
rho = max(rho, sgm_eps)
G_local = max(G_local, sgm_eps)

$:GPU_LOOP(parallelism='[seq]')
do i = 1, b_size - 1
q_prim_vf(eqn_idx%stress%beg + i - 1)%sf(j, k, l) = 0._wp
q_cons_vf(eqn_idx%stress%beg + i - 1)%sf(j, k, l) = 0._wp
end do
q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l) = 0._wp
q_cons_vf(eqn_idx%xi%end + 1)%sf(j, k, l) = 0._wp

if (G_local > verysmall) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, tensor_size
Expand Down Expand Up @@ -153,25 +161,21 @@ contains
! 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 3: tensora stores F^T, so b = F F^T = transpose(tensora)*tensora
tensorb(1) = tensora(1)**2 + tensora(4)**2 + tensora(7)**2
tensorb(5) = tensora(2)**2 + tensora(5)**2 + tensora(8)**2
tensorb(9) = tensora(3)**2 + tensora(6)**2 + tensora(9)**2
tensorb(2) = tensora(1)*tensora(2) + tensora(4)*tensora(5) + tensora(7)*tensora(8)
tensorb(3) = tensora(1)*tensora(3) + tensora(4)*tensora(6) + tensora(7)*tensora(9)
tensorb(6) = tensora(2)*tensora(3) + tensora(5)*tensora(6) + tensora(8)*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
call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G_local, j, k, l)
! 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
Expand All @@ -181,6 +185,7 @@ contains
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
q_cons_vf(eqn_idx%xi%end + 1)%sf(j, k, l) = rho*q_prim_vf(eqn_idx%xi%end + 1)%sf(j, k, l)
end if
end if
end do
Expand Down Expand Up @@ -218,35 +223,6 @@ contains

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()

Expand Down
Loading
Loading