Skip to content

Commit 372ddf0

Browse files
danieljvickersdanieljvickers
committed
Fixed compiler issue with airfoil changes. All tests pass locally
1 parent 0cabfff commit 372ddf0

File tree

2 files changed

+35
-37
lines changed

2 files changed

+35
-37
lines changed

‎src/common/m_compute_levelset.fpp‎

Lines changed: 31 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -65,36 +65,35 @@ contains
6565

6666
pure subroutine s_airfoil_levelset(ib_patch_id, levelset, levelset_norm)
6767

68-
type(levelset_field), intent(INOUT), optional :: levelset
69-
type(levelset_norm_field), intent(INOUT), optional :: levelset_norm
70-
integer, intent(IN) :: ib_patch_id
68+
type(levelset_field), intent(inout), optional :: levelset
69+
type(levelset_norm_field), intent(inout), optional :: levelset_norm
70+
integer, intent(in) :: ib_patch_id
7171

7272
real(wp) :: dist, global_dist
7373
integer :: global_id
74-
real(wp) :: x_centroid, y_centroid, x_act, y_act, theta
74+
real(wp) :: x_centroid, y_centroid
7575
real(wp), dimension(3) :: dist_vec
7676

77+
real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
78+
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
79+
7780
integer :: i, j, k !< Loop index variables
7881

7982
x_centroid = patch_ib(ib_patch_id)%x_centroid
8083
y_centroid = patch_ib(ib_patch_id)%y_centroid
81-
theta = pi*patch_ib(ib_patch_id)%theta/180._wp
84+
inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
85+
rotation(:, :) = patch_ib(patch_id)%rotation_matrix(:, :)
8286

8387
do i = 0, m
8488
do j = 0, n
89+
xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp] ! get coordinate frame centered on IB
90+
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordiantes
8591
86-
if (.not. f_is_default(patch_ib(ib_patch_id)%theta)) then
87-
x_act = (x_cc(i) - x_centroid)*cos(theta) - (y_cc(j) - y_centroid)*sin(theta) + x_centroid
88-
y_act = (x_cc(i) - x_centroid)*sin(theta) + (y_cc(j) - y_centroid)*cos(theta) + y_centroid
89-
else
90-
x_act = x_cc(i)
91-
y_act = y_cc(j)
92-
end if
93-
94-
if (y_act >= y_centroid) then
92+
if (xy_local(2) >= 0._wp) then
9593
do k = 1, Np
96-
dist_vec(1) = x_cc(i) - airfoil_grid_u(k)%x
97-
dist_vec(2) = y_cc(j) - airfoil_grid_u(k)%y
94+
! TODO :: I think that witht he IBM implementation that airfoil_grid_u(k).y = 0._wp for all indices. Delete?
95+
dist_vec(1) = xy_local(1) - airfoil_grid_u(k)%x
96+
dist_vec(2) = xy_local(2) - airfoil_grid_u(k)%y
9897
dist_vec(3) = 0
9998
dist = sqrt(sum(dist_vec**2))
10099
if (k == 1) then
@@ -107,14 +106,14 @@ contains
107106
end if
108107
end if
109108
end do
110-
dist_vec(1) = x_cc(i) - airfoil_grid_u(global_id)%x
111-
dist_vec(2) = y_cc(j) - airfoil_grid_u(global_id)%y
109+
dist_vec(1) = xy_local(1) - airfoil_grid_u(global_id)%x
110+
dist_vec(2) = xy_local(2) - airfoil_grid_u(global_id)%y
112111
dist_vec(3) = 0
113112
dist = global_dist
114113
else
115114
do k = 1, Np
116-
dist_vec(1) = x_cc(i) - airfoil_grid_l(k)%x
117-
dist_vec(2) = y_cc(j) - airfoil_grid_l(k)%y
115+
dist_vec(1) = xy_local(1) - airfoil_grid_l(k)%x
116+
dist_vec(2) = xy_local(2) - airfoil_grid_l(k)%y
118117
dist_vec(3) = 0
119118
dist = sqrt(sum(dist_vec**2))
120119
if (k == 1) then
@@ -127,8 +126,8 @@ contains
127126
end if
128127
end if
129128
end do
130-
dist_vec(1) = x_cc(i) - airfoil_grid_l(global_id)%x
131-
dist_vec(2) = y_cc(j) - airfoil_grid_l(global_id)%y
129+
dist_vec(1) = xy_local(1) - airfoil_grid_l(global_id)%x
130+
dist_vec(2) = xy_local(2) - airfoil_grid_l(global_id)%y
132131
dist_vec(3) = 0
133132
dist = global_dist
134133
end if
@@ -138,7 +137,7 @@ contains
138137
levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
139138
else
140139
levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
141-
dist_vec(:)/dist
140+
matmul(rotation, dist_vec(:)/dist) ! convert the normal vector back to global grid coordinates
142141
end if
143142
144143
end do
@@ -154,9 +153,12 @@ contains
154153
155154
real(wp) :: dist, dist_surf, dist_side, global_dist
156155
integer :: global_id
157-
real(wp) :: x_centroid, y_centroid, z_centroid, lz, z_max, z_min, x_act, y_act, theta
156+
real(wp) :: x_centroid, y_centroid, z_centroid, lz, z_max, z_min
158157
real(wp), dimension(3) :: dist_vec
159158
159+
real(wp), dimension(1:3) :: xyz_local !< x, y, z coordinates in local IB frame
160+
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
161+
160162
real(wp) :: length_z
161163
162164
integer :: i, j, k, l !< Loop index variables
@@ -165,7 +167,8 @@ contains
165167
y_centroid = patch_ib(ib_patch_id)%y_centroid
166168
z_centroid = patch_ib(ib_patch_id)%z_centroid
167169
lz = patch_ib(ib_patch_id)%length_z
168-
theta = pi*patch_ib(ib_patch_id)%theta/180._wp
170+
inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
171+
rotation(:, :) = patch_ib(patch_id)%rotation_matrix(:, :)
169172
170173
z_max = z_centroid + lz/2
171174
z_min = z_centroid - lz/2
@@ -174,15 +177,10 @@ contains
174177
do j = 0, n
175178
do i = 0, m
176179
177-
if (.not. f_is_default(patch_ib(ib_patch_id)%theta)) then
178-
x_act = (x_cc(i) - x_centroid)*cos(theta) - (y_cc(j) - y_centroid)*sin(theta) + x_centroid
179-
y_act = (x_cc(i) - x_centroid)*sin(theta) + (y_cc(j) - y_centroid)*cos(theta) + y_centroid
180-
else
181-
x_act = x_cc(i)
182-
y_act = y_cc(j)
183-
end if
180+
xyz_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(l) - z_centroid] ! get coordinate frame centered on IB
181+
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordiantes
184182

