88module m_hypoelastic
99
1010 use m_derived_types !< Definitions of the derived types
11-
1211 use m_global_parameters !< Definitions of the global parameters
13-
1412 use m_finite_differences
15-
1613 use m_helper
17-
1814 use m_mpi_proxy !< Message passing interface (MPI) module proxy
1915
2016 implicit none
2117
2218 private; public :: s_initialize_hypoelastic_module, &
19+ s_finalize_hypoelastic_module, &
2320 s_compute_hypoelastic_rhs
2421
2522 real (wp), allocatable, dimension (:) :: Gs
@@ -33,11 +30,16 @@ module m_hypoelastic
3330 real (wp), allocatable, dimension (:, :, :) :: rho_K_field, G_K_field
3431 !$acc declare create(rho_K_field, G_K_field)
3532
33+ real (wp), allocatable, dimension (:, :) :: fd_coeff_x_h
34+ real (wp), allocatable, dimension (:, :) :: fd_coeff_y_h
35+ real (wp), allocatable, dimension (:, :) :: fd_coeff_z_h
36+ !$acc declare create(fd_coeff_x_h,fd_coeff_y_h,fd_coeff_z_h)
37+
3638contains
3739
3840 subroutine s_initialize_hypoelastic_module
3941
40- integer :: i
42+ integer :: i, k, r
4143
4244 @:ALLOCATE(Gs(1 :num_fluids))
4345 @:ALLOCATE(rho_K_field(0 :m,0 :n,0 :p), G_K_field(0 :m,0 :n,0 :p))
@@ -55,6 +57,29 @@ contains
5557 end do
5658 !$acc update device(Gs)
5759
60+ @:ALLOCATE(fd_coeff_x_h(- fd_number:fd_number, 0 :m))
61+ if (n > 0 ) then
62+ @:ALLOCATE(fd_coeff_y_h(- fd_number:fd_number, 0 :n))
63+ end if
64+ if (p > 0 ) then
65+ @:ALLOCATE(fd_coeff_z_h(- fd_number:fd_number, 0 :p))
66+ end if
67+
68+ ! Computing centered finite difference coefficients
69+ call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_h, buff_size, &
70+ fd_number, fd_order)
71+ !$acc update device(fd_coeff_x_h)
72+ if (n > 0 ) then
73+ call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_h, buff_size, &
74+ fd_number, fd_order)
75+ !$acc update device(fd_coeff_y_h)
76+ end if
77+ if (p > 0 ) then
78+ call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_h, buff_size, &
79+ fd_number, fd_order)
80+ !$acc update device(fd_coeff_z_h)
81+ end if
82+
5883 end subroutine s_initialize_hypoelastic_module
5984
6085 !> The purpose of this procedure is to compute the source terms
@@ -70,7 +95,7 @@ contains
7095
7196 real (wp) :: rho_K, G_K
7297
73- integer :: i, k, l, q !< Loop variables
98+ integer :: i, k, l, q, r !< Loop variables
7499 integer :: ndirs !< Number of coordinate directions
75100
76101 ndirs = 1 ; if (n > 0 ) ndirs = 2 ; if (p > 0 ) ndirs = 3
@@ -83,82 +108,91 @@ contains
83108 do q = 0 , p
84109 do l = 0 , n
85110 do k = 0 , m
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))
111+ du_dx(k, l, q) = 0._wp
112+ end do
113+ end do
114+ end do
115+ !$acc end parallel loop
116+
117+ !$acc parallel loop collapse(3 ) gang vector default(present)
118+ do q = 0 , p
119+ do l = 0 , n
120+ do k = 0 , m
121+ !$acc loop seq
122+ do r = - fd_number, fd_number
123+ du_dx(k, l, q) = du_dx(k, l, q) &
124+ + q_prim_vf(momxb)%sf(k + r, l, q)* fd_coeff_x_h(r, k)
125+ end do
126+
92127 end do
93128 end do
94129 end do
130+ !$acc end parallel loop
95131
96132 if (ndirs > 1 ) then
97133 !$acc parallel loop collapse(3 ) gang vector default(present)
98134 do q = 0 , p
99135 do l = 0 , n
100136 do k = 0 , m
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))
137+ du_dy(k, l, q) = 0._wp ; dv_dx(k, l, q) = 0._wp ; dv_dy(k, l, q) = 0._wp
138+ end do
139+ end do
140+ end do
141+ !$acc end parallel loop
142+
143+ !$acc parallel loop collapse(3 ) gang vector default(present)
144+ do q = 0 , p
145+ do l = 0 , n
146+ do k = 0 , m
147+ !$acc loop seq
148+ do r = - fd_number, fd_number
149+ du_dy(k, l, q) = du_dy(k, l, q) &
150+ + q_prim_vf(momxb)%sf(k, l + r, q)* fd_coeff_y_h(r, l)
151+ dv_dx(k, l, q) = dv_dx(k, l, q) &
152+ + q_prim_vf(momxb + 1 )%sf(k + r, l, q)* fd_coeff_x_h(r, k)
153+ dv_dy(k, l, q) = dv_dy(k, l, q) &
154+ + q_prim_vf(momxb + 1 )%sf(k, l + r, q)* fd_coeff_y_h(r, l)
155+ end do
119156 end do
120157 end do
121158 end do
159+ !$acc end parallel loop
122160
123161 ! 3D
124162 if (ndirs == 3 ) then
163+
164+ !$acc parallel loop collapse(3 ) gang vector default(present)
165+ do q = 0 , p
166+ do l = 0 , n
167+ do k = 0 , m
168+ du_dz(k, l, q) = 0_wp ; dv_dz(k, l, q) = 0_wp ; dw_dx(k, l, q) = 0_wp ;
169+ dw_dy(k, l, q) = 0_wp ; dw_dz(k, l, q) = 0_wp ;
170+ end do
171+ end do
172+ end do
173+ !$acc end parallel loop
174+
125175 !$acc parallel loop collapse(3 ) gang vector default(present)
126176 do q = 0 , p
127177 do l = 0 , n
128178 do k = 0 , m
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))
179+ !$acc loop seq
180+ do r = - fd_number, fd_number
181+ du_dz(k, l, q) = du_dz(k, l, q) &
182+ + q_prim_vf(momxb)%sf(k, l, q + r)* fd_coeff_z_h(r, q)
183+ dv_dz(k, l, q) = dv_dz(k, l, q) &
184+ + q_prim_vf(momxb + 1 )%sf(k, l, q + r)* fd_coeff_z_h(r, q)
185+ dw_dx(k, l, q) = dw_dx(k, l, q) &
186+ + q_prim_vf(momxe)%sf(k + r, l, q)* fd_coeff_x_h(r, k)
187+ dw_dy(k, l, q) = dw_dy(k, l, q) &
188+ + q_prim_vf(momxe)%sf(k, l + r, q)* fd_coeff_y_h(r, l)
189+ dw_dz(k, l, q) = dw_dz(k, l, q) &
190+ + q_prim_vf(momxe)%sf(k, l, q + r)* fd_coeff_z_h(r, q)
191+ end do
159192 end do
160193 end do
161194 end do
195+ !$acc end parallel loop
162196 end if
163197 end if
164198
@@ -175,7 +209,7 @@ contains
175209 G_K_field(k, l, q) = G_K
176210
177211 !TODO: take this out if not needed
178- if (G_K < 1000 ) then
212+ if (G_K < verysmall ) then
179213 G_K_field(k, l, q) = 0
180214 end if
181215 end do
@@ -300,4 +334,21 @@ contains
300334
301335 end subroutine s_compute_hypoelastic_rhs
302336
337+ subroutine s_finalize_hypoelastic_module ()
338+
339+ @:DEALLOCATE(Gs)
340+ @:DEALLOCATE(rho_K_field, G_K_field)
341+ @:DEALLOCATE(du_dx)
342+ @:DEALLOCATE(fd_coeff_x_h)
343+ if (n > 0 ) then
344+ @:DEALLOCATE(du_dy,dv_dx,dv_dy)
345+ @:DEALLOCATE(fd_coeff_y_h)
346+ if (p > 0 ) then
347+ @:DEALLOCATE(du_dz, dv_dz, dw_dx, dw_dy, dw_dz)
348+ @:DEALLOCATE(fd_coeff_z_h)
349+ end if
350+ end if
351+
352+ end subroutine s_finalize_hypoelastic_module
353+
303354end module m_hypoelastic
0 commit comments