From 0dca297b488d760c6f699b472c6f9400258aa513 Mon Sep 17 00:00:00 2001 From: Mauro Rodriguez Date: Wed, 29 Jan 2025 17:10:57 -0500 Subject: [PATCH] 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