Skip to content

Commit bd6096b

Browse files
danieljvickersdanieljvickers
committed
Passed test suite of IBM cases
1 parent 77c1432 commit bd6096b

File tree

2 files changed

+43
-58
lines changed

2 files changed

+43
-58
lines changed

‎src/common/m_compute_levelset.fpp‎

Lines changed: 28 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -187,7 +187,7 @@ contains
187187
xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
188188
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
189189

190-
if (xyz_local(2) >= y_centroid) then
190+
if (xyz_local(2) >= center(2)) then
191191
do k = 1, Np
192192
dist_vec(1) = xyz_local(1) - airfoil_grid_u(k)%x
193193
dist_vec(2) = xyz_local(2) - airfoil_grid_u(k)%y
@@ -268,7 +268,7 @@ contains
268268
real(wp) :: side_dists(4)
269269

270270
real(wp) :: length_x, length_y
271-
real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
271+
real(wp), dimension(1:3) :: xy_local, dist_vec !< x and y coordinates in local IB frame
272272
real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
273273
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
274274

@@ -287,7 +287,7 @@ contains
287287
bottom_left(1) = -length_x/2
288288
bottom_left(2) = -length_y/2
289289

290-
$:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local]', copy='[levelset,levelset_norm]',&
290+
$:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset,levelset_norm]',&
291291
& copyin='[ib_patch_id,center,bottom_left,top_right,initial_distance_buffer,inverse_rotation,rotation]', collapse=2)
292292
do i = 0, m
293293
do j = 0, n
@@ -312,20 +312,17 @@ contains
312312
end do
313313

314314
levelset%sf(i, j, 0, ib_patch_id) = side_dists(idx)
315-
levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp
315+
dist_vec = 0._wp
316316
if (.not. f_approx_equal(side_dists(idx), 0._wp)) then
317317
if (idx == 1 .or. idx == 2) then
318318
! vector points along the x axis
319-
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(idx)/ &
320-
abs(side_dists(idx))
319+
dist_vec(1) = side_dists(idx)/abs(side_dists(idx))
321320
else
322321
! vector points along the y axis
323-
levelset_norm%sf(i, j, 0, ib_patch_id, 2) = side_dists(idx)/ &
324-
abs(side_dists(idx))
322+
dist_vec(2) = side_dists(idx)/abs(side_dists(idx))
325323
end if
326324
! convert the normal vector back into the global coordinate system
327-
levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
328-
matmul(rotation, levelset_norm%sf(i, j, 0, ib_patch_id, :))
325+
levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, dist_vec)
329326
else
330327
levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp
331328
end if
@@ -347,7 +344,7 @@ contains
347344

348345
real(wp), dimension(3) :: center
349346
real(wp) :: length_x, length_y, length_z
350-
real(wp), dimension(1:3) :: xyz_local !< x and y coordinates in local IB frame
347+
real(wp), dimension(1:3) :: xyz_local, dist_vec !< x and y coordinates in local IB frame
351348
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
352349

353350
integer :: i, j, k !< Loop index variables
@@ -370,7 +367,7 @@ contains
370367
Front = length_z/2
371368
Back = -length_z/2
372369

