Skip to content

Commit 10f3361

Browse files
prathi-windprathi-wind
committed
Fixed case when num_gps is 0 for OpenMP
1 parent d69a62b commit 10f3361

File tree

1 file changed

+128
-127
lines changed

1 file changed

+128
-127
lines changed

‎src/simulation/m_ibm.fpp‎

Lines changed: 128 additions & 127 deletions
Original file line numberDiff line numberDiff line change
@@ -168,153 +168,154 @@ contains
168168
real(wp) :: buf
169169
type(ghost_point) :: gp
170170
type(ghost_point) :: innerp
171-
172-
#:call GPU_PARALLEL_LOOP(private='[physical_loc,dyn_pres,alpha_rho_IP, alpha_IP,pres_IP,vel_IP,vel_g,vel_norm_IP,r_IP, v_IP,pb_IP,mv_IP,nmom_IP,presb_IP,massv_IP,rho, gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, j,k,l,q]')
173-
do i = 1, num_gps
174-
175-
gp = ghost_points(i)
176-
j = gp%loc(1)
177-
k = gp%loc(2)
178-
l = gp%loc(3)
179-
patch_id = ghost_points(i)%ib_patch_id
180-
181-
! Calculate physical location of GP
182-
if (p > 0) then
183-
physical_loc = [x_cc(j), y_cc(k), z_cc(l)]
184-
else
185-
physical_loc = [x_cc(j), y_cc(k), 0._wp]
186-
end if
187-
188-
!Interpolate primitive variables at image point associated w/ GP
189-
if (bubbles_euler .and. .not. qbmm) then
190-
call s_interpolate_image_point(q_prim_vf, gp, &
191-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
192-
r_IP, v_IP, pb_IP, mv_IP)
193-
else if (qbmm .and. polytropic) then
194-
call s_interpolate_image_point(q_prim_vf, gp, &
195-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
196-
r_IP, v_IP, pb_IP, mv_IP, nmom_IP)
197-
else if (qbmm .and. .not. polytropic) then
198-
call s_interpolate_image_point(q_prim_vf, gp, &
199-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
200-
r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
201-
else
202-
call s_interpolate_image_point(q_prim_vf, gp, &
203-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
204-
end if
205-
dyn_pres = 0._wp
206-
207-
! Set q_prim_vf params at GP so that mixture vars calculated properly
208-
$:GPU_LOOP(parallelism='[seq]')
209-
do q = 1, num_fluids
210-
q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
211-
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
212-
end do
213-
214-
if (surface_tension) then
215-
q_prim_vf(c_idx)%sf(j, k, l) = c_IP
216-
end if
217-
if (model_eqns /= 4) then
218-
! If in simulation, use acc mixture subroutines
219-
if (elasticity) then
220-
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
221-
alpha_rho_IP, Re_K, G_K, Gs)
222-
else if (bubbles_euler) then
223-
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
224-
alpha_rho_IP, Re_K)
171+
if (num_gps > 0) then
172+
#:call GPU_PARALLEL_LOOP(private='[physical_loc,dyn_pres,alpha_rho_IP, alpha_IP,pres_IP,vel_IP,vel_g,vel_norm_IP,r_IP, v_IP,pb_IP,mv_IP,nmom_IP,presb_IP,massv_IP,rho, gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, j,k,l,q]')
173+
do i = 1, num_gps
174+
175+
gp = ghost_points(i)
176+
j = gp%loc(1)
177+
k = gp%loc(2)
178+
l = gp%loc(3)
179+
patch_id = ghost_points(i)%ib_patch_id
180+
181+
! Calculate physical location of GP
182+
if (p > 0) then
183+
physical_loc = [x_cc(j), y_cc(k), z_cc(l)]
225184
else
226-
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
227-
alpha_rho_IP, Re_K)
185+
physical_loc = [x_cc(j), y_cc(k), 0._wp]
228186
end if
229-
end if
230187

