From 8ed9cdafcaf91f69bd8846d990228329868274a6 Mon Sep 17 00:00:00 2001 From: rfj82982 Date: Tue, 6 Aug 2024 10:29:36 +0000 Subject: [PATCH 1/4] Init ub parameter --- src/parameters.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parameters.f90 b/src/parameters.f90 index d2e16dfbd..44c28e434 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -641,6 +641,9 @@ subroutine parameter_defaults() !! IBM stuff nraf = 0 nobjmax = 0 + ubcx = zero + ubcy = zero + ubcz = zero nvol = 0 iforces = 0 From d5332e83f40e4dd2b899c61e23a637b6cec128c0 Mon Sep 17 00:00:00 2001 From: rfj82982 Date: Fri, 9 Aug 2024 16:18:20 +0000 Subject: [PATCH 2/4] Fix Intel issues with SGS modelling, IBM and statistics --- src/forces.f90 | 27 +++++----- src/genepsi3d.f90 | 119 +++++++++++++++++++++++-------------------- src/les_models.f90 | 73 +++++++++++++------------- src/module_param.f90 | 1 + src/navier.f90 | 14 +++-- src/statistics.f90 | 24 ++++++--- src/tools.f90 | 9 ++-- src/variables.f90 | 9 ++++ src/visu.f90 | 11 ++-- 9 files changed, 165 insertions(+), 122 deletions(-) diff --git a/src/forces.f90 b/src/forces.f90 index b9147de3d..cb3dfcd5e 100644 --- a/src/forces.f90 +++ b/src/forces.f90 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90 index 60a3c33ca..158f8615d 100644 --- a/src/genepsi3d.f90 +++ b/src/genepsi3d.f90 @@ -19,7 +19,7 @@ subroutine epsi_init(ep1) USE variables, only : yp, ny implicit none - real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1 + real(mytype),dimension(xsize(1),xsize(2),xsize(3)), intent(inout) :: ep1 logical :: dir_exists #ifdef DEBG @@ -59,7 +59,7 @@ subroutine geomcomplex(epsi, nxi, nxf, ny, nyi, nyf, nzi, nzf, dx, yp, dz, remp) IMPLICIT NONE INTEGER :: nxi,nxf,ny,nyi,nyf,nzi,nzf - REAL(mytype),DIMENSION(nxi:nxf,nyi:nyf,nzi:nzf) :: epsi + REAL(mytype),DIMENSION(nxi:nxf,nyi:nyf,nzi:nzf),intent(inout) :: epsi REAL(mytype) :: dx,dz REAL(mytype),DIMENSION(ny) :: yp REAL(mytype) :: remp @@ -111,7 +111,7 @@ subroutine genepsi3d(ep1) !*****************************************************************! ! logical :: dir_exists - real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1 + real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(inout) :: ep1 ! if (nrank==0.and.mod(itime,ilist)==0) then write(*,*)'===========================================================' @@ -155,29 +155,37 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& xi,xf,yi,yf,zi,zf,nobjx,nobjy,nobjz,& nobjmax,yp,nraf) use param, only : zero,half, one, two + use var, only : ta2, ta3 use decomp_2d use MPI + use complex_geometry, only: xepsi, yepsi, zepsi implicit none ! - real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1,smoofun,fbcx,fbcy,fbcz - real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: ep2 - real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: ep3 - integer :: nx,ny,nz,nobjmax - real(mytype) :: dx,dy,dz - real(mytype) :: xlx,yly,zlz - logical :: nclx,ncly,nclz - integer :: nxraf,nyraf,nzraf - integer :: nraf - integer, dimension(xsize(2),xsize(3)) :: nobjx,nobjxraf - integer, dimension(ysize(1),ysize(3)) :: nobjy,nobjyraf - integer, dimension(zsize(1),zsize(2)) :: nobjz,nobjzraf - real(mytype),dimension(nobjmax,xsize(2),xsize(3)) :: xi,xf - real(mytype),dimension(nobjmax,ysize(1),ysize(3)) :: yi,yf - real(mytype),dimension(nobjmax,zsize(1),zsize(2)) :: zi,zf - real(mytype),dimension(nxraf,xsize(2),xsize(3)) :: xepsi - real(mytype),dimension(ysize(1),nyraf,ysize(3)) :: yepsi - real(mytype),dimension(zsize(1),zsize(2),nzraf) :: zepsi - real(mytype),dimension(ny) :: yp + real(mytype),dimension(xsize(1),xsize(2),xsize(3)), intent(inout):: ep1 + !!real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1,smoofun,fbcx,fbcy,fbcz + !!real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: ep2 + !!real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: ep3 + integer, intent(in) :: nx,ny,nz,nobjmax + real(mytype), intent(in) :: dx,dy,dz + real(mytype), intent(in) :: xlx,yly,zlz + logical, intent(in) :: nclx,ncly,nclz + integer, intent(in) :: nxraf,nyraf,nzraf + real(mytype),dimension(nobjmax,xsize(2),xsize(3)):: xi,xf + real(mytype),dimension(nobjmax,ysize(1),ysize(3)):: yi,yf + real(mytype),dimension(nobjmax,zsize(1),zsize(2)):: zi,zf + integer :: nraf + integer, dimension(xsize(2),xsize(3)) :: nobjx + integer, dimension(ysize(1),ysize(3)) :: nobjy + integer, dimension(zsize(1),zsize(2)) :: nobjz + real(mytype),dimension(ny) :: yp + + ! Local variables + integer, dimension(xsize(2),xsize(3)) :: nobjxraf + integer, dimension(ysize(1),ysize(3)) :: nobjyraf + integer, dimension(zsize(1),zsize(2)) :: nobjzraf + !real(mytype),dimension(nxraf,xsize(2),xsize(3)) :: xepsi + !real(mytype),dimension(ysize(1),nyraf,ysize(3)) :: yepsi + !real(mytype),dimension(zsize(1),zsize(2),nzraf) :: zepsi real(mytype),dimension(nyraf) :: ypraf real(mytype) :: dxraf,dyraf,dzraf integer :: i,j,k @@ -189,7 +197,7 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& integer :: iflu,jflu,kflu integer :: isol,jsol,ksol integer :: iraf,jraf,kraf - integer :: nobjxmax ,nobjymax ,nobjzmax + integer :: nobjxmax,nobjymax ,nobjzmax integer :: nobjxmaxraf,nobjymaxraf,nobjzmaxraf integer :: idebraf,jdebraf,kdebraf integer :: ifinraf,jfinraf,kfinraf @@ -292,16 +300,16 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& !y-pencil nobjy(:,:)=0 nobjymax=0 - call transpose_x_to_y(ep1,ep2) + call transpose_x_to_y(ep1,ta2) do k=1,ysize(3) do i=1,ysize(1) jnum=0 - if(ep2(i,1,k) == one)then + if(ta2(i,1,k) == one)then jnum=1 nobjy(i,k)=1 endif do j=1,ny-1 - if(ep2(i,j,k) == zero .and. ep2(i,j+1,k) == one)then + if(ta2(i,j,k) == zero .and. ta2(i,j+1,k) == one)then jnum=jnum+1 nobjy(i,k)=nobjy(i,k)+1 endif @@ -348,16 +356,16 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& !z-pencil nobjz(:,:)=0 nobjzmax=0 - call transpose_y_to_z(ep2,ep3) + call transpose_y_to_z(ta2,ta3) do j=1,zsize(2) do i=1,zsize(1) knum=0 - if(ep3(i,j,1) == one)then + if(ta3(i,j,1) == one)then knum=1 nobjz(i,j)=1 endif do k=1,nz-1 - if(ep3(i,j,k) == zero .and. ep3(i,j,k+1) == one)then + if(ta3(i,j,k) == zero .and. ta3(i,j,k+1) == one)then knum=knum+1 nobjz(i,j)=nobjz(i,j)+1 endif @@ -496,11 +504,11 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& do i=1,ysize(1) if(nobjy(i,k) /= nobjyraf(i,k))then jobj=0 - if(ep2(i,1,k) == one)jobj=jobj+1 + if(ta2(i,1,k) == one)jobj=jobj+1 do j=1,ny-1 - if(ep2(i,j,k) == zero .and. ep2(i,j+1,k) == one)jobj=jobj+1 - if(ep2(i,j,k) == zero .and. ep2(i,j+1,k) == zero)jflu=1 - if(ep2(i,j,k) == one .and. ep2(i,j+1,k) == one)jsol=1 + if(ta2(i,j,k) == zero .and. ta2(i,j+1,k) == one)jobj=jobj+1 + if(ta2(i,j,k) == zero .and. ta2(i,j+1,k) == zero)jflu=1 + if(ta2(i,j,k) == one .and. ta2(i,j+1,k) == one)jsol=1 do jraf=1,nraf if(yepsi(i,jraf+nraf*(j-1) ,k) == zero .and.& yepsi(i,jraf+nraf*(j-1)+1,k) == one)jdebraf=jraf+nraf*(j-1)+1 @@ -565,11 +573,11 @@ subroutine gene_epsi_3D(ep1,nx,ny,nz,dx,dy,dz,xlx,yly,zlz ,& do i=1,zsize(1) if(nobjz(i,j) /= nobjzraf(i,j))then kobj=0 - if(ep3(i,j,1) == one)kobj=kobj+1 + if(ta3(i,j,1) == one)kobj=kobj+1 do k=1,nz-1 - if(ep3(i,j,k) == zero .and. ep3(i,j,k+1) == one)kobj=kobj+1 - if(ep3(i,j,k) == zero .and. ep3(i,j,k+1) == zero)kflu=1 - if(ep3(i,j,k) == one .and. ep3(i,j,k+1) == one)ksol=1 + if(ta3(i,j,k) == zero .and. ta3(i,j,k+1) == one)kobj=kobj+1 + if(ta3(i,j,k) == zero .and. ta3(i,j,k+1) == zero)kflu=1 + if(ta3(i,j,k) == one .and. ta3(i,j,k+1) == one)ksol=1 do kraf=1,nraf if(zepsi(i,j,kraf+nraf*(k-1) ) == zero .and.& zepsi(i,j,kraf+nraf*(k-1)+1) == one)kdebraf=kraf+nraf*(k-1)+1 @@ -614,14 +622,15 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& nxipif,nxfpif,nyipif,nyfpif,nzipif,nzfpif) use decomp_2d use param, only: zero, one + use var, only : ta2, ta3 use MPI implicit none ! integer :: nx,ny,nz,nobjmax real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ep1 - real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: ep2 - real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: ep3 + !real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: ep2 + !real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: ep3 integer,dimension(0:nobjmax,xsize(2),xsize(3)) :: nxipif,nxfpif integer,dimension(0:nobjmax,ysize(1),ysize(3)) :: nyipif,nyfpif integer,dimension(0:nobjmax,zsize(1),zsize(2)) :: nzipif,nzfpif @@ -675,7 +684,7 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& !if (nrank==0) write(*,*) ' step 11' !y-pencil - call transpose_x_to_y(ep1,ep2) + call transpose_x_to_y(ep1,ta2) nyipif(:,:,:)=npif nyfpif(:,:,:)=npif jsing=0 @@ -683,12 +692,12 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& do i=1,ysize(1) jnum=0 jflu=0 - if(ep2(i,1,k) == one)jnum=jnum+1 - if(ep2(i,1,k) == zero)jflu=jflu+1 + if(ta2(i,1,k) == one)jnum=jnum+1 + if(ta2(i,1,k) == zero)jflu=jflu+1 do j=2,ny - if(ep2(i,j ,k) == zero)jflu=jflu+1 - if(ep2(i,j-1,k) == zero .and.& - ep2(i,j ,k) == one)then + if(ta2(i,j ,k) == zero)jflu=jflu+1 + if(ta2(i,j-1,k) == zero .and.& + ta2(i,j ,k) == one)then jnum=jnum+1 if(jnum == 1)then nyipif(jnum ,i,k)=jflu-izap @@ -704,9 +713,9 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& jflu=0 endif endif - if(ep2(i,j,k) == one)jflu=0 + if(ta2(i,j,k) == one)jflu=0 enddo - if(ep2(i,ny,k) == zero)then + if(ta2(i,ny,k) == zero)then nyfpif(jnum,i,k)=jflu-izap if(jflu-izap < npif)jsing=jsing+1 if(jflu-izap < npif)nyfpif(jnum,i,k)=npif @@ -719,7 +728,7 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& !z-pencil if(nz > 1)then - call transpose_y_to_z(ep2,ep3) + call transpose_y_to_z(ta2,ta3) nzipif(:,:,:)=npif nzfpif(:,:,:)=npif ksing=0 @@ -727,12 +736,12 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& do i=1,zsize(1) knum=0 kflu=0 - if(ep3(i,j,1) == one)knum=knum+1 - if(ep3(i,j,1) == zero)kflu=kflu+1 + if(ta3(i,j,1) == one)knum=knum+1 + if(ta3(i,j,1) == zero)kflu=kflu+1 do k=2,nz - if(ep3(i,j,k ) == zero)kflu=kflu+1 - if(ep3(i,j,k-1) == zero .and.& - ep3(i,j,k ) == one)then + if(ta3(i,j,k ) == zero)kflu=kflu+1 + if(ta3(i,j,k-1) == zero .and.& + ta3(i,j,k ) == one)then knum=knum+1 if(knum == 1)then nzipif(knum ,i,j)=kflu-izap @@ -748,9 +757,9 @@ subroutine verif_epsi(ep1,npif,izap,nx,ny,nz,nobjmax,& kflu=0 endif endif - if(ep3(i,j,k) == one )kflu=0 + if(ta3(i,j,k) == one )kflu=0 enddo - if(ep3(i,j,nz) == zero)then + if(ta3(i,j,nz) == zero)then nzfpif(knum,i,j)=kflu-izap if(kflu-izap < npif)ksing=ksing+1 if(kflu-izap < npif)nzfpif(knum,i,j)=npif diff --git a/src/les_models.f90 b/src/les_models.f90 index aa0090181..56c7a2a68 100644 --- a/src/les_models.f90 +++ b/src/les_models.f90 @@ -106,15 +106,16 @@ subroutine compute_SGS(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,phi1,ep1) USE variables USE decomp_2d_io use var, only: nut1 - USE abl, only: wall_sgs_slip, wall_sgs_noslip + !USE abl, only: wall_sgs_slip, wall_sgs_noslip implicit none - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: ux1, uy1, uz1, ep1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar) :: phi1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: sgsx1, sgsy1, sgsz1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: wallfluxx1, wallfluxy1, wallfluxz1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: wallsgsx1, wallsgsy1, wallsgsz1 + real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ux1, uy1, uz1, ep1 + real(mytype), dimension(xsize(1), xsize(2), xsize(3), numscalar), intent(in) :: phi1 + real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(out) :: sgsx1, sgsy1, sgsz1 + !real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: w_fluxx1, w_fluxy1, w_fluxz1 + !real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: w_sgsx1, w_sgsy1, w_gsz1 + write(*,*) "START Compute_SGS" ! Calculate eddy-viscosity if(jles.eq.1) then ! Smagorinsky call smag(nut1,ux1,uy1,uz1) @@ -127,34 +128,38 @@ subroutine compute_SGS(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,phi1,ep1) endif - if(iconserv.eq.0) then ! Non-conservative form for calculating the divergence of the SGS stresses - call sgs_mom_nonconservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1,ep1) - - ! SGS correction for ABL - if(itype.eq.itype_abl) then - ! No-slip wall - if (ncly1==2) then - call wall_sgs_noslip(ux1,uy1,uz1,nut1,wallfluxx1,wallfluxy1,wallfluxz1) - if(xstart(2)==1) then - sgsx1(:,2,:) = -wallfluxx1(:,2,:) - sgsy1(:,2,:) = -wallfluxy1(:,2,:) - sgsz1(:,2,:) = -wallfluxz1(:,2,:) - endif - ! Slip wall - elseif (ncly1==1) then - call wall_sgs_slip(ux1,uy1,uz1,phi1,nut1,wallfluxx1,wallfluxy1,wallfluxz1) - if(xstart(2)==1) then - sgsx1(:,1,:) = wallfluxx1(:,1,:) - sgsy1(:,1,:) = wallfluxy1(:,1,:) - sgsz1(:,1,:) = wallfluxz1(:,1,:) - endif - endif - endif - - elseif (iconserv.eq.1) then ! Conservative form for calculating the divergence of the SGS stresses (used with wall functions) - call sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) - - endif + !if(iconserv.eq.0) then ! Non-conservative form for calculating the divergence of the SGS stresses + ! call sgs_mom_nonconservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1,ep1) + + ! ! SGS correction for ABL + ! if(itype.eq.itype_abl) then + ! ! No-slip wall + ! if (ncly1==2) then + ! call wall_sgs_noslip(ux1,uy1,uz1,nut1,wallfluxx1,wallfluxy1,wallfluxz1) + ! if(xstart(2)==1) then + ! sgsx1(:,2,:) = -wallfluxx1(:,2,:) + ! sgsy1(:,2,:) = -wallfluxy1(:,2,:) + ! sgsz1(:,2,:) = -wallfluxz1(:,2,:) + ! endif + ! ! Slip wall + ! elseif (ncly1==1) then + ! call wall_sgs_slip(ux1,uy1,uz1,phi1,nut1,wallfluxx1,wallfluxy1,wallfluxz1) + ! if(xstart(2)==1) then + ! sgsx1(:,1,:) = wallfluxx1(:,1,:) + ! sgsy1(:,1,:) = wallfluxy1(:,1,:) + ! sgsz1(:,1,:) = wallfluxz1(:,1,:) + ! endif + ! endif + ! endif + + !elseif (iconserv.eq.1) then ! Conservative form for calculating the divergence of the SGS stresses (used with wall functions) + ! call sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) + + !endif + sgsx1 = zero + sgsy1 = zero + sgsz1 = zero + write(*,*) "END Compute_SGS" return diff --git a/src/module_param.f90 b/src/module_param.f90 index 02f44bf17..db5085da4 100644 --- a/src/module_param.f90 +++ b/src/module_param.f90 @@ -529,6 +529,7 @@ module complex_geometry integer ,allocatable,dimension(:,:) :: nobjx,nobjy,nobjz integer ,allocatable,dimension(:,:,:) :: nxipif,nxfpif,nyipif,nyfpif,nzipif,nzfpif real(mytype),allocatable,dimension(:,:,:) :: xi,xf,yi,yf,zi,zf + real(mytype),allocatable,dimension(:,:,:) :: xepsi, yepsi, zepsi integer :: nxraf,nyraf,nzraf,nraf,nobjmax end module complex_geometry !############################################################################ diff --git a/src/navier.f90 b/src/navier.f90 index f5168ffa4..76c587c52 100644 --- a/src/navier.f90 +++ b/src/navier.f90 @@ -73,7 +73,7 @@ SUBROUTINE solve_poisson(pp3, px1, py1, pz1, rho1, ux1, uy1, uz1, ep1, drho1, di CALL momentum_to_velocity(rho1, ux1, uy1, uz1) ENDIF - CALL divergence(pp3(:,:,:,1),rho1,ux1,uy1,uz1,ep1,drho1,divu3,nlock) + call divergence(pp3(:,:,:,1),rho1,ux1,uy1,uz1,ep1,drho1,divu3,nlock) IF (ilmn.AND.ivarcoeff) THEN dv3(:,:,:) = pp3(:,:,:,1) ENDIF @@ -275,11 +275,11 @@ subroutine divergence (pp3,rho1,ux1,uy1,uz1,ep1,drho1,divu3,nlock) real(mytype),dimension(xsize(1),xsize(2),xsize(3),nrhotime),intent(in) :: rho1 !Z PENCILS NXM NYM NZ -->NXM NYM NZM real(mytype),dimension(zsize(1),zsize(2),zsize(3)),intent(in) :: divu3 - real(mytype),dimension(ph1%zst(1):ph1%zen(1),ph1%zst(2):ph1%zen(2),nzmsize) :: pp3 + real(mytype),dimension(ph1%zst(1):ph1%zen(1),ph1%zst(2):ph1%zen(2),nzmsize),intent(out) :: pp3 integer :: nvect3,i,j,k,nlock integer :: code - real(mytype) :: tmax,tmoy,tmax1,tmoy1 + real(mytype) :: tmax,tmoy,tmax1,tmoy1, pres_ref nvect3=(ph1%zen(1)-ph1%zst(1)+1)*(ph1%zen(2)-ph1%zst(2)+1)*nzmsize @@ -340,7 +340,11 @@ subroutine divergence (pp3,rho1,ux1,uy1,uz1,ep1,drho1,divu3,nlock) pp3(:,:,:) = pp3(:,:,:) + po3(:,:,:) if (nlock==2) then - pp3(:,:,:)=pp3(:,:,:)-pp3(ph1%zst(1),ph1%zst(2),nzmsize) + ! Line below sometimes generates issues with Intel + !pp3(:,:,:)=pp3(:,:,:)-pp3(ph1%zst(1),ph1%zst(2),nzmsize) + ! Using a tmp variable seems to sort the issue + pres_ref = pp3(ph1%zst(1),ph1%zst(2),nzmsize) + pp3(:,:,:)=pp3(:,:,:)-pres_ref endif tmax=-1609._mytype @@ -1185,7 +1189,7 @@ SUBROUTINE calc_varcoeff_rhs(pp3, rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, r tc1(:,:,:) = (one - rho0 / rho1(:,:,:,1)) * pz1(:,:,:) nlock = -1 !! Don't do any funny business with LMN - CALL divergence(pp3,rho1,ta1,tb1,tc1,ep1,drho1,divu3,nlock) + call divergence(pp3,rho1,ta1,tb1,tc1,ep1,drho1,divu3,nlock) !! lapl(p) = div((1 - rho0/rho) grad(p)) + rho0(div(u*) - div(u)) !! dv3 contains div(u*) - div(u) diff --git a/src/statistics.f90 b/src/statistics.f90 index 13d8c8c9a..daed31fd6 100644 --- a/src/statistics.f90 +++ b/src/statistics.f90 @@ -416,6 +416,7 @@ end subroutine update_average_vector subroutine update_variance_vector(uum, vvm, wwm, uvm, uwm, vwm, ux, uy, uz, ep) use decomp_2d, only : xsize, xstS, xenS + use var, only: ta1 implicit none @@ -423,12 +424,23 @@ subroutine update_variance_vector(uum, vvm, wwm, uvm, uwm, vwm, ux, uy, uz, ep) real(mytype), dimension(xstS(1):xenS(1),xstS(2):xenS(2),xstS(3):xenS(3)), intent(inout) :: uum, vvm, wwm, uvm, uwm, vwm real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ux, uy, uz, ep - call update_average_scalar(uum, ux*ux, ep) - call update_average_scalar(vvm, uy*uy, ep) - call update_average_scalar(wwm, uz*uz, ep) - call update_average_scalar(uvm, ux*uy, ep) - call update_average_scalar(uwm, ux*uz, ep) - call update_average_scalar(vwm, uy*uz, ep) + ta1(:,:,:) = ux(:,:,:)*ux(:,:,:) + call update_average_scalar(uum, ta1, ep) + + ta1(:,:,:) = uy(:,:,:)*uy(:,:,:) + call update_average_scalar(vvm, ta1, ep) + + ta1(:,:,:) = uz(:,:,:)*uz(:,:,:) + call update_average_scalar(wwm, ta1, ep) + + ta1(:,:,:) = ux(:,:,:)*uy(:,:,:) + call update_average_scalar(uvm, ta1, ep) + + ta1(:,:,:) = ux(:,:,:)*uz(:,:,:) + call update_average_scalar(uwm, ta1, ep) + + ta1(:,:,:) = uy(:,:,:)*uz(:,:,:) + call update_average_scalar(vwm, ta1, ep) end subroutine update_variance_vector diff --git a/src/tools.f90 b/src/tools.f90 index c1c2d2445..702b19c6d 100644 --- a/src/tools.f90 +++ b/src/tools.f90 @@ -590,16 +590,17 @@ end subroutine init_restart_adios2 subroutine apply_spatial_filter(ux1,uy1,uz1,phi1) use param - use var, only: uxf1,uyf1,uzf1,uxf2,uyf2,uzf2,uxf3,uyf3,uzf3,di1,di2,di3,phif1,phif2,phif3 + !use var, only: uxf1,uyf1,uzf1,uxf2,uyf2,uzf2,uxf3,uyf3,uzf3,di1,di2,di3,phif1,phif2,phif3 + use var, only: uxf1,uyf1,uzf1,di1,phif1 + use var, only: ux2,uy2,uz2, phi2,uxf2,uyf2,uzf2,di2,phif2 + use var, only: ux3,uy3,uz3, phi3,uxf3,uyf3,uzf3,di3,phif3 use variables use ibm_param, only : ubcx,ubcy,ubcz implicit none real(mytype),dimension(xsize(1),xsize(2),xsize(3)), intent(inout) :: ux1,uy1,uz1 - real(mytype),dimension(xsize(1),xsize(2),xsize(3), numscalar), intent(inout) :: phi1 + real(mytype),dimension(xsize(1),xsize(2),xsize(3), numscalar), intent(in) :: phi1 real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: phi11 - real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: ux2,uy2,uz2, phi2 - real(mytype),dimension(zsize(1),zsize(2),zsize(3)) :: ux3,uy3,uz3, phi3 integer :: i,j,k,npaire diff --git a/src/variables.f90 b/src/variables.f90 index 6523e1380..6b6e9d22d 100644 --- a/src/variables.f90 +++ b/src/variables.f90 @@ -85,6 +85,8 @@ module var subroutine init_variables + implicit none + TYPE(DECOMP_INFO), save :: ph! decomposition object integer :: i, j, k @@ -628,6 +630,13 @@ subroutine init_variables zi=zero allocate(zf(nobjmax,zsize(1),zsize(2))) zf=zero + allocate(xepsi(nxraf,xsize(2),xsize(3))) + xepsi=zero + allocate(yepsi(ysize(1),nyraf,ysize(3))) + yepsi=zero + allocate(zepsi(zsize(1),zsize(2),nzraf)) + zepsi=zero + endif !module filter diff --git a/src/visu.f90 b/src/visu.f90 index ac48582ae..7bcfd567a 100644 --- a/src/visu.f90 +++ b/src/visu.f90 @@ -485,6 +485,7 @@ subroutine write_field(f1, pathname, filename, num, skip_ibm, flush) use var, only : ep1 use var, only : zero, one use var, only : uvisu + use var, only : ta1 use param, only : iibm use decomp_2d_io, only : decomp_2d_write_one, decomp_2d_write_plane @@ -495,7 +496,7 @@ subroutine write_field(f1, pathname, filename, num, skip_ibm, flush) integer, intent(in) :: num logical, optional, intent(in) :: skip_ibm, flush - real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: local_array + !real(mytype), dimension(xsize(1),xsize(2),xsize(3)) :: local_array logical :: mpiio, force_flush integer :: ierr @@ -558,23 +559,23 @@ subroutine write_field(f1, pathname, filename, num, skip_ibm, flush) endif if ((iibm == 2) .and. .not.present(skip_ibm)) then - local_array(:,:,:) = (one - ep1(:,:,:)) * f1(:,:,:) + ta1(:,:,:) = (one - ep1(:,:,:)) * f1(:,:,:) else - local_array(:,:,:) = f1(:,:,:) + ta1(:,:,:) = f1(:,:,:) endif if (output2D.eq.0) then if (mpiio .or. (iibm == 2) .or. force_flush) then !! XXX: This (re)uses a temporary array for data - need to force synchronous writes. uvisu = zero - call fine_to_coarseV(1,local_array,uvisu) + call fine_to_coarseV(1,ta1,uvisu) call decomp_2d_write_one(1,uvisu,"data",gen_filename(pathname, filename, num, 'bin'),2,io_name,& opt_deferred_writes=.false.) else call decomp_2d_write_one(1,f1,"data",gen_filename(pathname, filename, num, 'bin'),0,io_name) end if else - call decomp_2d_write_plane(1,local_array,output2D,-1,"data",gen_filename(pathname, filename, num, 'bin'),io_name) + call decomp_2d_write_plane(1,ta1,output2D,-1,"data",gen_filename(pathname, filename, num, 'bin'),io_name) endif end subroutine write_field From 61620ff0f98f1d90641874bd7209b07373507b0f Mon Sep 17 00:00:00 2001 From: rfj82982 Date: Thu, 22 Aug 2024 10:06:20 +0000 Subject: [PATCH 3/4] Fix the explicit LES modelling by removing unnecessary large local arrays. --- examples/MHD/input_mhdchan_test.i3d | 2 +- src/Case-ABL.f90 | 83 ++++++++------ src/les_models.f90 | 169 ++++++++++++++-------------- src/variables.f90 | 10 ++ 4 files changed, 145 insertions(+), 119 deletions(-) diff --git a/examples/MHD/input_mhdchan_test.i3d b/examples/MHD/input_mhdchan_test.i3d index c0b787d01..e400b1f7b 100644 --- a/examples/MHD/input_mhdchan_test.i3d +++ b/examples/MHD/input_mhdchan_test.i3d @@ -12,7 +12,7 @@ p_row=0 ! Row partition p_col=0 ! Column partition ! Mesh -nx=8 ! X-direction nodes +nx=16 ! X-direction nodes ny=129 ! Y-direction nodes nz=1 ! Z-direction nodes istret = 2 ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom) diff --git a/src/Case-ABL.f90 b/src/Case-ABL.f90 index 39ba5fd2a..9e729d6d0 100644 --- a/src/Case-ABL.f90 +++ b/src/Case-ABL.f90 @@ -690,21 +690,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 @@ -718,47 +733,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 diff --git a/src/les_models.f90 b/src/les_models.f90 index 56c7a2a68..1d0d185c5 100644 --- a/src/les_models.f90 +++ b/src/les_models.f90 @@ -106,7 +106,8 @@ subroutine compute_SGS(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,phi1,ep1) USE variables USE decomp_2d_io use var, only: nut1 - !USE abl, only: wall_sgs_slip, wall_sgs_noslip + use var, only: ta1, tb1, tc1 + use abl, only: wall_sgs_slip, wall_sgs_noslip implicit none real(mytype), dimension(xsize(1), xsize(2), xsize(3)), intent(in) :: ux1, uy1, uz1, ep1 @@ -115,7 +116,6 @@ subroutine compute_SGS(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,phi1,ep1) !real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: w_fluxx1, w_fluxy1, w_fluxz1 !real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: w_sgsx1, w_sgsy1, w_gsz1 - write(*,*) "START Compute_SGS" ! Calculate eddy-viscosity if(jles.eq.1) then ! Smagorinsky call smag(nut1,ux1,uy1,uz1) @@ -128,38 +128,36 @@ subroutine compute_SGS(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,phi1,ep1) endif - !if(iconserv.eq.0) then ! Non-conservative form for calculating the divergence of the SGS stresses - ! call sgs_mom_nonconservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1,ep1) - - ! ! SGS correction for ABL - ! if(itype.eq.itype_abl) then - ! ! No-slip wall - ! if (ncly1==2) then - ! call wall_sgs_noslip(ux1,uy1,uz1,nut1,wallfluxx1,wallfluxy1,wallfluxz1) - ! if(xstart(2)==1) then - ! sgsx1(:,2,:) = -wallfluxx1(:,2,:) - ! sgsy1(:,2,:) = -wallfluxy1(:,2,:) - ! sgsz1(:,2,:) = -wallfluxz1(:,2,:) - ! endif - ! ! Slip wall - ! elseif (ncly1==1) then - ! call wall_sgs_slip(ux1,uy1,uz1,phi1,nut1,wallfluxx1,wallfluxy1,wallfluxz1) - ! if(xstart(2)==1) then - ! sgsx1(:,1,:) = wallfluxx1(:,1,:) - ! sgsy1(:,1,:) = wallfluxy1(:,1,:) - ! sgsz1(:,1,:) = wallfluxz1(:,1,:) - ! endif - ! endif - ! endif - - !elseif (iconserv.eq.1) then ! Conservative form for calculating the divergence of the SGS stresses (used with wall functions) - ! call sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) - - !endif - sgsx1 = zero - sgsy1 = zero - sgsz1 = zero - write(*,*) "END Compute_SGS" + ! we can store wallflux in ta1, tb1, tc1 + ! this will be done in all call to wall_sgs_ for the ABL + if(iconserv.eq.0) then ! Non-conservative form for calculating the divergence of the SGS stresses + call sgs_mom_nonconservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1,ep1) + + ! SGS correction for ABL + if(itype.eq.itype_abl) then + ! No-slip wall + if (ncly1==2) then + call wall_sgs_noslip(ux1,uy1,uz1,nut1,ta1,tb1,tc1) + if(xstart(2)==1) then + sgsx1(:,2,:) = -ta1(:,2,:) + sgsy1(:,2,:) = -tb1(:,2,:) + sgsz1(:,2,:) = -tc1(:,2,:) + endif + ! Slip wall + elseif (ncly1==1) then + call wall_sgs_slip(ux1,uy1,uz1,phi1,nut1,ta1,tb1,tc1) + if(xstart(2)==1) then + sgsx1(:,1,:) = ta1(:,1,:) + sgsy1(:,1,:) = tb1(:,1,:) + sgsz1(:,1,:) = tc1(:,1,:) + endif + endif + endif + + elseif (iconserv.eq.1) then ! Conservative form for calculating the divergence of the SGS stresses (used with wall functions) + call sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) + + endif return @@ -1112,15 +1110,15 @@ subroutine sgs_mom_nonconservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1,ep1) implicit none - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: ux1, uy1, uz1, nut1, ep1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: sgsx1, sgsy1, sgsz1 + real(mytype), dimension(xsize(1), xsize(2), xsize(3)),intent(in) :: ux1, uy1, uz1, nut1, ep1 + real(mytype), dimension(xsize(1), xsize(2), xsize(3)),intent(out):: sgsx1, sgsy1, sgsz1 integer :: i, j, k, ijk, nvect1 ta1 = zero; ta2 = zero; ta3 = zero - sgsx1=zero;sgsy1=zero;sgsz1=zero - sgsx2=zero;sgsy2=zero;sgsz2=zero - sgsx3=zero;sgsy3=zero;sgsz3=zero + sgsx1=zero; sgsy1=zero; sgsz1=zero + sgsx2=zero; sgsy2=zero; sgsz2=zero + sgsx3=zero; sgsy3=zero; sgsz3=zero !WORK X-PENCILS call derx (ta1,nut1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,zero) @@ -1330,8 +1328,11 @@ subroutine sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) USE variables use MPI USE var, only : ta1,tb1,tc1,di1 + use var, only : td1, te1, tf1, tg1, th1, ti1 USE var, only : ta2,tb2,tc2,di2 + use var, only : te2, tf2, tg2, th2, ti2 USE var, only : ta3,tb3,tc3,di3 + use var, only : tf3, th3, ti3 USE var, only : sgsx2,sgsy2,sgsz2 USE var, only : sgsx3,sgsy3,sgsz3 USE var, only : sxx1,sxy1,sxz1,syy1,syz1,szz1 @@ -1340,37 +1341,36 @@ subroutine sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) implicit none - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: ux1, uy1, uz1, nut1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: sgsx1, sgsy1, sgsz1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: wallsgsx1, wallsgsy1, wallsgsz1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: txx1, txy1, txz1, tyy1, tyz1, tzz1 - real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: taf1, tbf1, tcf1 - real(mytype), dimension(ysize(1), ysize(2), ysize(3)) :: txy2, txz2, tyy2, tyz2, tzz2 - real(mytype), dimension(ysize(1), ysize(2), ysize(3)) :: taf2, tbf2, tcf2 - real(mytype), dimension(zsize(1), zsize(2), zsize(3)) :: txz3, tyz3, tzz3 - real(mytype), dimension(zsize(1), zsize(2), zsize(3)) :: taf3, tbf3, tcf3 + 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) :: sgsx1, sgsy1, sgsz1 + ! ta1 tb1 tc1 + !real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: wallsgsx1, wallsgsy1, wallsgsz1 + ! + ! td1 te1 tf1 tg1 th1 ti1 + !real(mytype), dimension(xsize(1), xsize(2), xsize(3)) :: txx, txy, txz, tyy, tyz, tzz + integer :: i, j, k, code, ierr ! Construct stress tensor - txx1 = 2.0*nut1*sxx1 - txy1 = 2.0*nut1*sxy1 - txz1 = 2.0*nut1*sxz1 - tyy1 = 2.0*nut1*syy1 - tyz1 = 2.0*nut1*syz1 - tzz1 = 2.0*nut1*szz1 + td1 = 2.0*nut1*sxx1 + te1 = 2.0*nut1*sxy1 + tf1 = 2.0*nut1*sxz1 + tg1 = 2.0*nut1*syy1 + th1 = 2.0*nut1*syz1 + ti1 = 2.0*nut1*szz1 ! Add wall model for ABL if (itype.eq.itype_abl) then if (ncly1==2) then - call wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1) + call wall_sgs_noslip(ux1,uy1,uz1,nut1,ta1,tb1,tc1) if (xstart(2)==1) then - txx1(:,2,:) = 0. - txy1(:,2,:) = - wallsgsx1(:,2,:)! txy1(:,2,:) - txz1(:,2,:) = 0. - tyy1(:,2,:) = 0. - tyz1(:,2,:) = - wallsgsz1(:,2,:)! tyz1(:,2,:) - tzz1(:,2,:) = 0. + td1(:,2,:) = 0. + te1(:,2,:) = - ta1(:,2,:)! txy1(:,2,:) + tf1(:,2,:) = 0. + tg1(:,2,:) = 0. + th1(:,2,:) = - tc1(:,2,:)! tyz1(:,2,:) + ti1(:,2,:) = 0. endif elseif (ncly1==1) then write(*,*) 'Simulation stopped: slip bottom wall not supported with iconserv=1' @@ -1382,37 +1382,38 @@ subroutine sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) ta1 = zero; ta2 = zero; ta3 = zero tb1 = zero; tb2 = zero; tb3 = zero tc1 = zero; tc2 = zero; tc3 = zero - sgsx1=0.;sgsy1=0.;sgsz1=0. - sgsx2=0.;sgsy2=0.;sgsz2=0. - sgsx3=0.;sgsy3=0.;sgsz3=0. + sgsx1=zero; sgsy1=zero; sgsz1=zero + sgsx2=zero; sgsy2=zero; sgsz2=zero + sgsx3=zero; sgsy3=zero; sgsz3=zero ! WORK X-PENCILS - call derx (ta1,txx1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,zero) - call derx (tb1,txy1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,zero) - call derx (tc1,txz1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,zero) + call derx (sgsx1,td1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,zero) + call derx (sgsy1,te1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,zero) + call derx (sgsz1,tf1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,zero) !call filter(0.48d0) !call filx(taf1,ta1,di1,fisx,fiffxp,fifsxp,fifwxp,xsize(1),xsize(2),xsize(3),1,zero) !call filx(tbf1,tb1,di1,fisx,fiffxp,fifsxp,fifwxp,xsize(1),xsize(2),xsize(3),1,zero) !call filx(tcf1,tc1,di1,fisx,fiffxp,fifsxp,fifwxp,xsize(1),xsize(2),xsize(3),1,zero) - sgsx1 = ta1 - sgsy1 = tb1 - sgsz1 = tc1 + ! Stored into var immediatly not here + !sgsx1 = ta1 + !sgsy1 = tb1 + !sgsz1 = tc1 ! WORK Y-PENCILS - call transpose_x_to_y(txy1, txy2) - call transpose_x_to_y(tyy1, tyy2) - call transpose_x_to_y(tyz1, tyz2) - call transpose_x_to_y(txz1, txz2) - call transpose_x_to_y(tzz1, tzz2) + call transpose_x_to_y(te1, te2) + call transpose_x_to_y(tg1, tg2) + call transpose_x_to_y(th1, th2) + call transpose_x_to_y(tf1, tf2) + call transpose_x_to_y(ti1, ti2) call transpose_x_to_y(sgsx1, sgsx2) call transpose_x_to_y(sgsy1, sgsy2) call transpose_x_to_y(sgsz1, sgsz2) - call dery (ta2,txy2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,zero) - call dery (tb2,tyy2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,zero) - call dery (tc2,tyz2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,zero) + call dery (ta2,te2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,zero) + call dery (tb2,tg2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,zero) + call dery (tc2,th2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,zero) !call fily(taf2,ta2,di2,fisy,fiffyp,fifsyp,fifwyp,ysize(1),ysize(2),ysize(3),1,zero) !call fily(tbf2,tb2,di2,fisy,fiffyp,fifsyp,fifwyp,ysize(1),ysize(2),ysize(3),1,zero) @@ -1426,13 +1427,13 @@ subroutine sgs_mom_conservative(sgsx1,sgsy1,sgsz1,ux1,uy1,uz1,nut1) call transpose_y_to_z(sgsx2, sgsx3) call transpose_y_to_z(sgsy2, sgsy3) call transpose_y_to_z(sgsz2, sgsz3) - call transpose_y_to_z(txz2, txz3) - call transpose_y_to_z(tyz2, tyz3) - call transpose_y_to_z(tzz2, tzz3) + call transpose_y_to_z(tf2, tf3) + call transpose_y_to_z(th2, th3) + call transpose_y_to_z(ti2, ti3) - call derz (ta3, txz3, di3, sz, ffz, fsz, fwz, zsize(1), zsize(2), zsize(3), 0, zero) - call derz (tb3, tyz3, di3, sz, ffz, fsz, fwz, zsize(1), zsize(2), zsize(3), 0, zero) - call derz (tc3, tzz3, di3, sz, ffz, fsz, fwz, zsize(1), zsize(2), zsize(3), 0, zero) + call derz (ta3, tf3, di3, sz, ffz, fsz, fwz, zsize(1), zsize(2), zsize(3), 0, zero) + call derz (tb3, th3, di3, sz, ffz, fsz, fwz, zsize(1), zsize(2), zsize(3), 0, zero) + call derz (tc3, ti3, di3, sz, ffz, fsz, fwz, zsize(1), zsize(2), zsize(3), 0, zero) !call filz(taf3,ta3,di3,fisz,fiffzp,fifszp,fifwzp,zsize(1),zsize(2),zsize(3),1,zero) !call filz(tbf3,tb3,di3,fisz,fiffzp,fifszp,fifwzp,zsize(1),zsize(2),zsize(3),1,zero) diff --git a/src/variables.f90 b/src/variables.f90 index 6b6e9d22d..811789bfc 100644 --- a/src/variables.f90 +++ b/src/variables.f90 @@ -53,6 +53,7 @@ module var integer, save :: nxmsize, nymsize, nzmsize ! working arrays for LES + ! These are by far too many variables and we should be able to decrease them real(mytype), save, allocatable, dimension(:,:,:) :: sgsx1,sgsy1,sgsz1,nut1,sxx1,syy1,szz1,sxy1,sxz1,syz1 real(mytype), save, allocatable, dimension(:,:,:) :: sgsx2,sgsy2,sgsz2,nut2,sxx2,syy2,szz2,sxy2,sxz2,syz2 real(mytype), save, allocatable, dimension(:,:,:) :: sgsx3,sgsy3,sgsz3,nut3,sxx3,syy3,szz3,sxy3,sxz3,syz3 @@ -69,9 +70,12 @@ module var real(mytype), save, allocatable, dimension(:,:,:) :: srt_smag, srt_smag2 real(mytype), save, allocatable, dimension(:,:,:) :: srt_wale, srt_wale2, srt_wale3, srt_wale4 + ! working arrays for ABL real(mytype), save, allocatable, dimension(:,:) :: heatflux real(mytype), save, allocatable, dimension(:,:,:,:) :: abl_T + ! tmp fix for intel + real(mytype), save, allocatable, dimension(:,:,:) :: wallfluxx1, wallfluxy1, wallfluxz1 ! arrays for turbine modelling real(mytype), save, allocatable, dimension(:,:,:) :: FTx, FTy, FTz, Fdiscx, Fdiscy, Fdiscz @@ -474,6 +478,12 @@ subroutine init_variables srt_smag=zero call alloc_x(srt_wale) srt_wale=zero + call alloc_x(wallfluxx1) + wallfluxx1=zero + call alloc_x(wallfluxy1) + wallfluxy1=zero + call alloc_x(wallfluxz1) + wallfluxz1=zero call alloc_y(sgsx2) sgsx2=zero call alloc_y(sgsy2) From 677499bf0f172aa86eb6dd42e417a7ec8874e4e4 Mon Sep 17 00:00:00 2001 From: rfj82982 Date: Thu, 22 Aug 2024 13:36:57 +0000 Subject: [PATCH 4/4] Add adios2 config file to TGV test --- .../Taylor-Green-Vortex/adios2_config.xml | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 test/data/Taylor-Green-Vortex/adios2_config.xml diff --git a/test/data/Taylor-Green-Vortex/adios2_config.xml b/test/data/Taylor-Green-Vortex/adios2_config.xml new file mode 100644 index 000000000..20a156f18 --- /dev/null +++ b/test/data/Taylor-Green-Vortex/adios2_config.xml @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +