Skip to content

Commit

Permalink
Merge pull request #289 from xcompact3d/fix_intel
Browse files Browse the repository at this point in the history 
Fix intel compiler bugs (caused by local 3D arrays)
  • Loading branch information
@slaizet
slaizet committed Sep 5, 2024
2 parents 5565268 + 1ac666d commit 46e874e
Show file tree
Hide file tree
Showing 12 changed files with 309 additions and 186 deletions.
83 changes: 49 additions & 34 deletions src/Case-ABL.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -675,21 +675,36 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1)
use param
use variables
use var, only: di1, di2, di3
use var, only: sxy1, syz1, tb1, ta2, tb2
use var, only: sxy1, syz1
use var, only: te1, th1, ti1
use var, only: ta2, tb2, tc2, td2, te2, th2, ti2
use var, only: th3, ti3
use ibm_param, only : ubcx, ubcy, ubcz

implicit none
real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1,uy1,uz1,nut1
real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(out) :: wallsgsx1,wallsgsy1,wallsgsz1

! Local variables
real(mytype),dimension(ysize(1),ysize(3)) :: tauwallxy2, tauwallzy2
integer :: i,j,k,code,j0
integer :: nxc, nyc, nzc, xsize1, xsize2, xsize3
real(mytype) :: delta
real(mytype) :: ux_HAve_local, uz_HAve_local
real(mytype) :: ux_HAve, uz_HAve, S_HAve, ux_delta, uz_delta, S_delta
real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: txy1,tyz1,dtwxydx
real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: txy2,tyz2,wallsgsx2,wallsgsz2
real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: tyz3,dtwyzdz
!real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: txy1,tyz1,dtwxydx
!real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: txy2,tyz2,wallsgsx2,wallsgsz2
!real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: tyz3,dtwyzdz

! tc -> ux
! td -> uz
! txy -> te
! tyz -> th
! wallsgsx -> ta
! wallsgsz -> tb
! dtwxydx -> ti1
! dtwyzdz -> ti3
! tauwallxy2 -> sy

real(mytype) :: y_sampling

Expand All @@ -703,47 +718,47 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1)

if(iconserv==0) then
! Construct Smag SGS stress tensor
txy1 = -2.0*nut1*sxy1
tyz1 = -2.0*nut1*syz1
call transpose_x_to_y(txy1,txy2)
call transpose_x_to_y(tyz1,tyz2)
te1 = -2.0*nut1*sxy1
th1 = -2.0*nut1*syz1
call transpose_x_to_y(te1,te2)
call transpose_x_to_y(th1,th2)
endif

! Work on Y-pencil
call transpose_x_to_y(ux1,ta2)
call transpose_x_to_y(uz1,tb2)
call transpose_x_to_y(ux1,tc2)
call transpose_x_to_y(uz1,td2)

! Apply BCs locally
do k=1,ysize(3)
do i=1,ysize(1)
!sampling at dsampling*dy from wall
j0=floor(delta/dy)
y_sampling=delta-real(j0,mytype)*dy
ux_delta=(1-y_sampling/dy)*ta2(i,j0+1,k)+(y_sampling/dy)*ta2(i,j0+2,k)
uz_delta=(1-y_sampling/dy)*tb2(i,j0+1,k)+(y_sampling/dy)*tb2(i,j0+2,k)
S_delta=sqrt(ux_delta**2.+uz_delta**2.)
tauwallxy2(i,k)=-(k_roughness/(log(delta/z_zero)))**two*ux_delta*S_delta
tauwallzy2(i,k)=-(k_roughness/(log(delta/z_zero)))**two*uz_delta*S_delta
txy2(i,2,k) = tauwallxy2(i,k)
tyz2(i,2,k) = tauwallzy2(i,k)
enddo
do i=1,ysize(1)
!sampling at dsampling*dy from wall
j0=floor(delta/dy)
y_sampling=delta-real(j0,mytype)*dy
ux_delta=(1-y_sampling/dy)*tc2(i,j0+1,k)+(y_sampling/dy)*tc2(i,j0+2,k)
uz_delta=(1-y_sampling/dy)*td2(i,j0+1,k)+(y_sampling/dy)*td2(i,j0+2,k)
S_delta=sqrt(ux_delta**2.+uz_delta**2.)
tauwallxy2(i,k)=-(k_roughness/(log(delta/z_zero)))**two*ux_delta*S_delta
tauwallzy2(i,k)=-(k_roughness/(log(delta/z_zero)))**two*uz_delta*S_delta
te2(i,2,k) = tauwallxy2(i,k)
th2(i,2,k) = tauwallzy2(i,k)
enddo
enddo

if (iconserv==0) then
! Derivative of wallmodel-corrected SGS stress tensor
call dery_22(wallsgsx2,txy2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),1,ubcy)
call dery_22(wallsgsz2,tyz2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),1,ubcy)
call transpose_y_to_x(wallsgsx2,wallsgsx1)
call transpose_y_to_x(wallsgsz2,wallsgsz1)
call derx(dtwxydx,txy1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx)
call transpose_y_to_z(tyz2,tyz3)
call derz(dtwyzdz,tyz3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,ubcz)
call transpose_z_to_y(dtwyzdz,tb2)
call transpose_y_to_x(tb2,tb1)
wallsgsy1 = dtwxydx + tb1
call dery_22(ta2,te2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),1,ubcy)
call dery_22(tb2,th2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),1,ubcy)
call transpose_y_to_x(ta2,wallsgsx1)
call transpose_y_to_x(tb2,wallsgsz1)
call derx(wallsgsy1,te1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx)
call transpose_y_to_z(th2,th3)
call derz(ti1,th3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,ubcz)
call transpose_z_to_y(ti3,tb2)
call transpose_y_to_x(tb2,ti1)
wallsgsy1 = wallsgsy1 + ti1
elseif (iconserv==1) then
call transpose_y_to_x(txy2,wallsgsx1)
call transpose_y_to_x(tyz2,wallsgsz1)
call transpose_y_to_x(te2,wallsgsx1)
call transpose_y_to_x(th2,wallsgsz1)
endif


Expand Down
27 changes: 14 additions & 13 deletions src/forces.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,13 +121,13 @@ subroutine init_forces
elseif (i2dsim==0) then
write(*,*) '========================Forces============================='
write(*,*) ' (icvlf) (icvrt) (kcvbk) (kcvfr)'
write(*,*) ' (jcvup) B____________C B`_____________B'
write(*,*) ' | | | | |'
write(*,*) ' | __ | | | ____ |'
write(*,*) ' | \__\ | | | \___\ |'
write(*,*) ' | | | | |'
write(*,*) ' | CV | | | (Front) |'
write(*,*) ' (jcvlw) A____________D | A`_____________A'
write(*,*) ' (jcvup) B____________C B`_____________B '
write(*,*) ' \ \ | \ \ '
write(*,*) ' \ __ \ | \ ____ \ '
write(*,*) ' \ \__\ \ | \ \___\ \ '
write(*,*) ' \ \ | \ \ '
write(*,*) ' \ CV \ | \ (Front) \ '
write(*,*) ' (jcvlw) A____________D | A`_____________A '
do iv=1,nvol
write(*,"(' Control Volume : #',I1)") iv
write(*,"(' xld, icvlf : (',F6.2,',',I6,')')") xld(iv), icvlf(iv)
Expand Down Expand Up @@ -244,6 +244,7 @@ subroutine force(ux1,uy1,uz1,ep1)

use var, only : ta1, tb1, tc1, td1, te1, di1, tg1, tg2, tg3, th1, th2, th3, tf2, tf1
use var, only : ux2, ux3, uy2, uy3, uz2, ta2, tb2, td2, te2, di2, di3
use var, only : tc2

implicit none
character(len=30) :: filename, filename2
Expand All @@ -255,7 +256,7 @@ subroutine force(ux1,uy1,uz1,ep1)
real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1, uy1, uz1
real(mytype), dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ep1

real(mytype), dimension(ysize(1),ysize(2),ysize(3)) :: ppi2
!real(mytype), dimension(ysize(1),ysize(2),ysize(3)) :: ppi2 ! we'll use tc2

real(mytype), dimension(nz) :: yLift,xDrag
real(mytype) :: yLift_mean,xDrag_mean
Expand Down Expand Up @@ -315,7 +316,7 @@ subroutine force(ux1,uy1,uz1,ep1)

call transpose_x_to_y(ux1,ux2)
call transpose_x_to_y(uy1,uy2)
call transpose_x_to_y(ppi1,ppi2)
call transpose_x_to_y(ppi1,tc2)

call dery (td2,ux2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) ! dudy
call dery (te2,uy2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy) ! dvdy
Expand Down Expand Up @@ -506,7 +507,7 @@ subroutine force(ux1,uy1,uz1,ep1)
fcvy= fcvy -uxmid*uymid*del_y(j)

!pressure
prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
prmid=half*(tc2(i,j,k)+tc2(i,j+1,k))
fprx = fprx +prmid*del_y(j)

!viscous term
Expand Down Expand Up @@ -542,7 +543,7 @@ subroutine force(ux1,uy1,uz1,ep1)
fcvy= fcvy +uxmid*uymid*del_y(j)

!pressure
prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
prmid=half*(tc2(i,j,k)+tc2(i,j+1,k))
fprx = fprx -prmid*del_y(j)

!viscous term
Expand Down Expand Up @@ -680,7 +681,7 @@ subroutine force(ux1,uy1,uz1,ep1)
fcvy = fcvy + uxmid*uymid*del_y(j)*dz

!pressure
prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
prmid=half*(tc2(i,j,k)+tc2(i,j+1,k))
fprx = fprx -prmid*del_y(j)*dz

!viscous term
Expand Down Expand Up @@ -716,7 +717,7 @@ subroutine force(ux1,uy1,uz1,ep1)
fcvy = fcvy - uxmid*uymid*del_y(j)*dz

!pressure
prmid=half*(ppi2(i,j,k)+ppi2(i,j+1,k))
prmid=half*(tc2(i,j,k)+tc2(i,j+1,k))
fprx = fprx + prmid*del_y(j)*dz

!viscous term
Expand Down
Loading

0 comments on commit 46e874e

Please sign in to comment.