Skip to content

Commit 630695c

Browse files
ChrisZYJChrisZYJ
authored
Fix hypoelastic instability (#773)
1 parent 8c288fd commit 630695c

File tree

8 files changed

+564
-167
lines changed

8 files changed

+564
-167
lines changed

‎src/simulation/m_hypoelastic.fpp‎

Lines changed: 59 additions & 114 deletions
Original file line numberDiff line numberDiff line change
@@ -13,14 +13,13 @@ module m_hypoelastic
1313

1414
use m_global_parameters !< Definitions of the global parameters
1515

16-
use m_finite_differences
17-
use m_helper
16+
use m_mpi_proxy !< Message passing interface (MPI) module proxy
17+
1818
! ==========================================================================
1919

2020
implicit none
2121

2222
private; public :: s_initialize_hypoelastic_module, &
23-
s_finalize_hypoelastic_module, &
2423
s_compute_hypoelastic_rhs
2524

2625
real(wp), allocatable, dimension(:) :: Gs
@@ -34,16 +33,11 @@ module m_hypoelastic
3433
real(wp), allocatable, dimension(:, :, :) :: rho_K_field, G_K_field
3534
!$acc declare create(rho_K_field, G_K_field)
3635

37-
real(wp), allocatable, dimension(:, :) :: fd_coeff_x_h
38-
real(wp), allocatable, dimension(:, :) :: fd_coeff_y_h
39-
real(wp), allocatable, dimension(:, :) :: fd_coeff_z_h
40-
!$acc declare create(fd_coeff_x_h,fd_coeff_y_h,fd_coeff_z_h)
41-
4236
contains
4337

4438
subroutine s_initialize_hypoelastic_module
4539

46-
integer :: i, k, r
40+
integer :: i
4741

4842
@:ALLOCATE(Gs(1:num_fluids))
4943
@:ALLOCATE(rho_K_field(0:m,0:n,0:p), G_K_field(0:m,0:n,0:p))
@@ -61,29 +55,6 @@ contains
6155
end do
6256
!$acc update device(Gs)
6357

64-
@:ALLOCATE(fd_coeff_x_h(-fd_number:fd_number, 0:m))
65-
if (n > 0) then
66-
@:ALLOCATE(fd_coeff_y_h(-fd_number:fd_number, 0:n))
67-
end if
68-
if (p > 0) then
69-
@:ALLOCATE(fd_coeff_z_h(-fd_number:fd_number, 0:p))
70-
end if
71-
72-
! Computing centered finite difference coefficients
73-
call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_h, buff_size, &
74-
fd_number, fd_order)
75-
!$acc update device(fd_coeff_x_h)
76-
if (n > 0) then
77-
call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_h, buff_size, &
78-
fd_number, fd_order)
79-
!$acc update device(fd_coeff_y_h)
80-
end if
81-
if (p > 0) then
82-
call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_h, buff_size, &
83-
fd_number, fd_order)
84-
!$acc update device(fd_coeff_z_h)
85-
end if
86-
8758
end subroutine s_initialize_hypoelastic_module
8859

8960
!> The purpose of this procedure is to compute the source terms
@@ -99,7 +70,7 @@ contains
9970

10071
real(wp) :: rho_K, G_K
10172

102-
integer :: i, k, l, q, r !< Loop variables
73+
integer :: i, k, l, q !< Loop variables
10374
integer :: ndirs !< Number of coordinate directions
10475

10576
ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3
@@ -112,91 +83,82 @@ contains
11283
do q = 0, p
11384
do l = 0, n
11485
do k = 0, m
115-
du_dx(k, l, q) = 0._wp
86+
du_dx(k, l, q) = &
87+
(q_prim_vf(momxb)%sf(k - 2, l, q) &
88+
- 8._wp*q_prim_vf(momxb)%sf(k - 1, l, q) &
89+
+ 8._wp*q_prim_vf(momxb)%sf(k + 1, l, q) &
90+
- q_prim_vf(momxb)%sf(k + 2, l, q)) &
91+
/(12._wp*dx(k))
11692
end do
11793
end do
11894
end do
119-
!$acc end parallel loop
120-
121-
!$acc parallel loop collapse(3) gang vector default(present)
122-
do q = 0, p
123-
do l = 0, n
124-
do k = 0, m
125-
!$acc loop seq
126-
do r = -fd_number, fd_number
127-
du_dx(k, l, q) = du_dx(k, l, q) &
128-
+ q_prim_vf(momxb)%sf(k + r, l, q)*fd_coeff_x_h(r, k)
129-
end do
130-
131-
end do
132-
end do
133-
end do
134-
!$acc end parallel loop
13595

13696
if (ndirs > 1) then
13797
!$acc parallel loop collapse(3) gang vector default(present)
13898
do q = 0, p
13999
do l = 0, n
140100
do k = 0, m
141-
du_dy(k, l, q) = 0._wp; dv_dx(k, l, q) = 0._wp; dv_dy(k, l, q) = 0._wp
101+
du_dy(k, l, q) = &
102+
(q_prim_vf(momxb)%sf(k, l - 2, q) &
103+
- 8._wp*q_prim_vf(momxb)%sf(k, l - 1, q) &
104+
+ 8._wp*q_prim_vf(momxb)%sf(k, l + 1, q) &
105+
- q_prim_vf(momxb)%sf(k, l + 2, q)) &
106+
/(12._wp*dy(l))
107+
dv_dx(k, l, q) = &
108+
(q_prim_vf(momxb + 1)%sf(k - 2, l, q) &
109+
- 8._wp*q_prim_vf(momxb + 1)%sf(k - 1, l, q) &
110+
+ 8._wp*q_prim_vf(momxb + 1)%sf(k + 1, l, q) &
111+
- q_prim_vf(momxb + 1)%sf(k + 2, l, q)) &
112+
/(12._wp*dx(k))
113+
dv_dy(k, l, q) = &
114+
(q_prim_vf(momxb + 1)%sf(k, l - 2, q) &
115+
- 8._wp*q_prim_vf(momxb + 1)%sf(k, l - 1, q) &
116+
+ 8._wp*q_prim_vf(momxb + 1)%sf(k, l + 1, q) &
117+
- q_prim_vf(momxb + 1)%sf(k, l + 2, q)) &
118+
/(12._wp*dy(l))
142119
end do
143120
end do
144121
end do
145-
!$acc end parallel loop
146-
147-
!$acc parallel loop collapse(3) gang vector default(present)
148-
do q = 0, p
149-
do l = 0, n
150-
do k = 0, m
151-
!$acc loop seq
152-
do r = -fd_number, fd_number
153-
du_dy(k, l, q) = du_dy(k, l, q) &
154-
+ q_prim_vf(momxb)%sf(k, l + r, q)*fd_coeff_y_h(r, l)
155-
dv_dx(k, l, q) = dv_dx(k, l, q) &
156-
+ q_prim_vf(momxb + 1)%sf(k + r, l, q)*fd_coeff_x_h(r, k)
157-
dv_dy(k, l, q) = dv_dy(k, l, q) &
158-
+ q_prim_vf(momxb + 1)%sf(k, l + r, q)*fd_coeff_y_h(r, l)
159-
end do
160-
end do
161-
end do
162-
end do
163-
!$acc end parallel loop
164122