231-
! Calculate velocity of ghost cell
232-
if (gp%slip) then
233-
norm(1:3) = levelset_norm%sf(gp%loc(1), gp%loc(2), gp%loc(3), gp%ib_patch_id, 1:3)
234-
buf = sqrt(sum(norm**2))
235-
norm = norm/buf
236-
vel_norm_IP = sum(vel_IP*norm)*norm
237-
vel_g = vel_IP - vel_norm_IP
238-
else
239-
vel_g = 0._wp
240-
end if
241-
242-
! Set momentum
243-
$:GPU_LOOP(parallelism='[seq]')
244-
do q = momxb, momxe
245-
q_cons_vf(q)%sf(j, k, l) = rho*vel_g(q - momxb + 1)
246-
dyn_pres = dyn_pres + q_cons_vf(q)%sf(j, k, l)* &
247-
vel_g(q - momxb + 1)/2._wp
248-
end do
249-
250-
! Set continuity and adv vars
251-
$:GPU_LOOP(parallelism='[seq]')
252-
do q = 1, num_fluids
253-
q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
254-
q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
255-
end do
256-
257-
! Set color function
258-
if (surface_tension) then
259-
q_cons_vf(c_idx)%sf(j, k, l) = c_IP
260-
end if
188+
!Interpolate primitive variables at image point associated w/ GP
189+
if (bubbles_euler .and. .not. qbmm) then
190+
call s_interpolate_image_point(q_prim_vf, gp, &
191+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
192+
r_IP, v_IP, pb_IP, mv_IP)
193+
else if (qbmm .and. polytropic) then
194+
call s_interpolate_image_point(q_prim_vf, gp, &
195+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
196+
r_IP, v_IP, pb_IP, mv_IP, nmom_IP)
197+
else if (qbmm .and. .not. polytropic) then
198+
call s_interpolate_image_point(q_prim_vf, gp, &
199+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
200+
r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
201+
else
202+
call s_interpolate_image_point(q_prim_vf, gp, &
203+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
204+
end if
205+
dyn_pres = 0._wp
261206

262-
! Set Energy
263-
if (bubbles_euler) then
264-
q_cons_vf(E_idx)%sf(j, k, l) = (1 - alpha_IP(1))*(gamma*pres_IP + pi_inf + dyn_pres)
265-
else
266-
q_cons_vf(E_idx)%sf(j, k, l) = gamma*pres_IP + pi_inf + dyn_pres
267-
end if
268-
! Set bubble vars
269-
if (bubbles_euler .and. .not. qbmm) then
270-
call s_comp_n_from_prim(alpha_IP(1), r_IP, nbub, weight)
207+
! Set q_prim_vf params at GP so that mixture vars calculated properly
271208
$:GPU_LOOP(parallelism='[seq]')
272-
do q = 1, nb
273-
q_cons_vf(bubxb + (q - 1)*2)%sf(j, k, l) = nbub*r_IP(q)
274-
q_cons_vf(bubxb + (q - 1)*2 + 1)%sf(j, k, l) = nbub*v_IP(q)
275-
if (.not. polytropic) then
276-
q_cons_vf(bubxb + (q - 1)*4)%sf(j, k, l) = nbub*r_IP(q)
277-
q_cons_vf(bubxb + (q - 1)*4 + 1)%sf(j, k, l) = nbub*v_IP(q)
278-
q_cons_vf(bubxb + (q - 1)*4 + 2)%sf(j, k, l) = nbub*pb_IP(q)
279-
q_cons_vf(bubxb + (q - 1)*4 + 3)%sf(j, k, l) = nbub*mv_IP(q)
280-
end if
209+
do q = 1, num_fluids
210+
q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
211+
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
281212
end do
282-
end if
283213

284-
if (qbmm) then
214+
if (surface_tension) then
215+
q_prim_vf(c_idx)%sf(j, k, l) = c_IP
216+
end if
217+
if (model_eqns /= 4) then
218+
! If in simulation, use acc mixture subroutines
219+
if (elasticity) then
220+
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
221+
alpha_rho_IP, Re_K, G_K, Gs)
222+
else if (bubbles_euler) then
223+
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
224+
alpha_rho_IP, Re_K)
225+
else
226+
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
227+
alpha_rho_IP, Re_K)
228+
end if
229+
end if
285230

286-
nbub = nmom_IP(1)
231+
! Calculate velocity of ghost cell
232+
if (gp%slip) then
233+
norm(1:3) = levelset_norm%sf(gp%loc(1), gp%loc(2), gp%loc(3), gp%ib_patch_id, 1:3)
234+
buf = sqrt(sum(norm**2))
235+
norm = norm/buf
236+
vel_norm_IP = sum(vel_IP*norm)*norm
237+
vel_g = vel_IP - vel_norm_IP
238+
else
239+
vel_g = 0._wp
240+
end if
241+
242+
! Set momentum
287243
$:GPU_LOOP(parallelism='[seq]')
288-
do q = 1, nb*nmom
289-
q_cons_vf(bubxb + q - 1)%sf(j, k, l) = nbub*nmom_IP(q)
244+
do q = momxb, momxe
245+
q_cons_vf(q)%sf(j, k, l) = rho*vel_g(q - momxb + 1)
246+
dyn_pres = dyn_pres + q_cons_vf(q)%sf(j, k, l)* &
247+
vel_g(q - momxb + 1)/2._wp
290248
end do
291249

