Skip to content

Commit 44eb445

Browse files
mrodrig6mrodrig6
authored
1d/2d hyper (#6)
2 parents b6c2863 + 0dca297 commit 44eb445

File tree

1 file changed

+124
-65
lines changed

1 file changed

+124
-65
lines changed

‎src/simulation/m_hyperelastic.fpp‎

Lines changed: 124 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -108,112 +108,158 @@ contains
108108

109109
!$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, &
110110
!$acc rho, gamma, pi_inf, qv, G, Re, tensora, tensorb)
111-
do l = 0, p - 2
112-
do k = 0, n - 2
113-
do j = 2, m - 2
111+
do l = 0, p
112+
do k = 0, n
113+
do j = 0, m
114114
!$acc loop seq
115115
do i = 1, num_fluids
116116
alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l)
117117
alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l)
118118
end do
119+
119120
! If in simulation, use acc mixture subroutines
120121
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, &
121122
alpha_rho_k, Re, j, k, l, G, Gs)
122123
rho = max(rho, sgm_eps)
123124
G = max(G, sgm_eps)
124-
!if ( G <= verysmall ) G_K = 0_wp
125+
!if ( G <= verysmall ) G_K = 0._wp
125126

126127
if (G > verysmall) then
127128
!$acc loop seq
128-
do i = 1, tensor_size
129-
tensora(i) = 0_wp
129+
do i = 1, tensor_size - 1
130+
tensora(i) = tensorb(i)/tensorb(tensor_size)
130131
end do
131132

132-
! Step 1: computing the grad_xi tensor using finite differences
133+
! STEP 1: computing the grad_xi tensor using finite differences
133134
! grad_xi definition / organization
134-
! number for the tensor 1-3: dxix_dx, dxiy_dx, dxiz_dx
135-
! 4-6 : dxix_dy, dxiy_dy, dxiz_dy
136-
! 7-9 : dxix_dz, dxiy_dz, dxiz_dz
135+
! tensora(1,2,3): dxix_dx, dxiy_dx, dxiz_dx
136+
! tensora(4,5,6): dxix_dy, dxiy_dy, dxiz_dy
137+
! tensora(7,8,9): dxix_dz, dxiy_dz, dxiz_dz
137138
!$acc loop seq
138139
do r = -fd_number, fd_number
139140
! derivatives in the x-direction
140141
tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j)
141-
tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j)
142-
tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j)
143-
! derivatives in the y-direction
144-
tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k)
145-
tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k)
146-
tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k)
147-
! derivatives in the z-direction
148-
tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l)
149-
tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l)
150-
tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l)
142+
if (n > 0) then
143+
! derivatives in the x-direction
144+
tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j)
145+
! derivatives in the y-direction
146+
tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k)
147+
tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k)
148+
end if
149+
if (p > 0) then
150+
! derivatives in the x-direction
151+
tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j)
152+
tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j)
153+
! derivatives in the y-direction
154+
tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k)
155+
tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k)
156+
tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k)
157+
! derivatives in the z-direction
158+
tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l)
159+
tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l)
160+
tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l)
161+
end if
151162
end do
152-
! Step 2a: computing the adjoint of the grad_xi tensor for the inverse
153-
tensorb(1) = tensora(5)*tensora(9) - tensora(6)*tensora(8)
154-
tensorb(2) = -(tensora(2)*tensora(9) - tensora(3)*tensora(8))
155-
tensorb(3) = tensora(2)*tensora(6) - tensora(3)*tensora(5)
156-
tensorb(4) = -(tensora(4)*tensora(9) - tensora(6)*tensora(7))
157-
tensorb(5) = tensora(1)*tensora(9) - tensora(3)*tensora(7)
158-
tensorb(6) = -(tensora(1)*tensora(6) - tensora(4)*tensora(3))
159-
tensorb(7) = tensora(4)*tensora(8) - tensora(5)*tensora(7)
160-
tensorb(8) = -(tensora(1)*tensora(8) - tensora(2)*tensora(7))
161-
tensorb(9) = tensora(1)*tensora(5) - tensora(2)*tensora(4)
162-
163-
! Step 2b: computing the determinant of the grad_xi tensor
164-
tensorb(tensor_size) = tensora(1)*(tensora(5)*tensora(9) - tensora(6)*tensora(8)) &
165-
- tensora(2)*(tensora(4)*tensora(9) - tensora(6)*tensora(7)) &
166-
+ tensora(3)*(tensora(4)*tensora(8) - tensora(5)*tensora(7))
163+
164+
! STEP 2a: computing the determinant of the grad_xi tensor
165+
tensorb(tensor_size) = tensora(1)
166+
! STEP 2b: computing the inverse of the grad_xi tensor
167+
tensorb(1) = 1._wp/(tensora(1)**2)
168+
if (n > 0) then
169+
! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse
170+
tensorb(1) = tensora(4)
171+
tensorb(2) = -tensora(3)
172+
tensorb(3) = -tensora(2)
173+
tensorb(4) = tensora(1)
174+
! STEP 2b: computing the determinant of the grad_xi tensor
175+
tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3)
176+
end if
177+
if (p > 0) then
178+
! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse
179+
tensorb(1) = tensora(5)*tensora(9) - tensora(6)*tensora(8)
180+
tensorb(2) = -(tensora(2)*tensora(9) - tensora(3)*tensora(8))
181+
tensorb(3) = tensora(2)*tensora(6) - tensora(3)*tensora(5)
182+
tensorb(4) = -(tensora(4)*tensora(9) - tensora(6)*tensora(7))
183+
tensorb(5) = tensora(1)*tensora(9) - tensora(3)*tensora(7)
184+
tensorb(6) = -(tensora(1)*tensora(6) - tensora(4)*tensora(3))
185+
tensorb(7) = tensora(4)*tensora(8) - tensora(5)*tensora(7)
186+
tensorb(8) = -(tensora(1)*tensora(8) - tensora(2)*tensora(7))
187+
tensorb(9) = tensora(1)*tensora(5) - tensora(2)*tensora(4)
188+
! STEP 2b: computing the determinant of the grad_xi tensor
189+
tensorb(tensor_size) = tensora(1)*(tensora(5)*tensora(9) - tensora(6)*tensora(8)) &
190+
- tensora(2)*(tensora(4)*tensora(9) - tensora(6)*tensora(7)) &
191+
+ tensora(3)*(tensora(4)*tensora(8) - tensora(5)*tensora(7))
192+
end if
167193

