diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index a28a0191e..c005feeb7 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -108,101 +108,145 @@ 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 - 2 - do k = 0, n - 2 - do j = 2, m - 2 + 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 ) G_K = 0._wp 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 + ! 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 + ! 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) - 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) + 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 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)) + + ! 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 + ! 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: 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 + ! 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 5a: updating the Cauchy stress primitive scalar field + + ! 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 5b: updating the pressure field + + ! 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 5c: updating the Cauchy stress conservative scalar field + + ! 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) + 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 @@ -210,10 +254,12 @@ contains end do end do !$acc end parallel loop + 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 @@ -232,15 +278,27 @@ contains 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) - + trace = btensor(1)%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 + 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 + ! 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) = & @@ -252,7 +310,8 @@ contains end subroutine s_neoHookean_cauchy_solver - !> 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