373-
$:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local]', copy='[levelset,levelset_norm]',&
370+
$:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local,dist_vec]', copy='[levelset,levelset_norm]',&
374371
& copyin='[ib_patch_id,center,inverse_rotation,rotation,Right,Left,Top,Bottom,Front,Back]', collapse=3)
375372
do i = 0, m
376373
do j = 0, n
@@ -398,51 +395,44 @@ contains
398395
! leading to undesired behavior. This should be resolved
399396
! and this code should be cleaned up. It also means that
400397
! rotating the box 90 degrees will cause tests to fail.
401-
levelset_norm%sf(i, j, k, ib_patch_id, :) = 0._wp
398+
dist_vec = 0._wp
402399
if (f_approx_equal(min_dist, abs(side_dists(1)))) then
403400
levelset%sf(i, j, k, ib_patch_id) = side_dists(1)
404401
if (.not. f_approx_equal(side_dists(1), 0._wp)) then
405-
levelset_norm%sf(i, j, k, ib_patch_id, 1) = side_dists(1)/ &
406-
abs(side_dists(1))
402+
dist_vec(1) = side_dists(1)/abs(side_dists(1))
407403
end if
408404
409405
else if (f_approx_equal(min_dist, abs(side_dists(2)))) then
410406
levelset%sf(i, j, k, ib_patch_id) = side_dists(2)
411407
if (.not. f_approx_equal(side_dists(2), 0._wp)) then
412-
levelset_norm%sf(i, j, k, ib_patch_id, 1) = -side_dists(2)/ &
413-
abs(side_dists(2))
408+
dist_vec(1) = -side_dists(2)/abs(side_dists(2))
414409
end if
415410
416411
else if (f_approx_equal(min_dist, abs(side_dists(3)))) then
417412
levelset%sf(i, j, k, ib_patch_id) = side_dists(3)
418413
if (.not. f_approx_equal(side_dists(3), 0._wp)) then
419-
levelset_norm%sf(i, j, k, ib_patch_id, 2) = side_dists(3)/ &
420-
abs(side_dists(3))
414+
dist_vec(2) = side_dists(3)/abs(side_dists(3))
421415
end if
422416
423417
else if (f_approx_equal(min_dist, abs(side_dists(4)))) then
424418
levelset%sf(i, j, k, ib_patch_id) = side_dists(4)
425419
if (.not. f_approx_equal(side_dists(4), 0._wp)) then
426-
levelset_norm%sf(i, j, k, ib_patch_id, 2) = -side_dists(4)/ &
427-
abs(side_dists(4))
420+
dist_vec(2) = -side_dists(4)/abs(side_dists(4))
428421
end if
429422
430423
else if (f_approx_equal(min_dist, abs(side_dists(5)))) then
431424
levelset%sf(i, j, k, ib_patch_id) = side_dists(5)
432425
if (.not. f_approx_equal(side_dists(5), 0._wp)) then
433-
levelset_norm%sf(i, j, k, ib_patch_id, 3) = side_dists(5)/ &
434-
abs(side_dists(5))
426+
dist_vec(3) = side_dists(5)/abs(side_dists(5))
435427
end if
436428
437429
else if (f_approx_equal(min_dist, abs(side_dists(6)))) then
438430
levelset%sf(i, j, k, ib_patch_id) = side_dists(6)
439431
if (.not. f_approx_equal(side_dists(6), 0._wp)) then
440-
levelset_norm%sf(i, j, k, ib_patch_id, 3) = -side_dists(6)/ &
441-
abs(side_dists(6))
432+
dist_vec(3) = -side_dists(6)/abs(side_dists(6))
442433
end if
443434
end if
444-
levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
445-
matmul(rotation, levelset_norm%sf(i, j, 0, ib_patch_id, :))
435+
levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, dist_vec)
446436
end if
447437
end do
448438
end do
@@ -479,8 +469,7 @@ contains
479469
if (f_approx_equal(dist, 0._wp)) then
480470
levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
481471
else
482-
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
483-
dist_vec(:)/dist
472+
levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_vec(:)/dist
484473
end if
485474
end do
486475
end do
@@ -495,7 +484,7 @@ contains
495484
integer, intent(IN) :: ib_patch_id
496485
497486
real(wp) :: radius
498-
real(wp), dimension(3) :: centroid_vec, dist_sides_vec, dist_surface_vec, length
487+
real(wp), dimension(3) :: dist_sides_vec, dist_surface_vec, length
499488
real(wp), dimension(2) :: boundary
500489
real(wp) :: dist_side, dist_surface, side_pos
501490
integer :: i, j, k !< Loop index variables
@@ -514,24 +503,24 @@ contains
514503
inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
515504
rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
516505
517-
if (.not. f_approx_equal(length_x, 0._wp)) then
506+
if (.not. f_approx_equal(length(1), 0._wp)) then
518507
boundary(1) = -0.5_wp*length(1)
519508
boundary(2) = 0.5_wp*length(1)
520509
dist_sides_vec = (/1, 0, 0/)
521510
dist_surface_vec = (/0, 1, 1/)
522-
else if (.not. f_approx_equal(length_y, 0._wp)) then
511+
else if (.not. f_approx_equal(length(2), 0._wp)) then
523512
boundary(1) = -0.5_wp*length(2)
524513
boundary(2) = 0.5_wp*length(2)
525514
dist_sides_vec = (/0, 1, 0/)
526515
dist_surface_vec = (/1, 0, 1/)
527-
else if (.not. f_approx_equal(length_z, 0._wp)) then
516+
else if (.not. f_approx_equal(length(3), 0._wp)) then
528517
boundary(1) = -0.5_wp*length(3)
529518
boundary(2) = 0.5_wp*length(3)
530519
dist_sides_vec = (/0, 0, 1/)
531520
dist_surface_vec = (/1, 1, 0/)
532521
end if
533522
534-
$:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset,levelset_norm]',&
523+
$:GPU_PARALLEL_LOOP(private='[i,j,k,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset,levelset_norm]',&
535524
& copyin='[ib_patch_id,center,radius,inverse_rotation,rotation,dist_sides_vec,dist_surface_vec]', collapse=3)
536525
do i = 0, m
537526
do j = 0, n
@@ -551,21 +540,17 @@ contains
551540
! if the closest edge is flat
552541
levelset%sf(i, j, k, ib_patch_id) = -dist_side
553542
if (f_approx_equal(dist_side, abs(side_pos - boundary(1)))) then
554-
levelset_norm%sf(i, j, k, ib_patch_id, :) = -dist_sides_vec
543+
levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, -dist_sides_vec)
555544
else
556-
levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec
545+
levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, dist_sides_vec)
557546
end if
558547
else
559548
levelset%sf(i, j, k, ib_patch_id) = dist_surface
560549

