@@ -182,7 +182,7 @@ contains
182182
183183 coeffs = 0d0
184184
185- do i1 = 0 , 2 ; do i2 = 0 , 2
185+ do i2 = 0 , 2 ; do i1 = 0 , 2
186186 if ((i1 + i2) <= 2 ) then
187187 if (bubble_model == 3 ) then
188188 ! RPE
@@ -223,7 +223,8 @@ contains
223223 real (kind (0d0 )), dimension (nb) :: Rvec
224224 real (kind (0d0 )), dimension (nnode, nb) :: wght, abscX, abscY
225225 real (kind (0d0 )), dimension (nterms, 0 :2 , 0 :2 ) :: mom3d_terms, coeff
226- real (kind (0d0 )) :: pres, rho, nbub, c, alf
226+ real (kind (0d0 )) :: pres, rho, nbub, c, alf, R3, momsum
227+ real (kind (0d0 )) :: start, finish
227228 real (kind (0d0 )) :: n_tait, B_tait
228229
229230 integer :: j, k, l, q, r, s !< Loop variables
@@ -234,7 +235,8 @@ contains
234235
235236 !$acc update device(is1, is2, is3)
236237
237- !$acc parallel loop collapse(3 ) gang vector default(present) private(moms, Rvec, wght, abscX, abscY, mom3d_terms, coeff)
238+
239+ !$acc parallel loop collapse(3 ) gang vector default(present) private(moms, wght, abscX, abscY, coeff)
238240 do id3 = is3%beg, is3%end
239241 do id2 = is2%beg, is2%end
240242 do id1 = is1%beg, is1%end
@@ -261,12 +263,14 @@ contains
261263
262264 if (alf > small_alf) then
263265
266+ R3 = 0d0
267+
264268 !$acc loop seq
265269 do q = 1 , nb
266- Rvec(q) = q_prim_vf(bubrs(q))%sf(id1, id2, id3)
270+ R3 = R3 + weight(q) * q_prim_vf(bubrs(q))%sf(id1, id2, id3)** 3d0
267271 end do
268272
269- call s_comp_n_from_prim(alf, Rvec, nbub)
273+ nbub = ( 3.d0 / ( 4.d0 * pi)) * alf / R3
270274
271275 !$acc loop seq
272276 do q = 1 , nb
@@ -275,48 +279,29 @@ contains
275279 moms(r) = q_prim_vf(bubmoms(q, r))%sf(id1, id2, id3)
276280 end do
277281
278- ! IF (id1==0 ) THEN
279- ! PRINT * , ' pres: ' , pres
280- ! PRINT * , ' nb : ' , nbub
281- ! PRINT * , ' alf: ' , alf
282- ! DO s = 1 ,nmom
283- ! PRINT * , ' mom: ' , moms(s)
284- ! END DO
285- ! END IF
282+
286283
287284 call s_chyqmom(moms, wght(:, q), abscX(:, q), abscY(:, q))
288285
289- !$acc loop seq
290- do j = 1 , nterms
291- !$acc loop seq
292- do i2 = 0 , 2
293- !$acc loop seq
294- do i1 = 0 , 2
295- if ((i1 + i2) <= 2 ) then
296-
297- mom3d_terms(j, i1, i2) = coeff(j, i1, i2)* (R0(q)** momrhs(3 , i1, i2, j, q)) &
298- * f_quad2D(abscX(:, q), abscY(:, q), wght(:, q), momrhs(:, i1, i2, j, q))
299- end if
300- end do
301- end do
302- end do
303286
304287 !$acc loop seq
305- do i1 = 0 , 2
288+ do i2 = 0 , 2
306289 !$acc loop seq
307- do i2 = 0 , 2
290+ do i1 = 0 , 2
308291 if ((i1 + i2) <= 2 ) then
309- moms3d(i1, i2, q)%sf(id1, id2, id3) = nbub * sum (mom3d_terms(:, i1, i2))
310- ! IF (moms3d(i1,i2,q)%sf(id1,id2,id3) .NE. moms3d(i1,i2,q)%sf(id1,id2,id3)) THEN
311- ! PRINT * , ' nan in mom3d ' , i1,i2,id1
312- ! PRINT * , ' nbu: ' , nbub
313- ! PRINT * , ' alf: ' , alf
314- ! PRINT * , ' moms: ' , moms(:)
315- ! CALL s_mpi_abort()
316- ! END IF
292+ momsum = 0d0
293+ !$acc loop seq
294+ do j = 1 , nterms
295+ momsum = momsum + coeff(j, i1, i2) * (R0(q) ** momrhs( 3 , i1, i2, j, q)) &
296+ * f_quad2D(abscX(:, q), abscY(:, q), wght(:, q), momrhs(:, i1, i2, j, q))
297+ end do
298+ moms3d(i1, i2, q)%sf(id1, id2, id3) = nbub * momsum
299+
317300 end if
318301 end do
319302 end do
303+
304+
320305 end do
321306
322307 momsp(1 )%sf(id1, id2, id3) = f_quad(abscX, abscY, wght, 3d0 , 0d0 , 0d0 )
@@ -329,19 +314,7 @@ contains
329314 momsp(4)%sf(id1, id2, id3) = f_quad(abscX, abscY, wght, 3d0*(1d0 - gam), 0d0, 3d0*gam)
330315 end if
331316
332- !!$acc loop seq
333- !do i1 = 1, 4
334- ! if (momsp(i1)%sf(id1, id2, id3) /= momsp(i1)%sf(id1, id2, id3)) then
335- ! print *, ' NaN in sp moment' , i1, ' location' , id1, id2, id3
336- ! print *, ' Rs' , Rvec(:)
337- ! print *, ' alpha' , alf
338- ! print *, ' nbub' , nbub
339- ! print *, ' abscX' , abscX(:, :)
340- ! print *, ' abscY' , abscY(:, :)
341- ! print *, ' wght' , wght(:, :)
342- ! call s_mpi_abort()
343- !end if
344- !end do
317+
345318 else
346319 !$acc loop seq
347320 do q = 1, nb
@@ -365,6 +338,7 @@ contains
365338 end do
366339 end do
367340
341+
368342 end subroutine s_mom_inv
369343
370344 subroutine s_chyqmom(momin, wght, abscX, abscY)
0 commit comments