250+
! Set continuity and adv vars
292251
$:GPU_LOOP(parallelism='[seq]')
293-
do q = 1, nb
294-
q_cons_vf(bubxb + (q - 1)*nmom)%sf(j, k, l) = nbub
252+
do q = 1, num_fluids
253+
q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
254+
q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
295255
end do
296256

297-
if (.not. polytropic) then
257+
! Set color function
258+
if (surface_tension) then
259+
q_cons_vf(c_idx)%sf(j, k, l) = c_IP
260+
end if
261+
262+
! Set Energy
263+
if (bubbles_euler) then
264+
q_cons_vf(E_idx)%sf(j, k, l) = (1 - alpha_IP(1))*(gamma*pres_IP + pi_inf + dyn_pres)
265+
else
266+
q_cons_vf(E_idx)%sf(j, k, l) = gamma*pres_IP + pi_inf + dyn_pres
267+
end if
268+
! Set bubble vars
269+
if (bubbles_euler .and. .not. qbmm) then
270+
call s_comp_n_from_prim(alpha_IP(1), r_IP, nbub, weight)
271+
$:GPU_LOOP(parallelism='[seq]')
272+
do q = 1, nb
273+
q_cons_vf(bubxb + (q - 1)*2)%sf(j, k, l) = nbub*r_IP(q)
274+
q_cons_vf(bubxb + (q - 1)*2 + 1)%sf(j, k, l) = nbub*v_IP(q)
275+
if (.not. polytropic) then
276+
q_cons_vf(bubxb + (q - 1)*4)%sf(j, k, l) = nbub*r_IP(q)
277+
q_cons_vf(bubxb + (q - 1)*4 + 1)%sf(j, k, l) = nbub*v_IP(q)
278+
q_cons_vf(bubxb + (q - 1)*4 + 2)%sf(j, k, l) = nbub*pb_IP(q)
279+
q_cons_vf(bubxb + (q - 1)*4 + 3)%sf(j, k, l) = nbub*mv_IP(q)
280+
end if
281+
end do
282+
end if
283+
284+
if (qbmm) then
285+
286+
nbub = nmom_IP(1)
287+
$:GPU_LOOP(parallelism='[seq]')
288+
do q = 1, nb*nmom
289+
q_cons_vf(bubxb + q - 1)%sf(j, k, l) = nbub*nmom_IP(q)
290+
end do
291+
298292
$:GPU_LOOP(parallelism='[seq]')
299293
do q = 1, nb
294+
q_cons_vf(bubxb + (q - 1)*nmom)%sf(j, k, l) = nbub
295+
end do
296+
297+
if (.not. polytropic) then
300298
$:GPU_LOOP(parallelism='[seq]')
301-
do r = 1, nnode
302-
pb_in(j, k, l, r, q) = presb_IP((q - 1)*nnode + r)
303-
mv_in(j, k, l, r, q) = massv_IP((q - 1)*nnode + r)
299+
do q = 1, nb
300+
$:GPU_LOOP(parallelism='[seq]')
301+
do r = 1, nnode
302+
pb_in(j, k, l, r, q) = presb_IP((q - 1)*nnode + r)
303+
mv_in(j, k, l, r, q) = massv_IP((q - 1)*nnode + r)
304+
end do
304305
end do
305-
end do
306+
end if
306307
end if
307-
end if
308308

309-
if (model_eqns == 3) then
310-
$:GPU_LOOP(parallelism='[seq]')
311-
do q = intxb, intxe
312-
q_cons_vf(q)%sf(j, k, l) = alpha_IP(q - intxb + 1)*(gammas(q - intxb + 1)*pres_IP &
313-
+ pi_infs(q - intxb + 1))
314-
end do
315-
end if
316-
end do
317-
#:endcall GPU_PARALLEL_LOOP
309+
if (model_eqns == 3) then
310+
$:GPU_LOOP(parallelism='[seq]')
311+
do q = intxb, intxe
312+
q_cons_vf(q)%sf(j, k, l) = alpha_IP(q - intxb + 1)*(gammas(q - intxb + 1)*pres_IP &
313+
+ pi_infs(q - intxb + 1))
314+
end do
315+
end if
316+
end do
317+
#:endcall GPU_PARALLEL_LOOP
318+
end if
318319

319320
!Correct the state of the inner points in IBs
320321
if (num_inner_gps > 0) then

0 commit comments

Comments
 (0)