@@ -14,8 +14,9 @@ module m_hypoelastic
1414
1515 implicit none
1616
17- private; public :: s_initialize_hypoelastic_module, s_finalize_hypoelastic_module, s_compute_hypoelastic_rhs_legacy, &
18- & s_compute_hypoelastic_rhs_iface, s_compute_hypoelastic_rhs_axisym_geom_iface, s_compute_damage_state
17+ private; public :: s_initialize_hypoelastic_module, s_finalize_hypoelastic_module, &
18+ & s_compute_hypoelastic_rhs_finite_diff_per_sweep, s_compute_hypoelastic_rhs_iface, &
19+ & s_compute_hypoelastic_rhs_axisym_geom_iface, s_compute_damage_state
1920
2021 real (wp), allocatable, dimension (:) :: Gs_hypo
2122 $:GPU_DECLARE(create= ' [Gs_hypo]'
@@ -84,7 +85,7 @@ contains
8485 !! @param idir Dimension splitting index
8586 !! @param q_prim_vf Primitive variables
8687 !! @param rhs_vf rhs variables
87- subroutine s_compute_hypoelastic_rhs_legacy (idir , q_prim_vf , rhs_vf )
88+ subroutine s_compute_hypoelastic_rhs_finite_diff_per_sweep (idir , q_prim_vf , rhs_vf )
8889
8990 integer , intent (in ) :: idir
9091 type(scalar_field), dimension (sys_size), intent (in ) :: q_prim_vf
@@ -348,7 +349,7 @@ contains
348349 $:END_GPU_PARALLEL_LOOP()
349350 end if
350351
351- end subroutine s_compute_hypoelastic_rhs_legacy
352+ end subroutine s_compute_hypoelastic_rhs_finite_diff_per_sweep
352353
353354 !> Interface- consistent hypoelastic RHS (Mode 2 : HLL/ HLLC). Uses interface velocities from the Riemann solver to compute
354355 !! velocity gradients. Called once after all dimensional sweeps. Supports 1D , 2D Cartesian, 2D axisymmetric, and 3D Cartesian.
@@ -361,9 +362,7 @@ contains
361362 type(scalar_field), dimension (sys_size), intent (inout ) :: rhs_vf
362363 type(vector_field), dimension (:), intent (in ) :: nc_iface_vel_n
363364 real (wp) :: rho_K, G_K
364- real (wp) :: A_x, B_x, C_x, D_x, E_x, F_x, H_x, J1_x, J2_x
365- real (wp) :: A_y, B_y, C_y, D_y, E_y, F_y, H_y, J1_y, J2_y
366- real (wp) :: A_z, B_z, C_z, D_z, E_z, F_z, H_z, J1_z, J2_z
365+ real (wp) :: trace, shear, shear2, diag, diag_z, offdiag, cross1, cross2
367366 real (wp) :: txx, txy, tyy, txz, tyz, tzz
368367 integer :: i, k, l, q, r
369368 integer :: ndirs
@@ -477,8 +476,8 @@ contains
477476 end if
478477
479478 if (ndirs == 3 .and. .not. cyl_coord) then
480- $:GPU_PARALLEL_LOOP(collapse= 3 ,private= ' [txx, txy, tyy, txz, tyz, tzz, A_x, B_x, C_y, D_y, C_z, D_z, B_y, H_z, J1_z , &
481- & J2_z, H_y, J1_y, J2_y, B_z, C_x, D_x, A_y, E_z, F_z, E_x, F_x, E_y, F_y, A_z, H_x, J1_x, J2_x ]' )
479+ $:GPU_PARALLEL_LOOP(collapse= 3 ,private= ' [txx, txy, tyy, txz, tyz, tzz, trace, shear, shear2, diag, diag_z, offdiag , &
480+ & cross1, cross2 ]' )
482481 do q = 0 , p
483482 do l = 0 , n
484483 do k = 0 , m
@@ -490,52 +489,51 @@ contains
490489 tzz = q_prim_vf(eqn_idx%stress%beg + 5 )%sf(k, l, q)
491490
492491 ! z- direction contributions to tau_xx
493- C_z = - (2._wp / 3._wp * G_K_field(k, l, q) + txx)
494- D_z = 2._wp * txz
492+ trace = - (2._wp / 3._wp * G_K_field(k, l, q) + txx)
493+ shear = 2._wp * txz
495494 rhs_vf(eqn_idx%stress%beg)%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg)%sf(k, l, q) + rho_K_field(k, l, &
496- & q)* (C_z * dw_dz_hypo(k, l, q) + D_z * du_dz_hypo(k, l, q))
495+ & q)* (trace * dw_dz_hypo(k, l, q) + shear * du_dz_hypo(k, l, q))
497496
498497 ! z- direction contributions to tau_xy
499- H_z = - txy
500- J1_z = tyz
501- J2_z = txz
498+ offdiag = - txy
499+ cross1 = tyz
500+ cross2 = txz
502501 rhs_vf(eqn_idx%stress%beg + 1 )%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 1 )%sf(k, l, q) + rho_K_field(k, &
503- & l, q)* (H_z * dw_dz_hypo(k, l, q) + J1_z * du_dz_hypo(k, l, q) + J2_z * dv_dz_hypo(k, l, q))
502+ & l, q)* (offdiag * dw_dz_hypo(k, l, q) + cross1 * du_dz_hypo(k, l, q) + cross2 * dv_dz_hypo(k, l, q))
504503
505- ! tau_yy: z- direction contributions
506- E_z = - (2._wp / 3._wp * G_K_field(k, l, q) + tyy)
507- F_z = 2._wp * tyz
504+ ! z- direction contributions to tau_yy
505+ trace = - (2._wp / 3._wp * G_K_field(k, l, q) + tyy)
506+ shear = 2._wp * tyz
508507 rhs_vf(eqn_idx%stress%beg + 2 )%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 2 )%sf(k, l, q) + rho_K_field(k, &
509- & l, q)* (E_z * dw_dz_hypo(k, l, q) + F_z * dv_dz_hypo(k, l, q))
508+ & l, q)* (trace * dw_dz_hypo(k, l, q) + shear * dv_dz_hypo(k, l, q))
510509
511510 ! tau_xz (stress%beg+3 )
512- B_x = G_K_field(k, l, q) + txx
513- H_y = - txz
514- J1_y = tyz
515- J2_y = txy
516- B_z = G_K_field(k, l, q) + tzz
511+ diag = G_K_field(k, l, q) + txx
512+ offdiag = - txz
513+ cross1 = tyz
514+ cross2 = txy
515+ diag_z = G_K_field(k, l, q) + tzz
517516 rhs_vf(eqn_idx%stress%beg + 3 )%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 3 )%sf(k, l, q) + rho_K_field(k, &
518- & l, q)* (B_x * dw_dx_hypo(k, l, q) + H_y * dv_dy_hypo(k, l, q) + J1_y * du_dy_hypo(k, l, &
519- & q) + J2_y * dw_dy_hypo(k, l, q) + B_z * du_dz_hypo(k, l, q))
517+ & l, q)* (diag * dw_dx_hypo(k, l, q) + offdiag * dv_dy_hypo(k, l, q) + cross1 * du_dy_hypo(k, l, &
518+ & q) + cross2 * dw_dy_hypo(k, l, q) + diag_z * du_dz_hypo(k, l, q))
520519
521520 ! tau_yz (stress%beg+4 )
522- H_x = - tyz
523- J1_x = txz
524- J2_x = txy
525- B_y = G_K_field(k, l, q) + tyy
521+ offdiag = - tyz
522+ cross1 = txz
523+ cross2 = txy
524+ diag = G_K_field(k, l, q) + tyy
526525 rhs_vf(eqn_idx%stress%beg + 4 )%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 4 )%sf(k, l, q) + rho_K_field(k, &
527- & l, q)* (H_x * du_dx_hypo(k, l, q) + J1_x * dv_dx_hypo(k, l, q) + J2_x * dw_dx_hypo(k, l, &
528- & q) + B_y * dw_dy_hypo(k, l, q) + B_z * dv_dz_hypo(k, l, q))
526+ & l, q)* (offdiag * du_dx_hypo(k, l, q) + cross1 * dv_dx_hypo(k, l, q) + cross2 * dw_dx_hypo(k, l, &
527+ & q) + diag * dw_dy_hypo(k, l, q) + diag_z * dv_dz_hypo(k, l, q))
529528
530529 ! tau_zz (stress%beg+5 )
531- E_x = - (2._wp / 3._wp * G_K_field(k, l, q) + tzz)
532- F_x = 2._wp * txz
533- E_y = - (2._wp / 3._wp * G_K_field(k, l, q) + tzz)
534- F_y = 2._wp * tyz
535- A_z = 4._wp / 3._wp * G_K_field(k, l, q) + tzz
530+ trace = - (2._wp / 3._wp * G_K_field(k, l, q) + tzz)
531+ shear = 2._wp * txz
532+ shear2 = 2._wp * tyz
533+ diag = 4._wp / 3._wp * G_K_field(k, l, q) + tzz
536534 rhs_vf(eqn_idx%stress%beg + 5 )%sf(k, l, q) = rhs_vf(eqn_idx%stress%beg + 5 )%sf(k, l, q) + rho_K_field(k, &
537- & l, q)* (E_x * du_dx_hypo(k, l, q) + F_x * dw_dx_hypo(k, l, q) + E_y * dv_dy_hypo(k, l, &
538- & q) + F_y * dw_dy_hypo(k, l, q) + A_z * dw_dz_hypo(k, l, q))
535+ & l, q)* (trace * du_dx_hypo(k, l, q) + shear * dw_dx_hypo(k, l, q) + trace * dv_dy_hypo(k, l, &
536+ & q) + shear2 * dw_dy_hypo(k, l, q) + diag * dw_dz_hypo(k, l, q))
539537 end do
540538 end do
541539 end do
@@ -548,6 +546,14 @@ contains
548546
549547 end subroutine s_compute_hypoelastic_rhs_iface
550548
549+ !> Axisymmetric geometric source terms for the hypoelastic stress evolution, using interface velocities. Adds the v/ r and div(u)
550+ !! contributions that arise in cylindrical (r- z) coordinates: tau_xx, tau_xr, tau_rr get a - rho* (v/ r) source; tau_thetatheta
551+ !! gets a combined divergence and hoop- stress source. Called from s_compute_hypoelastic_rhs_iface when grid_geometry == 2 .
552+ !! @param q_prim_vf Primitive variables
553+ !! @param rhs_vf rhs variables
554+ !! @param nc_iface_vel_x_vf Interface velocities in x- direction
555+ !! @param nc_iface_vel_y_vf Interface velocities in y- direction
556+ !! @param weight Sub- step weighting factor
551557 subroutine s_compute_hypoelastic_rhs_axisym_geom_iface (q_prim_vf , rhs_vf , nc_iface_vel_x_vf , nc_iface_vel_y_vf , weight )
552558
553559 type(scalar_field), dimension (sys_size), intent (in ) :: q_prim_vf
0 commit comments