165123
! 3D
166124
if (ndirs == 3) then
167-
168-
!$acc parallel loop collapse(3) gang vector default(present)
169-
do q = 0, p
170-
do l = 0, n
171-
do k = 0, m
172-
du_dz(k, l, q) = 0_wp; dv_dz(k, l, q) = 0_wp; dw_dx(k, l, q) = 0_wp;
173-
dw_dy(k, l, q) = 0_wp; dw_dz(k, l, q) = 0_wp;
174-
end do
175-
end do
176-
end do
177-
!$acc end parallel loop
178-
179125
!$acc parallel loop collapse(3) gang vector default(present)
180126
do q = 0, p
181127
do l = 0, n
182128
do k = 0, m
183-
!$acc loop seq
184-
do r = -fd_number, fd_number
185-
du_dz(k, l, q) = du_dz(k, l, q) &
186-
+ q_prim_vf(momxb)%sf(k, l, q + r)*fd_coeff_z_h(r, q)
187-
dv_dz(k, l, q) = dv_dz(k, l, q) &
188-
+ q_prim_vf(momxb + 1)%sf(k, l, q + r)*fd_coeff_z_h(r, q)
189-
dw_dx(k, l, q) = dw_dx(k, l, q) &
190-
+ q_prim_vf(momxe)%sf(k + r, l, q)*fd_coeff_x_h(r, k)
191-
dw_dy(k, l, q) = dw_dy(k, l, q) &
192-
+ q_prim_vf(momxe)%sf(k, l + r, q)*fd_coeff_y_h(r, l)
193-
dw_dz(k, l, q) = dw_dz(k, l, q) &
194-
+ q_prim_vf(momxe)%sf(k, l, q + r)*fd_coeff_z_h(r, q)
195-
end do
129+
du_dz(k, l, q) = &
130+
(q_prim_vf(momxb)%sf(k, l, q - 2) &
131+
- 8._wp*q_prim_vf(momxb)%sf(k, l, q - 1) &
132+
+ 8._wp*q_prim_vf(momxb)%sf(k, l, q + 1) &
133+
- q_prim_vf(momxb)%sf(k, l, q + 2)) &
134+
/(12._wp*dz(q))
135+
dv_dz(k, l, q) = &
136+
(q_prim_vf(momxb + 1)%sf(k, l, q - 2) &
137+
- 8._wp*q_prim_vf(momxb + 1)%sf(k, l, q - 1) &
138+
+ 8._wp*q_prim_vf(momxb + 1)%sf(k, l, q + 1) &
139+
- q_prim_vf(momxb + 1)%sf(k, l, q + 2)) &
140+
/(12._wp*dz(q))
141+
dw_dx(k, l, q) = &
142+
(q_prim_vf(momxe)%sf(k - 2, l, q) &
143+
- 8._wp*q_prim_vf(momxe)%sf(k - 1, l, q) &
144+
+ 8._wp*q_prim_vf(momxe)%sf(k + 1, l, q) &
145+
- q_prim_vf(momxe)%sf(k + 2, l, q)) &
146+
/(12._wp*dx(k))
147+
dw_dy(k, l, q) = &
148+
(q_prim_vf(momxe)%sf(k, l - 2, q) &
149+
- 8._wp*q_prim_vf(momxe)%sf(k, l - 1, q) &
150+
+ 8._wp*q_prim_vf(momxe)%sf(k, l + 1, q) &
151+
- q_prim_vf(momxe)%sf(k, l + 2, q)) &
152+
/(12._wp*dy(l))
153+
dw_dz(k, l, q) = &
154+
(q_prim_vf(momxe)%sf(k, l, q - 2) &
155+
- 8._wp*q_prim_vf(momxe)%sf(k, l, q - 1) &
156+
+ 8._wp*q_prim_vf(momxe)%sf(k, l, q + 1) &
157+
- q_prim_vf(momxe)%sf(k, l, q + 2)) &
158+
/(12._wp*dz(q))
196159
end do
197160
end do
198161
end do
199-
!$acc end parallel loop
200162
end if
201163
end if
202164

@@ -213,7 +175,7 @@ contains
213175
G_K_field(k, l, q) = G_K
214176

215177
!TODO: take this out if not needed
216-
if (G_K < verysmall) then
178+
if (G_K < 1000) then
217179
G_K_field(k, l, q) = 0
218180
end if
219181
end do
@@ -338,21 +300,4 @@ contains
338300

339301
end subroutine s_compute_hypoelastic_rhs
340302

341-
subroutine s_finalize_hypoelastic_module() ! --------------------
342-
343-
@:DEALLOCATE(Gs)
344-
@:DEALLOCATE(rho_K_field, G_K_field)
345-
@:DEALLOCATE(du_dx)
346-
@:DEALLOCATE(fd_coeff_x_h)
347-
if (n > 0) then
348-
@:DEALLOCATE(du_dy,dv_dx,dv_dy)
349-
@:DEALLOCATE(fd_coeff_y_h)
350-
if (p > 0) then
351-
@:DEALLOCATE(du_dz, dv_dz, dw_dx, dw_dy, dw_dz)
352-
@:DEALLOCATE(fd_coeff_z_h)
353-
end if
354-
end if
355-
356-
end subroutine s_finalize_hypoelastic_module
357-
358303
end module m_hypoelastic

‎src/simulation/m_start_up.fpp‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1636,7 +1636,6 @@ contains
16361636
subroutine s_finalize_modules
16371637

16381638
call s_finalize_time_steppers_module()
1639-
if (hypoelasticity) call s_finalize_hypoelastic_module()
16401639
if (hyperelasticity) call s_finalize_hyperelastic_module()
16411640
call s_finalize_derived_variables_module()
16421641
call s_finalize_data_output_module()

0 commit comments

Comments
 (0)