185-
if (y_act >= y_centroid) then
183+
if (xzy_local(2) >= y_centroid) then
186184
do k = 1, Np
187185
dist_vec(1) = x_cc(i) - airfoil_grid_u(k)%x
188186
dist_vec(2) = y_cc(j) - airfoil_grid_u(k)%y

‎src/common/m_ib_patches.fpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -180,7 +180,7 @@ contains
180180
integer, intent(in) :: patch_id
181181
integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
182182

183-
real(wp) :: f, ca_in, pa, ma, ta, theta
183+
real(wp) :: f, ca_in, pa, ma, ta
184184
real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
185185
integer :: i, j, k
186186
integer :: Np1, Np2
@@ -194,7 +194,6 @@ contains
194194
pa = patch_ib(patch_id)%p
195195
ma = patch_ib(patch_id)%m
196196
ta = patch_ib(patch_id)%t
197-
theta = pi*patch_ib(patch_id)%theta/180._wp
198197
inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
199198

200199
! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors
@@ -340,6 +339,7 @@ contains
340339
pa = patch_ib(patch_id)%p
341340
ma = patch_ib(patch_id)%m
342341
ta = patch_ib(patch_id)%t
342+
inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
343343
344344
! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors
345345
#ifdef MFC_PRE_PROCESS
@@ -653,7 +653,7 @@ contains
653653
cart_z = z_cc(k)
654654
end if
655655
xyz_local = [x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid] ! get coordinate frame centered on IB
656-
xyz_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordiantes
656+
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordiantes
657657
658658
if (x_boundary%beg <= xyz_local(1) .and. &
659659
x_boundary%end >= xyz_local(1) .and. &
@@ -733,7 +733,7 @@ contains
733733
cart_z = z_cc(k)
734734
end if
735735
xyz_local = [x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid] ! get coordinate frame centered on IB
736-
xyz_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordiantes
736+
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordiantes
737737
738738
if (((.not. f_is_default(length_x) .and. &
739739
xyz_local(2)**2 &

0 commit comments

Comments
 (0)