168194
if (tensorb(tensor_size) > verysmall) then
169-
! Step 2c: computing the inverse of grad_xi tensor = F
195+
! STEP 2c: computing the inverse of grad_xi tensor = F
170196
! tensorb is the adjoint, tensora becomes F
171197
!$acc loop seq
172198
do i = 1, tensor_size - 1
173199
tensora(i) = tensorb(i)/tensorb(tensor_size)
174200
end do
175201

176-
! Step 2d: computing the J = det(F) = 1/det(\grad{\xi})
177-
tensorb(tensor_size) = 1_wp/tensorb(tensor_size)
178-
179-
! Step 3: computing F transpose F
180-
tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensora(3)**2
181-
tensorb(5) = tensora(4)**2 + tensora(5)**2 + tensora(6)**2
182-
tensorb(9) = tensora(7)**2 + tensora(8)**2 + tensora(9)**2
183-
tensorb(2) = tensora(1)*tensora(4) + tensora(2)*tensora(5) + tensora(3)*tensora(6)
184-
tensorb(3) = tensora(1)*tensora(7) + tensora(2)*tensora(8) + tensora(3)*tensora(9)
185-
tensorb(6) = tensora(4)*tensora(7) + tensora(5)*tensora(8) + tensora(6)*tensora(9)
186-
! Step 4: update the btensor, this is consistent with Riemann solvers
187-
#:for BIJ, TXY in [(1,1),(2,2),(3,5),(4,3),(5,6),(6,9)]
188-
btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$)
189-
#:endfor
190-
! store the determinant at the last entry of the btensor
202+
! STEP 2d: computing the J = det(F) = 1/det(\grad{\xi})
203+
tensorb(tensor_size) = 1._wp/tensorb(tensor_size)
204+
! STEP 3: update the btensor, this is consistent with Riemann solvers
205+
! \b_xx
206+
btensor%vf(1)%sf(j, k, l) = tensorb(1)
207+
if (n > 0) then
208+
! STEP 2e: override adjoint (tensorb) to be F transpose F
209+
tensorb(1) = tensora(1)**2 + tensora(2)**2
210+
tensorb(4) = tensora(3)**2 + tensora(4)**2
211+
tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4)
212+
tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2)
213+
! STEP 3: update the btensor, this is consistent with Riemann solvers
214+
#:for BIJ, TXY in [(1,1),(2,2),(3,4)]
215+
btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$)
216+
#:endfor
217+
end if
218+
if (p > 0) then
219+
! STEP 2e: computing F transpose F
220+
tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensora(3)**2
221+
tensorb(5) = tensora(4)**2 + tensora(5)**2 + tensora(6)**2
222+
tensorb(9) = tensora(7)**2 + tensora(8)**2 + tensora(9)**2
223+
tensorb(2) = tensora(1)*tensora(4) + tensora(2)*tensora(5) + tensora(3)*tensora(6)
224+
tensorb(3) = tensora(1)*tensora(7) + tensora(2)*tensora(8) + tensora(3)*tensora(9)
225+
tensorb(6) = tensora(4)*tensora(7) + tensora(5)*tensora(8) + tensora(6)*tensora(9)
226+
! STEP 3a: update the btensor, this is consistent with Riemann solvers
227+
#:for BIJ, TXY in [(1,1),(2,2),(3,5),(4,3),(5,6),(6,9)]
228+
btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$)
229+
#:endfor
230+
end if
231+
232+
!STEP 3b: store the determinant at the last entry of the btensor
191233
btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size)
192-
! Step 5a: updating the Cauchy stress primitive scalar field
234+
235+
! STEP 4a: updating the Cauchy stress primitive scalar field
193236
if (hyper_model == 1) then
194237
call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l)
195238
elseif (hyper_model == 2) then
196239
call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l)
197240
end if
198-
! Step 5b: updating the pressure field
241+
242+
! STEP 4b: updating the pressure field
199243
q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - &
200244
G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma
201-
! Step 5c: updating the Cauchy stress conservative scalar field
245+
246+
! STEP 4c: updating the Cauchy stress conservative scalar field
202247
!$acc loop seq
203248
do i = 1, b_size - 1
204-
q_cons_vf(strxb + i - 1)%sf(j, k, l) = &
205-
rho*q_prim_vf(strxb + i - 1)%sf(j, k, l)
249+
q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l)
206250
end do
207251
end if
208252
end if
209253
end do
210254
end do
211255
end do
212256
!$acc end parallel loop
257+
213258
end subroutine s_hyperelastic_rmt_stress_update
214259