561-
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
562-
xyz_local*dist_surface_vec
563-
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
564-
levelset_norm%sf(i, j, k, ib_patch_id, :)/ &
565-
norm2(levelset_norm%sf(i, j, k, ib_patch_id, :))
550+
xyz_local = xyz_local*dist_surface_vec
551+
xyz_local = xyz_local/norm2(xyz_local)
552+
levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, xyz_local)
566553
end if
567-
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
568-
matmul(rotation, levelset_norm%sf(i, j, k, ib_patch_id, :))
569554
end do
570555
end do
571556
end do

‎src/common/m_ib_patches.fpp‎

Lines changed: 15 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -565,11 +565,11 @@ contains
565565
566566
integer, intent(in) :: patch_id
567567
integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
568-
integer, dimension(1:3):: center
569568
570569
! Generic loop iterators
571570
integer :: i, j, k
572571
real(wp) :: radius
572+
real(wp), dimension(1:3):: center
573573
574574
!! Variables to initialize the pressure field that corresponds to the
575575
!! bubble-collapse test case found in Tiwari et al. (2013)
@@ -590,7 +590,7 @@ contains
590590
! and verifying whether the current patch has permission to write to
591591
! that cell. If both queries check out, the primitive variables of
592592
! the current patch are assigned to this cell.
593-
$:GPU_PARALLEL_LOOP(private='[i,j,k]', copy='[ib_markers_sf]',&
593+
$:GPU_PARALLEL_LOOP(private='[i,j,k,cart_y,cart_z]', copy='[ib_markers_sf]',&
594594
& copyin='[patch_id,center,radius]', collapse=3)
595595
do k = 0, p
596596
do j = 0, n
@@ -651,7 +651,7 @@ contains
651651
! and verifying whether the current patch has permission to write to
652652
! to that cell. If both queries check out, the primitive variables
653653
! of the current patch are assigned to this cell.
654-
$:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local]', copy='[ib_markers_sf]',&
654+
$:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local,cart_y,cart_z]', copy='[ib_markers_sf]',&
655655
& copyin='[patch_id,center,length,inverse_rotation]', collapse=3)
656656
do k = 0, p
657657
do j = 0, n
@@ -664,7 +664,7 @@ contains
664664
cart_y = y_cc(j)
665665
cart_z = z_cc(k)
666666
end if
667-
xyz_local = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! get coordinate frame centered on IB
667+
xyz_local = [x_cc(i), cart_y, cart_z] - center ! get coordinate frame centered on IB
668668
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
669669
670670
if (-0.5*length(1) <= xyz_local(1) .and. &
@@ -725,7 +725,7 @@ contains
725725
! domain and verifying whether the current patch has the permission
726726
! to write to that cell. If both queries check out, the primitive
727727
! variables of the current patch are assigned to this cell.
728-
$:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local]', copy='[ib_markers_sf]',&
728+
$:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local,cart_y,cart_z]', copy='[ib_markers_sf]',&
729729
& copyin='[patch_id,center,length,radius,inverse_rotation]', collapse=3)
730730
do k = 0, p
731731
do j = 0, n
@@ -737,26 +737,26 @@ contains
737737
cart_y = y_cc(j)
738738
cart_z = z_cc(k)
739739
end if
740-
xyz_local = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! get coordinate frame centered on IB
740+
xyz_local = [x_cc(i), cart_y, cart_z] - center ! get coordinate frame centered on IB
741741
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
742742
743-
if (((.not. f_is_default(length_x) .and. &
743+
if (((.not. f_is_default(length(1)) .and. &
744744
xyz_local(2)**2 &
745745
+ xyz_local(3)**2 <= radius**2 .and. &
746-
-0.5*length(1) <= xyz_local(1) .and. &
747-
0.5*length(1) >= xyz_local(1)) &
746+
-0.5_wp*length(1) <= xyz_local(1) .and. &
747+
0.5_wp*length(1) >= xyz_local(1)) &
748748
.or. &
749-
(.not. f_is_default(length_y) .and. &
749+
(.not. f_is_default(length(2)) .and. &
750750
xyz_local(1)**2 &
751751
+ xyz_local(3)**2 <= radius**2 .and. &
752-
-0.5*length(2) <= xyz_local(2) .and. &
753-
0.5*length(2) >= xyz_local(2)) &
752+
-0.5_wp*length(2) <= xyz_local(2) .and. &
753+
0.5_wp*length(2) >= xyz_local(2)) &
754754
.or. &
755-
(.not. f_is_default(length_z) .and. &
755+
(.not. f_is_default(length(3)) .and. &
756756
xyz_local(1)**2 &
757757
+ xyz_local(2)**2 <= radius**2 .and. &
758-
-0.5*length(3) <= xyz_local(3) .and. &
759-
0.5*length(3) >= xyz_local(3)))) then
758+
-0.5_wp*length(3) <= xyz_local(3) .and. &
759+
0.5_wp*length(3) >= xyz_local(3)))) then
760760
761761
! Updating the patch identities bookkeeping variable
762762
ib_markers_sf(i, j, k) = patch_id

0 commit comments

Comments
 (0)