From 7922dcc1b9f5ed2fb212e8f1d9851170e0dc7024 Mon Sep 17 00:00:00 2001 From: mcarcana Date: Thu, 23 Jan 2025 14:20:07 -0500 Subject: [PATCH 1/9] basic 1D/2D hyperelasticity scheme added --- src/simulation/m_hyperelastic.fpp | 254 +++++++++++++++++++++++++++-- src/simulation/m_rhs.fpp | 2 +- src/simulation/m_time_steppers.fpp | 2 +- 3 files changed, 244 insertions(+), 14 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index d01440257..e926d5b1b 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -94,8 +94,9 @@ contains !! calculate the inverse of grad_xi to obtain F, F is a nxn tensor !! calculate the FFtranspose to obtain the btensor, btensor is nxn tensor !! btensor is symmetric, save the data space - subroutine s_hyperelastic_rmt_stress_update(q_cons_vf, q_prim_vf) - + subroutine s_hyperelastic_rmt_stress_update(num_dims, q_cons_vf, q_prim_vf) + + integer, intent(in) :: num_dims type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf @@ -106,6 +107,170 @@ contains real(wp) :: G integer :: j, k, l, i, r + if (num_dims == 1) then + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & + !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) + do l = 0, p + do k = 0, n + do j = 0, m + !$acc loop seq + do i = 1, num_fluids + alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) + alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) + end do + ! 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, j, k, l, G, Gs) + rho = max(rho, sgm_eps) + G = max(G, sgm_eps) + !if ( G <= verysmall ) G_K = 0d0 + + if ( G > verysmall ) then + !$acc loop seq + do i = 1, tensor_size + tensora(i) = 0d0 + end do + ! STEP 1: computing the grad_xi tensor using finite differences + ! grad_xi definition / organization + ! number for the tensor 1-2: dxix_dxy + !$acc loop seq + do r = -fd_number, fd_number + ! derivatives in the x-direction + tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) + end do + ! STEP 2a: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1) + ! STEP 2b: computing the inverse of the grad_xi tensor + tensorb(1) = 1d0/(tensora(1)**2) + + if (tensorb(tensor_size) > verysmall) then + ! STEP 2c: computing the inverse of grad_xi tensor = F + ! tensorb is the adjoint, tensora becomes F + !$acc loop seq + do i = 1, tensor_size - 1 + tensora(i) = tensorb(i)/tensorb(tensor_size) + end do + ! STEP 3: update the btensor, this is consistent with Riemann solvers + ! \b_xx + btensor%vf(1)%sf(j, k, l) = tensorb(1) + ! 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_1D(btensor%vf, q_prim_vf, G, j, k, l) + !elseif (hyper_model == 2) then + ! call s_Mooney_Rivlin_cauchy_solver_1D(btensor%vf, q_prim_vf, G, j, k, l) + !end if + ! STEP 5b: updating the pressure field + q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & + G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma + ! STEP 5c: updating the Cauchy stress conservative scalar field + !$acc loop seq + do i = 1, b_size - 1 + q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) + end do + end if + end if + end do + end do + end do + !$acc end parallel loop + + elseif (num_dims == 2) then + + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & + !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) + do l = 0, p + do k = 0, n + do j = 0, m + !$acc loop seq + do i = 1, num_fluids + alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) + alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) + end do + ! 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, j, k, l, G, Gs) + rho = max(rho, sgm_eps) + G = max(G, sgm_eps) + !if ( G <= verysmall ) G_K = 0d0 + + if ( G > verysmall ) then + !$acc loop seq + do i = 1, tensor_size + tensora(i) = 0d0 + end do + ! STEP 1: computing the grad_xi tensor using finite differences + ! grad_xi definition / organization + ! number for the tensor 1-2: dxix_dx, dxiy_dx + ! 3-4: dxix_dy, dxiy_dy + !$acc loop seq + do r = -fd_number, fd_number + ! derivatives in the x-direction + tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + end do + !print *, 'j :: ',j,', tensor1 :: ',tensora(1),', tensor2 :: ',tensora(2),', tensora(3) :: ',tensora(3),', tensora(4) :: ',tensora(4) + ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse + tensorb(1) = tensora(4) + tensorb(2) = -tensora(3) + tensorb(3) = -tensora(2) + tensorb(4) = tensora(1) + ! STEP 2b: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3) + + if (tensorb(tensor_size) > verysmall) then + ! STEP 2c: computing the inverse of grad_xi tensor = F + ! tensorb is the adjoint, tensora becomes F + !$acc loop 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) = 1d0/tensorb(tensor_size) + ! STEP 3: override adjoint (tensorb) to be F transpose F + tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensorb(4) = tensora(3)**2 + tensora(4)**2 + tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4) + tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2) + + ! STEP 4: update the btensor, this is consistent with Riemann solvers + #:for BIJ, TXY in [(1,1),(2,2),(3,4)] + 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 + !print *,'hyper model :: ',hyper_model + !if (hyper_model == 1) then + call s_neoHookean_cauchy_solver_2D(btensor%vf, q_prim_vf, G, j, k, l) + !elseif (hyper_model == 2) then + ! call s_Mooney_Rivlin_cauchy_solver_2D(btensor%vf, q_prim_vf, G, j, k, l) + !end if + ! print *, 'before j :: ',j,', k :: ',k,', p :: ',q_prim_vf(E_idx)%sf(j,k,l),', val :: ',q_prim_vf(xiend + 1)%sf(j, k, l) + + ! STEP 5b: updating the pressure field + q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & + G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma + + ! print *, 'after j :: ',j,', k :: ',k,', p :: ',q_prim_vf(E_idx)%sf(j,k,l),' val :: ',q_prim_vf(xiend + 1)%sf(j, k, l) + + ! STEP 5c: updating the Cauchy stress conservative scalar field + !$acc loop seq + do i = 1, b_size - 1 + q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) + end do + end if + end if + end do + end do + end do + !$acc end parallel loop + elseif (num_dims == 3) then + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) do l = 0, p - 2 @@ -189,11 +354,11 @@ contains ! 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, j, k, l) - elseif (hyper_model == 2) then - call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) - end if + !if (hyper_model == 1) then + call s_neoHookean_cauchy_solver_3D(btensor%vf, q_prim_vf, G, j, k, l) + !elseif (hyper_model == 2) then + ! call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + !end if ! STEP 5b: updating the pressure field q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma @@ -209,17 +374,81 @@ contains end do end do !$acc end parallel loop + end if end subroutine s_hyperelastic_rmt_stress_update - !> The following subroutine handles the calculation of the btensor. - !! The calculation of the btensor takes qprimvf. + !> The following subroutine handles the calculation of the btensor + !! with a neo-Hookean material model. + !! The calculation of the btensor takes qprimvf. !! @param q_prim_vf Primitive variables !! @param btensor is the output !! calculate the grad_xi, grad_xi is a nxn tensor !! calculate the inverse of grad_xi to obtain F, F is a nxn tensor !! calculate the FFtranspose to obtain the btensor, btensor is nxn tensor !! btensor is symmetric, save the data space - subroutine s_neoHookean_cauchy_solver(btensor, q_prim_vf, G, j, k, l) + subroutine s_neoHookean_cauchy_solver_1D(btensor, q_prim_vf, G, j, k, l) + !$acc routine seq + type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf + type(scalar_field), dimension(b_size), intent(inout) :: btensor + real(wp), intent(in) :: G + integer, intent(in) :: j, k, l + + real(wp) :: trace + real(wp), parameter :: f13 = 1._wp/3._wp + integer :: i + + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace + + ! dividing by the jacobian for neo-Hookean model + ! setting the tensor to the stresses for riemann solver + !$acc loop seq + do i = 1, b_size - 1 + q_prim_vf(strxb + i - 1)%sf(j, k, l) = & + G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + end do + ! compute the invariant without the elastic modulus + q_prim_vf(xiend + 1)%sf(j, k, l) = & + 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) + + end subroutine s_neoHookean_cauchy_solver_1D + + subroutine s_neoHookean_cauchy_solver_2D(btensor, q_prim_vf, G, j, k, l) + !$acc routine seq + type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf + type(scalar_field), dimension(b_size), intent(inout) :: btensor + real(wp), intent(in) :: G + integer, intent(in) :: j, k, l + + real(wp) :: trace + real(wp), parameter :: f13 = 1._wp/3._wp + integer :: i + + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace + btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace + + ! dividing by the jacobian for neo-Hookean model + ! setting the tensor to the stresses for riemann solver + !$acc loop seq + do i = 1, b_size - 1 + q_prim_vf(strxb + i - 1)%sf(j, k, l) = & + G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + ! print *,'j :: ',j,', val :: ', G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + end do + ! compute the invariant without the elastic modulus + q_prim_vf(xiend + 1)%sf(j, k, l) = & + 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) + + end subroutine s_neoHookean_cauchy_solver_2D + + subroutine s_neoHookean_cauchy_solver_3D(btensor, q_prim_vf, G, j, k, l) !$acc routine seq type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf type(scalar_field), dimension(b_size), intent(inout) :: btensor @@ -249,9 +478,10 @@ contains q_prim_vf(xiend + 1)%sf(j, k, l) = & 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) - end subroutine s_neoHookean_cauchy_solver + end subroutine s_neoHookean_cauchy_solver_3D - !> The following subroutine handles the calculation of the btensor. + !> The following subroutine handles the calculation of the btensor + !! with a Mooney-Rivlin material model. !! The calculation of the btensor takes qprimvf. !! @param q_prim_vf Primitive variables !! @param btensor is the output diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index b00f4a0b0..bce7956e3 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -667,7 +667,7 @@ contains call nvtxStartRange("RHS-ELASTIC") if (hyperelasticity) then - call s_hyperelastic_rmt_stress_update(q_cons_qp%vf, q_prim_qp%vf) + call s_hyperelastic_rmt_stress_update(num_dims, q_cons_qp%vf, q_prim_qp%vf) call s_populate_variables_buffers(q_prim_qp%vf, pb, mv) end if call nvtxEndRange diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 7206d8f50..19e4ad653 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -846,7 +846,7 @@ contains if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) call nvtxStartRange("RHS-ELASTIC") - if (hyperelasticity) call s_hyperelastic_rmt_stress_update(q_cons_ts(1)%vf, q_prim_vf) + if (hyperelasticity) call s_hyperelastic_rmt_stress_update(num_dims, q_cons_ts(1)%vf, q_prim_vf) call nvtxEndRange if (ib) then From 08347163ca94c1a90f4f04bd0123da6ca577252e Mon Sep 17 00:00:00 2001 From: mcarcana Date: Thu, 23 Jan 2025 16:57:43 -0500 Subject: [PATCH 2/9] removed print statements in 2d hyper --- src/simulation/m_hyperelastic.fpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index e926d5b1b..5b6ece37c 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -213,7 +213,6 @@ contains tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) end do - !print *, 'j :: ',j,', tensor1 :: ',tensora(1),', tensor2 :: ',tensora(2),', tensora(3) :: ',tensora(3),', tensora(4) :: ',tensora(4) ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse tensorb(1) = tensora(4) tensorb(2) = -tensora(3) @@ -250,14 +249,11 @@ contains !elseif (hyper_model == 2) then ! call s_Mooney_Rivlin_cauchy_solver_2D(btensor%vf, q_prim_vf, G, j, k, l) !end if - ! print *, 'before j :: ',j,', k :: ',k,', p :: ',q_prim_vf(E_idx)%sf(j,k,l),', val :: ',q_prim_vf(xiend + 1)%sf(j, k, l) ! STEP 5b: updating the pressure field q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma - ! print *, 'after j :: ',j,', k :: ',k,', p :: ',q_prim_vf(E_idx)%sf(j,k,l),' val :: ',q_prim_vf(xiend + 1)%sf(j, k, l) - ! STEP 5c: updating the Cauchy stress conservative scalar field !$acc loop seq do i = 1, b_size - 1 From 32141de16b799e19cb270bd57350a557da94707e Mon Sep 17 00:00:00 2001 From: mcarcana Date: Thu, 23 Jan 2025 17:17:07 -0500 Subject: [PATCH 3/9] removed print statements from m_hyperelastic part 2 --- src/simulation/m_hyperelastic.fpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 5b6ece37c..9527e1cae 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -243,7 +243,6 @@ contains ! 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 - !print *,'hyper model :: ',hyper_model !if (hyper_model == 1) then call s_neoHookean_cauchy_solver_2D(btensor%vf, q_prim_vf, G, j, k, l) !elseif (hyper_model == 2) then @@ -436,7 +435,6 @@ contains do i = 1, b_size - 1 q_prim_vf(strxb + i - 1)%sf(j, k, l) = & G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) - ! print *,'j :: ',j,', val :: ', G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) end do ! compute the invariant without the elastic modulus q_prim_vf(xiend + 1)%sf(j, k, l) = & From 52ba7d5caab860e45f329f681031cab7db8c8a3c Mon Sep 17 00:00:00 2001 From: mcarcana Date: Thu, 23 Jan 2025 17:20:34 -0500 Subject: [PATCH 4/9] corrected double precision intrinsics in m_hyper --- src/simulation/m_hyperelastic.fpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 9527e1cae..9bd22f9f5 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -123,12 +123,12 @@ contains alpha_rho_k, Re, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0d0 + !if ( G <= verysmall ) G_K = 0_wp if ( G > verysmall ) then !$acc loop seq do i = 1, tensor_size - tensora(i) = 0d0 + tensora(i) = 0_wp end do ! STEP 1: computing the grad_xi tensor using finite differences ! grad_xi definition / organization @@ -193,12 +193,12 @@ contains alpha_rho_k, Re, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0d0 + !if ( G <= verysmall ) G_K = 0_wp if ( G > verysmall ) then !$acc loop seq do i = 1, tensor_size - tensora(i) = 0d0 + tensora(i) = 0_wp end do ! STEP 1: computing the grad_xi tensor using finite differences ! grad_xi definition / organization From aa9a0072c69ce5be8af2c7469c34e7369d92ad3b Mon Sep 17 00:00:00 2001 From: mcarcana Date: Thu, 23 Jan 2025 17:26:23 -0500 Subject: [PATCH 5/9] corrected double precision intrinsics in m_hyper part 2 --- src/simulation/m_hyperelastic.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 9bd22f9f5..8f6c2de06 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -141,7 +141,7 @@ contains ! STEP 2a: computing the determinant of the grad_xi tensor tensorb(tensor_size) = tensora(1) ! STEP 2b: computing the inverse of the grad_xi tensor - tensorb(1) = 1d0/(tensora(1)**2) + tensorb(1) = 1_wp/(tensora(1)**2) if (tensorb(tensor_size) > verysmall) then ! STEP 2c: computing the inverse of grad_xi tensor = F @@ -229,7 +229,7 @@ contains tensora(i) = tensorb(i)/tensorb(tensor_size) end do ! STEP 2d: computing the J = det(F) = 1/det(\grad{\xi}) - tensorb(tensor_size) = 1d0/tensorb(tensor_size) + tensorb(tensor_size) = 1_wp/tensorb(tensor_size) ! STEP 3: override adjoint (tensorb) to be F transpose F tensorb(1) = tensora(1)**2 + tensora(2)**2 tensorb(4) = tensora(3)**2 + tensora(4)**2 From 20d84ece6405d5986a8e655450d2172bb43661e1 Mon Sep 17 00:00:00 2001 From: mcarcana Date: Thu, 23 Jan 2025 17:34:33 -0500 Subject: [PATCH 6/9] actually corrected double precision intrinsics in m_hyper part 3 --- src/simulation/m_hyperelastic.fpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 8f6c2de06..36e651a03 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -123,12 +123,12 @@ contains alpha_rho_k, Re, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0_wp + !if ( G <= verysmall ) G_K = 0._wp if ( G > verysmall ) then !$acc loop seq do i = 1, tensor_size - tensora(i) = 0_wp + tensora(i) = 0._wp end do ! STEP 1: computing the grad_xi tensor using finite differences ! grad_xi definition / organization @@ -141,7 +141,7 @@ contains ! STEP 2a: computing the determinant of the grad_xi tensor tensorb(tensor_size) = tensora(1) ! STEP 2b: computing the inverse of the grad_xi tensor - tensorb(1) = 1_wp/(tensora(1)**2) + tensorb(1) = 1._wp/(tensora(1)**2) if (tensorb(tensor_size) > verysmall) then ! STEP 2c: computing the inverse of grad_xi tensor = F @@ -193,12 +193,12 @@ contains alpha_rho_k, Re, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0_wp + !if ( G <= verysmall ) G_K = 0._wp if ( G > verysmall ) then !$acc loop seq do i = 1, tensor_size - tensora(i) = 0_wp + tensora(i) = 0._wp end do ! STEP 1: computing the grad_xi tensor using finite differences ! grad_xi definition / organization @@ -229,7 +229,7 @@ contains 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) + tensorb(tensor_size) = 1._wp/tensorb(tensor_size) ! STEP 3: override adjoint (tensorb) to be F transpose F tensorb(1) = tensora(1)**2 + tensora(2)**2 tensorb(4) = tensora(3)**2 + tensora(4)**2 @@ -281,12 +281,12 @@ contains alpha_rho_k, Re, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0_wp + !if ( G <= verysmall ) G_K = 0._wp if (G > verysmall) then !$acc loop seq do i = 1, tensor_size - tensora(i) = 0_wp + tensora(i) = 0._wp end do ! STEP 1: computing the grad_xi tensor using finite differences ! grad_xi definition / organization @@ -333,7 +333,7 @@ contains end do ! STEP 2d: computing the J = det(F) = 1/det(\grad{\xi}) - tensorb(tensor_size) = 1_wp/tensorb(tensor_size) + 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 From 9bcf5e25493c420986ceb55b215891eb4536f4aa Mon Sep 17 00:00:00 2001 From: mcarcana Date: Wed, 29 Jan 2025 11:41:25 -0500 Subject: [PATCH 7/9] resolving conflicts in m_hyperelastic --- src/simulation/m_hyperelastic.fpp | 478 ++++++++++------------------- src/simulation/m_rhs.fpp | 2 +- src/simulation/m_time_steppers.fpp | 2 +- 3 files changed, 158 insertions(+), 324 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 2611f7fa9..ba304004f 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -94,9 +94,8 @@ contains !! calculate the inverse of grad_xi to obtain F, F is a nxn tensor !! calculate the FFtranspose to obtain the btensor, btensor is nxn tensor !! btensor is symmetric, save the data space - subroutine s_hyperelastic_rmt_stress_update(num_dims, q_cons_vf, q_prim_vf) - - integer, intent(in) :: num_dims + 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 @@ -107,271 +106,155 @@ contains real(wp) :: G integer :: j, k, l, i, r - if (num_dims == 1) then - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & - !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) - do l = 0, p - do k = 0, n - do j = 0, m - !$acc loop seq - do i = 1, num_fluids - alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) - alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - ! 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, j, k, l, G, Gs) - rho = max(rho, sgm_eps) - G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0._wp - - if ( G > verysmall ) then - !$acc loop 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-2: dxix_dxy - !$acc loop seq - do r = -fd_number, fd_number - ! derivatives in the x-direction - tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) - end do - ! STEP 2a: computing the determinant of the grad_xi tensor - tensorb(tensor_size) = tensora(1) - ! STEP 2b: computing the inverse of the grad_xi tensor - tensorb(1) = 1._wp/(tensora(1)**2) - - if (tensorb(tensor_size) > verysmall) then - ! STEP 2c: computing the inverse of grad_xi tensor = F - ! tensorb is the adjoint, tensora becomes F - !$acc loop seq - do i = 1, tensor_size - 1 - tensora(i) = tensorb(i)/tensorb(tensor_size) - end do - ! STEP 3: update the btensor, this is consistent with Riemann solvers - ! \b_xx - btensor%vf(1)%sf(j, k, l) = tensorb(1) - ! 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_1D(btensor%vf, q_prim_vf, G, j, k, l) - !elseif (hyper_model == 2) then - ! call s_Mooney_Rivlin_cauchy_solver_1D(btensor%vf, q_prim_vf, G, j, k, l) - !end if - ! STEP 5b: updating the pressure field - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma - ! STEP 5c: updating the Cauchy stress conservative scalar field - !$acc loop seq - do i = 1, b_size - 1 - q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) - end do - end if - end if - end do - end do - end do - !$acc end parallel loop - - elseif (num_dims == 2) then - - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & - !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) - do l = 0, p - do k = 0, n - do j = 0, m - !$acc loop seq - do i = 1, num_fluids - alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) - alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - ! 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, j, k, l, G, Gs) + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & + !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) + do l = 0, p + do k = 0, n + do j = 0, m + !$acc loop seq + do i = 1, num_fluids + alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) + alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) + end do + + ! 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, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) !if ( G <= verysmall ) G_K = 0._wp - if ( G > verysmall ) then - !$acc loop 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-2: dxix_dx, dxiy_dx - ! 3-4: dxix_dy, dxiy_dy - !$acc loop seq - do r = -fd_number, fd_number - ! derivatives in the x-direction - tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) - tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) - ! derivatives in the y-direction - tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) - end do - ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse - tensorb(1) = tensora(4) - tensorb(2) = -tensora(3) - tensorb(3) = -tensora(2) - tensorb(4) = tensora(1) - ! STEP 2b: computing the determinant of the grad_xi tensor - tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3) - - if (tensorb(tensor_size) > verysmall) then - ! STEP 2c: computing the inverse of grad_xi tensor = F - ! tensorb is the adjoint, tensora becomes F - !$acc loop 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: override adjoint (tensorb) to be F transpose F - tensorb(1) = tensora(1)**2 + tensora(2)**2 - tensorb(4) = tensora(3)**2 + tensora(4)**2 - tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4) - tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2) - - ! STEP 4: update the btensor, this is consistent with Riemann solvers - #:for BIJ, TXY in [(1,1),(2,2),(3,4)] - 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_2D(btensor%vf, q_prim_vf, G, j, k, l) - !elseif (hyper_model == 2) then - ! call s_Mooney_Rivlin_cauchy_solver_2D(btensor%vf, q_prim_vf, G, j, k, l) - !end if - - ! STEP 5b: updating the pressure field - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma - - ! STEP 5c: updating the Cauchy stress conservative scalar field - !$acc loop seq - do i = 1, b_size - 1 - q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) - end do - end if - end if - end do - end do - end do - !$acc end parallel loop - elseif (num_dims == 3) then - - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & - !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) - do l = 0, p - 2 - do k = 0, n - 2 - do j = 2, m - 2 - !$acc loop seq - do i = 1, num_fluids - alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) - alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - ! 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, j, k, l, G, Gs) - rho = max(rho, sgm_eps) - G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0._wp - - if (G > verysmall) then - !$acc loop 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 - !$acc loop seq - do r = -fd_number, fd_number - ! derivatives in the x-direction - tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) - tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) - tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j) - ! derivatives in the y-direction - tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k) - ! derivatives in the z-direction - tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l) - tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l) - tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(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 + if (G > verysmall) then !$acc loop seq - do i = 1, tensor_size - 1 - tensora(i) = tensorb(i)/tensorb(tensor_size) + do i = 1, tensor_size + tensora(i) = 0._wp 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, j, k, l) - elseif (hyper_model == 2) then - call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) - end if - ! Step 5b: updating the pressure field - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma - ! Step 5c: updating the Cauchy stress conservative scalar field + ! STEP 1: computing the grad_xi tensor using finite differences + ! grad_xi definition / organization + ! tensora(1,2,3): dxix_dx, dxiy_dx, dxiz_dx + ! tensora(4,5,6): dxix_dy, dxiy_dy, dxiz_dy + ! tensora(7,8,9): dxix_dz, dxiy_dz, dxiz_dz !$acc loop seq - do i = 1, b_size - 1 - q_cons_vf(strxb + i - 1)%sf(j, k, l) = & - rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) + do r = -fd_number, fd_number + ! derivatives in the x-direction + tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) + if (n > 0) then + ! derivatives in the x-direction + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + end if + if (p > 0) then + ! derivatives in the x-direction + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k) + ! derivatives in the z-direction + tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l) + tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l) + tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l) + end if end do + + ! STEP 2a: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1) + ! STEP 2b: computing the inverse of the grad_xi tensor + tensorb(1) = 1._wp/(tensora(1)**2) + if (n > 0) then + ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse + tensorb(1) = tensora(4) + tensorb(2) = -tensora(3) + tensorb(3) = -tensora(2) + tensorb(4) = tensora(1) + ! STEP 2b: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3) + end if + if (p > 0) then + ! 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)) + end if + + if (tensorb(tensor_size) > verysmall) then + ! STEP 2c: computing the inverse of grad_xi tensor = F + ! tensorb is the adjoint, tensora becomes F + !$acc loop 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: update the btensor, this is consistent with Riemann solvers + ! \b_xx + btensor%vf(1)%sf(j, k, l) = tensorb(1) + if (n > 0) then + ! STEP 2e: override adjoint (tensorb) to be F transpose F + tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensorb(4) = tensora(3)**2 + tensora(4)**2 + tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4) + tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2) + ! STEP 3: update the btensor, this is consistent with Riemann solvers + #:for BIJ, TXY in [(1,1),(2,2),(3,4)] + btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) + #:endfor + end if + if (p > 0) then + ! STEP 2e: 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 3a: 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 + end if + + !STEP 3b: store the determinant at the last entry of the btensor + btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) + + ! STEP 4a: updating the Cauchy stress primitive scalar field + if (hyper_model == 1) then + call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + elseif (hyper_model == 2) then + call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + end if + + ! STEP 4b: updating the pressure field + q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & + G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma + + ! STEP 4c: updating the Cauchy stress conservative scalar field + !$acc loop seq + do i = 1, b_size - 1 + q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) + end do + end if end if - end if + end do end do end do - end do - !$acc end parallel loop - end if + !$acc end parallel loop + end subroutine s_hyperelastic_rmt_stress_update !> The following subroutine handles the calculation of the btensor @@ -383,7 +266,7 @@ contains !! calculate the inverse of grad_xi to obtain F, F is a nxn tensor !! calculate the FFtranspose to obtain the btensor, btensor is nxn tensor !! btensor is symmetric, save the data space - subroutine s_neoHookean_cauchy_solver_1D(btensor, q_prim_vf, G, j, k, l) + subroutine s_neoHookean_cauchy_solver(btensor, q_prim_vf, G, j, k, l) !$acc routine seq type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf type(scalar_field), dimension(b_size), intent(inout) :: btensor @@ -396,75 +279,26 @@ contains ! tensor is the symmetric tensor & calculate the trace of the tensor trace = btensor(1)%sf(j, k, l) - - ! calculate the deviatoric of the tensor - btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace - - ! dividing by the jacobian for neo-Hookean model - ! setting the tensor to the stresses for riemann solver - !$acc loop seq - do i = 1, b_size - 1 - q_prim_vf(strxb + i - 1)%sf(j, k, l) = & - G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) - end do - ! compute the invariant without the elastic modulus - q_prim_vf(xiend + 1)%sf(j, k, l) = & - 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) - - end subroutine s_neoHookean_cauchy_solver_1D - - subroutine s_neoHookean_cauchy_solver_2D(btensor, q_prim_vf, G, j, k, l) - !$acc routine seq - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(b_size), intent(inout) :: btensor - real(wp), intent(in) :: G - integer, intent(in) :: j, k, l - - real(wp) :: trace - real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i - - ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) - ! calculate the deviatoric of the tensor btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace - btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace - - ! dividing by the jacobian for neo-Hookean model - ! setting the tensor to the stresses for riemann solver - !$acc loop seq - do i = 1, b_size - 1 - q_prim_vf(strxb + i - 1)%sf(j, k, l) = & - G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) - end do - ! compute the invariant without the elastic modulus - q_prim_vf(xiend + 1)%sf(j, k, l) = & - 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) - - end subroutine s_neoHookean_cauchy_solver_2D - - subroutine s_neoHookean_cauchy_solver_3D(btensor, q_prim_vf, G, j, k, l) - !$acc routine seq - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(b_size), intent(inout) :: btensor - real(wp), intent(in) :: G - integer, intent(in) :: j, k, l - - real(wp) :: trace - real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i - - ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) + if (n > 0) then + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace + btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace + end if + if (p > 0) then + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) + ! calculate the deviatoric of the tensor + #:for IJ in [1,3,6] + btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace + #:endfor + end if - ! calculate the deviatoric of the tensor - #:for IJ in [1,3,6] - btensor(${IJ}$)%sf(j, k, l) = btensor(${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 - !$acc loop seq do i = 1, b_size - 1 q_prim_vf(strxb + i - 1)%sf(j, k, l) = & @@ -474,7 +308,7 @@ contains q_prim_vf(xiend + 1)%sf(j, k, l) = & 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) - end subroutine s_neoHookean_cauchy_solver_3D + end subroutine s_neoHookean_cauchy_solver !> The following subroutine handles the calculation of the btensor !! with a Mooney-Rivlin material model. diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index bce7956e3..b00f4a0b0 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -667,7 +667,7 @@ contains call nvtxStartRange("RHS-ELASTIC") if (hyperelasticity) then - call s_hyperelastic_rmt_stress_update(num_dims, q_cons_qp%vf, q_prim_qp%vf) + call s_hyperelastic_rmt_stress_update(q_cons_qp%vf, q_prim_qp%vf) call s_populate_variables_buffers(q_prim_qp%vf, pb, mv) end if call nvtxEndRange diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 8c71f2474..40d907e95 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -846,7 +846,7 @@ contains if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) call nvtxStartRange("RHS-ELASTIC") - if (hyperelasticity) call s_hyperelastic_rmt_stress_update(num_dims, q_cons_ts(1)%vf, q_prim_vf) + if (hyperelasticity) call s_hyperelastic_rmt_stress_update(q_cons_ts(1)%vf, q_prim_vf) call nvtxEndRange if (ib) then From be859c5f4f87f8e9cae07f399e49412d1d68feda Mon Sep 17 00:00:00 2001 From: mcarcana Date: Wed, 29 Jan 2025 11:49:04 -0500 Subject: [PATCH 8/9] returning an accidental deletion from conflict resolution --- src/simulation/m_hyperelastic.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index ba304004f..d0f5630eb 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -126,8 +126,8 @@ contains if (G > verysmall) then !$acc loop seq - do i = 1, tensor_size - tensora(i) = 0._wp + do i = 1, tensor_size - 1 + tensora(i) = tensorb(i)/tensorb(tensor_size) end do ! STEP 1: computing the grad_xi tensor using finite differences From 0dca297b488d760c6f699b472c6f9400258aa513 Mon Sep 17 00:00:00 2001 From: Mauro Rodriguez Date: Wed, 29 Jan 2025 17:10:57 -0500 Subject: [PATCH 9/9] further cleaning to the code --- src/simulation/m_hyperelastic.fpp | 294 +++++++++++++++--------------- 1 file changed, 147 insertions(+), 147 deletions(-) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index d0f5630eb..c005feeb7 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -108,152 +108,152 @@ contains !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, & !$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb) - do l = 0, p + do l = 0, p do k = 0, n - do j = 0, m - !$acc loop seq - do i = 1, num_fluids - alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) - alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - - ! 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, j, k, l, G, Gs) - rho = max(rho, sgm_eps) - G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0._wp - - if (G > verysmall) then + do j = 0, m + !$acc loop seq + do i = 1, num_fluids + alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) + alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) + end do + + ! 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, j, k, l, G, Gs) + rho = max(rho, sgm_eps) + G = max(G, sgm_eps) + !if ( G <= verysmall ) G_K = 0._wp + + if (G > verysmall) then + !$acc loop seq + do i = 1, tensor_size - 1 + tensora(i) = tensorb(i)/tensorb(tensor_size) + end do + + ! STEP 1: computing the grad_xi tensor using finite differences + ! grad_xi definition / organization + ! tensora(1,2,3): dxix_dx, dxiy_dx, dxiz_dx + ! tensora(4,5,6): dxix_dy, dxiy_dy, dxiz_dy + ! tensora(7,8,9): dxix_dz, dxiy_dz, dxiz_dz + !$acc loop seq + do r = -fd_number, fd_number + ! derivatives in the x-direction + tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) + if (n > 0) then + ! derivatives in the x-direction + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + end if + if (p > 0) then + ! derivatives in the x-direction + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k) + ! derivatives in the z-direction + tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l) + tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l) + tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l) + end if + end do + + ! STEP 2a: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1) + ! STEP 2b: computing the inverse of the grad_xi tensor + tensorb(1) = 1._wp/(tensora(1)**2) + if (n > 0) then + ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse + tensorb(1) = tensora(4) + tensorb(2) = -tensora(3) + tensorb(3) = -tensora(2) + tensorb(4) = tensora(1) + ! STEP 2b: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3) + end if + if (p > 0) then + ! 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)) + end if + + if (tensorb(tensor_size) > verysmall) then + ! STEP 2c: computing the inverse of grad_xi tensor = F + ! tensorb is the adjoint, tensora becomes F !$acc loop seq do i = 1, tensor_size - 1 tensora(i) = tensorb(i)/tensorb(tensor_size) end do - ! STEP 1: computing the grad_xi tensor using finite differences - ! grad_xi definition / organization - ! tensora(1,2,3): dxix_dx, dxiy_dx, dxiz_dx - ! tensora(4,5,6): dxix_dy, dxiy_dy, dxiz_dy - ! tensora(7,8,9): dxix_dz, dxiy_dz, dxiz_dz - !$acc loop seq - do r = -fd_number, fd_number - ! derivatives in the x-direction - tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) - if (n > 0) then - ! derivatives in the x-direction - tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) - ! derivatives in the y-direction - tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) - end if - if (p > 0) then - ! derivatives in the x-direction - tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) - tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j) - ! derivatives in the y-direction - tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k) - ! derivatives in the z-direction - tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l) - tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l) - tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l) - end if - end do - - ! STEP 2a: computing the determinant of the grad_xi tensor - tensorb(tensor_size) = tensora(1) - ! STEP 2b: computing the inverse of the grad_xi tensor - tensorb(1) = 1._wp/(tensora(1)**2) + ! STEP 2d: computing the J = det(F) = 1/det(\grad{\xi}) + tensorb(tensor_size) = 1._wp/tensorb(tensor_size) + ! STEP 3: update the btensor, this is consistent with Riemann solvers + ! \b_xx + btensor%vf(1)%sf(j, k, l) = tensorb(1) if (n > 0) then - ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse - tensorb(1) = tensora(4) - tensorb(2) = -tensora(3) - tensorb(3) = -tensora(2) - tensorb(4) = tensora(1) - ! STEP 2b: computing the determinant of the grad_xi tensor - tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3) + ! STEP 2e: override adjoint (tensorb) to be F transpose F + tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensorb(4) = tensora(3)**2 + tensora(4)**2 + tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4) + tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2) + ! STEP 3: update the btensor, this is consistent with Riemann solvers + #:for BIJ, TXY in [(1,1),(2,2),(3,4)] + btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) + #:endfor end if if (p > 0) then - ! 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)) - end if - - if (tensorb(tensor_size) > verysmall) then - ! STEP 2c: computing the inverse of grad_xi tensor = F - ! tensorb is the adjoint, tensora becomes F - !$acc loop 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: update the btensor, this is consistent with Riemann solvers - ! \b_xx - btensor%vf(1)%sf(j, k, l) = tensorb(1) - if (n > 0) then - ! STEP 2e: override adjoint (tensorb) to be F transpose F - tensorb(1) = tensora(1)**2 + tensora(2)**2 - tensorb(4) = tensora(3)**2 + tensora(4)**2 - tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4) - tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2) - ! STEP 3: update the btensor, this is consistent with Riemann solvers - #:for BIJ, TXY in [(1,1),(2,2),(3,4)] - btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) - #:endfor - end if - if (p > 0) then - ! STEP 2e: 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 3a: 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 - end if - - !STEP 3b: store the determinant at the last entry of the btensor - btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) - - ! STEP 4a: updating the Cauchy stress primitive scalar field - if (hyper_model == 1) then - call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) - elseif (hyper_model == 2) then - call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) - end if - - ! STEP 4b: updating the pressure field - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma - - ! STEP 4c: updating the Cauchy stress conservative scalar field - !$acc loop seq - do i = 1, b_size - 1 - q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) - end do + ! STEP 2e: 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 3a: 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 end if + + !STEP 3b: store the determinant at the last entry of the btensor + btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) + + ! STEP 4a: updating the Cauchy stress primitive scalar field + if (hyper_model == 1) then + call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + elseif (hyper_model == 2) then + call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + end if + + ! STEP 4b: updating the pressure field + q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & + G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma + + ! STEP 4c: updating the Cauchy stress conservative scalar field + !$acc loop seq + do i = 1, b_size - 1 + q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) + end do end if - end do + end if end do end do - !$acc end parallel loop + end do + !$acc end parallel loop end subroutine s_hyperelastic_rmt_stress_update @@ -282,20 +282,20 @@ contains ! calculate the deviatoric of the tensor btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace if (n > 0) then - ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) - ! calculate the deviatoric of the tensor - btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace - btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace - end if - if (p > 0) then - ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) - ! calculate the deviatoric of the tensor - #:for IJ in [1,3,6] - btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace - #:endfor - end if + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace + btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace + end if + if (p > 0) then + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) + ! calculate the deviatoric of the tensor + #:for IJ in [1,3,6] + btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace + #:endfor + end if ! dividing by the jacobian for neo-Hookean model ! setting the tensor to the stresses for riemann solver