215-
!> The following subroutine handles the calculation of the btensor.
216-
!! The calculation of the btensor takes qprimvf.
260+
!> The following subroutine handles the calculation of the btensor
261+
!! with a neo-Hookean material model.
262+
!! The calculation of the btensor takes qprimvf.
217263
!! @param q_prim_vf Primitive variables
218264
!! @param btensor is the output
219265
!! calculate the grad_xi, grad_xi is a nxn tensor
@@ -232,15 +278,27 @@ contains
232278
integer :: i
233279

234280
! tensor is the symmetric tensor & calculate the trace of the tensor
235-
trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l)
236-
281+
trace = btensor(1)%sf(j, k, l)
237282
! calculate the deviatoric of the tensor
238-
#:for IJ in [1,3,6]
239-
btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace
240-
#:endfor
283+
btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace
284+
if (n > 0) then
285+
! tensor is the symmetric tensor & calculate the trace of the tensor
286+
trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l)
287+
! calculate the deviatoric of the tensor
288+
btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace
289+
btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace
290+
end if
291+
if (p > 0) then
292+
! tensor is the symmetric tensor & calculate the trace of the tensor
293+
trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l)
294+
! calculate the deviatoric of the tensor
295+
#:for IJ in [1,3,6]
296+
btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace
297+
#:endfor
298+
end if
299+
241300
! dividing by the jacobian for neo-Hookean model
242301
! setting the tensor to the stresses for riemann solver
243-
244302
!$acc loop seq
245303
do i = 1, b_size - 1
246304
q_prim_vf(strxb + i - 1)%sf(j, k, l) = &
@@ -252,7 +310,8 @@ contains
252310

253311
end subroutine s_neoHookean_cauchy_solver
254312

255-
!> The following subroutine handles the calculation of the btensor.
313+
!> The following subroutine handles the calculation of the btensor
314+
!! with a Mooney-Rivlin material model.
256315
!! The calculation of the btensor takes qprimvf.
257316
!! @param q_prim_vf Primitive variables
258317
!! @param btensor is the output

0 commit comments

Comments
 (0)