1d/2d hyper (#6)

src/simulation/m_hyperelastic.fpp

@@ -108,112 +108,158 @@ 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
                 end do
             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