Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove dbg schemes #267

Merged
merged 2 commits into from
Apr 21, 2024
Merged

Remove dbg schemes #267

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
909090c
Remove kind-specific functions
mathrack Apr 21, 2024
0693919
Move the subroutines rl, iy and cx in a dedicated module
mathrack Apr 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ add_executable(xcompact3d
tools.f90
transeq.f90
turbine.f90
utilities.f90
variables.f90
visu.f90
xcompact3d.f90)
Expand Down
53 changes: 24 additions & 29 deletions src/Case-ABL.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,6 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz)
use var, only: di1, di3
use var, only: sxy1, syz1, heatflux, ta2, tb2, ta3, tb3
use ibm_param, only : ubcx, ubcz
use dbg_schemes, only: log_prec, tanh_prec, sqrt_prec, abs_prec, atan_prec

implicit none
real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux,uy,uz, nut1
Expand Down Expand Up @@ -512,19 +511,19 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz)
if (itherm==0) then
Q_HAve = TempFlux
else if (itherm==1) then
Q_HAve =-k_roughness**two*S_HAve*(Phi_HAve-(T_wall+TempRate*t))/((log_prec(delta/z_zero)-PsiM_HAve)*(log_prec(delta/z_zero)-PsiH_HAve))
Q_HAve =-k_roughness**two*S_HAve*(Phi_HAve-(T_wall+TempRate*t))/((log(delta/z_zero)-PsiM_HAve)*(log(delta/z_zero)-PsiH_HAve))
endif
L_HAve=-(k_roughness*S_HAve/(log_prec(delta/z_zero)-PsiM_HAve))**three*Phi_HAve/(k_roughness*gravv*Q_HAve)
L_HAve=-(k_roughness*S_HAve/(log(delta/z_zero)-PsiM_HAve))**three*Phi_HAve/(k_roughness*gravv*Q_HAve)
if (istrat==0) then
PsiM_HAve=-4.8_mytype*delta/L_HAve
PsiH_HAve=-7.8_mytype*delta/L_HAve
else if (istrat==1) then
zeta_HAve=(one-sixteen*delta/L_HAve)**zptwofive
PsiM_HAve=two*log_prec(half*(one+zeta_HAve))+log_prec(zpfive*(one+zeta_HAve**two))-two*atan_prec(zeta_HAve)+pi/two
PsiH_HAve=two*log_prec(half*(one+zeta_HAve**two))
PsiM_HAve=two*log(half*(one+zeta_HAve))+log(zpfive*(one+zeta_HAve**two))-two*atan(zeta_HAve)+pi/two
PsiH_HAve=two*log(half*(one+zeta_HAve**two))
endif
ii = ii + 1
OL_diff = abs_prec((L_HAve - Lold)/Lold)
OL_diff = abs((L_HAve - Lold)/Lold)
Lold = L_HAve
if (ii==50) exit
enddo
Expand All @@ -548,30 +547,30 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz)
do i=1,xsize(1)
! Horizontally-averaged formulation
if(iwallmodel==1) then
tauwallxy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM_HAve))**two*ux_HAve*S_HAve
tauwallzy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM_HAve))**two*uz_HAve*S_HAve
tauwallxy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM_HAve))**two*ux_HAve*S_HAve
tauwallzy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM_HAve))**two*uz_HAve*S_HAve
! Local formulation
else
ux12=half*(uxf1(i,1,k)+uxf1(i,2,k))
uz12=half*(uzf1(i,1,k)+uzf1(i,2,k))
S12=sqrt_prec(ux12**2.+uz12**2.)
S12=sqrt(ux12**2.+uz12**2.)
if (iscalar==1) then
Phi12= half*(phif1(i,1,k)+ phif1(i,2,k)) + Tstat12
do ii=1,10
if (itherm==1) heatflux(i,k)=-k_roughness**two*S12*(Phi12-(T_wall+TempRate*t))/((log_prec(delta/z_zero)-PsiM(i,k))*(log_prec(delta/z_zero)-PsiH(i,k)))
Obukhov(i,k)=-(k_roughness*S12/(log_prec(delta/z_zero)-PsiM(i,k)))**three*Phi12/(k_roughness*gravv*heatflux(i,k))
if (itherm==1) heatflux(i,k)=-k_roughness**two*S12*(Phi12-(T_wall+TempRate*t))/((log(delta/z_zero)-PsiM(i,k))*(log(delta/z_zero)-PsiH(i,k)))
Obukhov(i,k)=-(k_roughness*S12/(log(delta/z_zero)-PsiM(i,k)))**three*Phi12/(k_roughness*gravv*heatflux(i,k))
if (istrat==0) then
PsiM(i,k)=-4.8_mytype*delta/Obukhov(i,k)
PsiH(i,k)=-7.8_mytype*delta/Obukhov(i,k)
else if (istrat==1) then
zeta(i,k)=(one-sixteen*delta/Obukhov(i,k))**zptwofive
PsiM(i,k)=two*log_prec(half*(one+zeta(i,k)))+log_prec(zpfive*(one+zeta(i,k)**2.))-two*atan_prec(zeta(i,k))+pi/two
PsiH(i,k)=two*log_prec(half*(one+zeta(i,k)**two))
PsiM(i,k)=two*log(half*(one+zeta(i,k)))+log(zpfive*(one+zeta(i,k)**2.))-two*atan(zeta(i,k))+pi/two
PsiH(i,k)=two*log(half*(one+zeta(i,k)**two))
endif
enddo
endif
tauwallxy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM(i,k)))**two*ux12*S12
tauwallzy(i,k)=-(k_roughness/(log_prec(delta/z_zero)-PsiM(i,k)))**two*uz12*S12
tauwallxy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM(i,k)))**two*ux12*S12
tauwallzy(i,k)=-(k_roughness/(log(delta/z_zero)-PsiM(i,k)))**two*uz12*S12
endif
! Apply second-order upwind scheme for the near wall
! Below should change for non-uniform grids, same for wall_sgs_slip_scalar
Expand Down Expand Up @@ -617,7 +616,7 @@ subroutine wall_sgs_slip(ux,uy,uz,phi,nut1,wallfluxx,wallfluxy,wallfluxz)
endif

! Compute friction velocity u_shear and boundary layer height
u_shear=k_roughness*S_HAve/(log_prec(delta/z_zero)-PsiM_HAve)
u_shear=k_roughness*S_HAve/(log(delta/z_zero)-PsiM_HAve)
if (iheight==1) call boundary_height(ux,uy,uz,dBL)
if (iscalar==1) zL=dBL/L_HAve

Expand Down Expand Up @@ -692,7 +691,6 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1)
use variables
use var, only: di1, di2, di3
use var, only: sxy1, syz1, tb1, ta2, tb2
use dbg_schemes, only: log_prec, sqrt_prec
use ibm_param, only : ubcx, ubcy, ubcz

implicit none
Expand Down Expand Up @@ -738,9 +736,9 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1)
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_prec(ux_delta**2.+uz_delta**2.)
tauwallxy2(i,k)=-(k_roughness/(log_prec(delta/z_zero)))**two*ux_delta*S_delta
tauwallzy2(i,k)=-(k_roughness/(log_prec(delta/z_zero)))**two*uz_delta*S_delta
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
Expand Down Expand Up @@ -809,8 +807,8 @@ subroutine wall_sgs_noslip(ux1,uy1,uz1,nut1,wallsgsx1,wallsgsy1,wallsgsz1)
call MPI_ALLREDUCE(uz_HAve_local,uz_HAve,1,real_type,MPI_SUM,MPI_COMM_WORLD,code)
ux_HAve=ux_HAve/(nxc*nzc)
uz_HAve=uz_HAve/(nxc*nzc)
S_HAve=sqrt_prec(ux_HAve**2.+uz_HAve**2.)
u_shear=k_roughness*S_HAve/log_prec(5*dy/z_zero)
S_HAve=sqrt(ux_HAve**2.+uz_HAve**2.)
u_shear=k_roughness*S_HAve/log(5*dy/z_zero)

if (mod(itime,ilist)==0.and.nrank==0) then
! Write u_shear in file
Expand All @@ -835,7 +833,6 @@ subroutine forceabl (ux) ! Routine to force constant flow rate
use param
use var
use MPI
use dbg_schemes, only: log_prec

implicit none

Expand Down Expand Up @@ -995,7 +992,6 @@ subroutine damping_zone (dux1,duy1,duz1,ux1,uy1,uz1) ! Damping zone for ABL

use param
use var, only: yp
use dbg_schemes, only: sin_prec, cos_prec, log_prec

implicit none

Expand All @@ -1022,7 +1018,7 @@ subroutine damping_zone (dux1,duy1,duz1,ux1,uy1,uz1) ! Damping zone for ABL
! Damping for non-neutral ABL
if (ibuoyancy==1) then
if (y>=damp_lo) then
lambda=sin_prec(half*pi*(y-damp_lo)/(yly-damp_lo))**two
lambda=sin(half*pi*(y-damp_lo)/(yly-damp_lo))**two
else
lambda=zero
endif
Expand All @@ -1034,13 +1030,13 @@ subroutine damping_zone (dux1,duy1,duz1,ux1,uy1,uz1) ! Damping zone for ABL
if (y>=(dBL+half*dheight)) then
lambda=one
elseif (y>=(dBL-half*dheight).and.y<(dBL+half*dheight)) then
lambda=half*(one-cos_prec(pi*(y-(dBL-half*dheight))/dheight))
lambda=half*(one-cos(pi*(y-(dBL-half*dheight))/dheight))
else
lambda=zero
endif
xloc=real(i-1,mytype)*dx
if (iconcprec.eq.1.and.xloc.ge.pdl) lambda=0.
dux1(i,j,k,1)=dux1(i,j,k,1)-coeff*lambda*(ux1(i,j,k)-ustar/k_roughness*log_prec(dBL/z_zero))
dux1(i,j,k,1)=dux1(i,j,k,1)-coeff*lambda*(ux1(i,j,k)-ustar/k_roughness*log(dBL/z_zero))
duy1(i,j,k,1)=duy1(i,j,k,1)-coeff*lambda*(uy1(i,j,k)-UG(2))
duz1(i,j,k,1)=duz1(i,j,k,1)-coeff*lambda*(uz1(i,j,k)-UG(3))
endif
Expand All @@ -1058,7 +1054,6 @@ subroutine damping_zone_scalar (dphi1,phi1)

use param
use var, only: yp
use dbg_schemes, only: sin_prec

implicit none

Expand All @@ -1077,7 +1072,7 @@ subroutine damping_zone_scalar (dphi1,phi1)
if (istret==0) y=real(j+xstart(2)-1-1,mytype)*dy
if (istret/=0) y=yp(j+xstart(2)-1)
if (y>=damp_lo) then
lambda=sin_prec(half*pi*(y-damp_lo)/(yly-damp_lo))**two
lambda=sin(half*pi*(y-damp_lo)/(yly-damp_lo))**two
else
lambda=zero
endif
Expand Down
7 changes: 2 additions & 5 deletions src/Case-Cylinder-wake.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ subroutine geomcomplex_cyl(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)

use param, only : one, two, ten
use ibm_param
use dbg_schemes, only: sqrt_prec

implicit none

Expand Down Expand Up @@ -69,7 +68,7 @@ subroutine geomcomplex_cyl(epsi,nxi,nxf,ny,nyi,nyf,nzi,nzf,dx,yp,remp)
ym=yp(j)
do i=nxi,nxf
xm=real(i-1,mytype)*dx
r=sqrt_prec((xm-cexx)**two+(ym-ceyy)**two)
r=sqrt((xm-cexx)**two+(ym-ceyy)**two)
if (r-ra.gt.zeromach) then
cycle
endif
Expand Down Expand Up @@ -209,7 +208,6 @@ subroutine init_cyl (ux1,uy1,uz1,phi1)
USE variables
USE param
USE MPI
use dbg_schemes, only: exp_prec

implicit none

Expand Down Expand Up @@ -252,7 +250,7 @@ subroutine init_cyl (ux1,uy1,uz1,phi1)
do j=1,xsize(2)
if (istret.eq.0) y=(j+xstart(2)-1-1)*dy-yly/2.
if (istret.ne.0) y=yp(j+xstart(2)-1)-yly/2.
um=exp_prec(-zptwo*y*y)
um=exp(-zptwo*y*y)
do i=1,xsize(1)
ux1(i,j,k)=um*ux1(i,j,k)
uy1(i,j,k)=um*uy1(i,j,k)
Expand Down Expand Up @@ -290,7 +288,6 @@ subroutine postprocess_cyl(ux1,uy1,uz1,ep1)
USE var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1
USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
USE ibm_param
use dbg_schemes, only: sqrt_prec

real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) :: ux1, uy1, uz1, ep1

Expand Down
21 changes: 10 additions & 11 deletions src/Case-Mixing-layer.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1)
use param, only : u1, u2, dens1, dens2
use param, only : half, one, two, four, eight, sixteen
use param, only : ntime, nrhotime
use dbg_schemes, only: sin_prec, cos_prec, exp_prec, tanh_prec, sqrt_prec
use MPI

implicit none
Expand Down Expand Up @@ -57,8 +56,8 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1)
rhomax = max(dens1, dens2)
T1 = pressure0 / dens1
T2 = pressure0 / dens2
u1 = sqrt_prec(dens2 / dens1) / (sqrt_prec(dens2 / dens1) + one)
u2 = -sqrt_prec(dens1 / dens2) / (one + sqrt_prec(dens1 / dens2))
u1 = sqrt(dens2 / dens1) / (sqrt(dens2 / dens1) + one)
u2 = -sqrt(dens1 / dens2) / (one + sqrt(dens1 / dens2))
M = 0.2_mytype
rspech = 1.4_mytype
heatcap = (one / (T2 * (rspech - one))) * ((u1 - u2) / M)**2
Expand All @@ -71,7 +70,7 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1)

!! Set mean field
ux1(i, j, k) = ux1(i, j, k) + half * (u1 + u2) &
+ half * (u1 - u2) * tanh_prec(two * y)
+ half * (u1 - u2) * tanh(two * y)
uy1(i, j, k) = zero
uz1(i, j, k) = zero

Expand All @@ -84,14 +83,14 @@ subroutine init_mixlayer (rho1,ux1,uy1,uz1)
rho1(i, j, k, 1) = MIN(rho1(i, j, k, 1), rhomax)

! Calculate disturbance field (as given in Fortune2004)
disturb_decay = 0.025_mytype * (u1 - u2) * exp_prec(-0.05_mytype * (y**2))
u_disturb = disturb_decay * (sin_prec(eight * PI * x / xlx) &
+ sin_prec(four * PI * x / xlx) / eight &
+ sin_prec(two * PI * x / xlx) / sixteen)
disturb_decay = 0.025_mytype * (u1 - u2) * exp(-0.05_mytype * (y**2))
u_disturb = disturb_decay * (sin(eight * PI * x / xlx) &
+ sin(four * PI * x / xlx) / eight &
+ sin(two * PI * x / xlx) / sixteen)
u_disturb = (0.05_mytype * y * xlx / PI) * u_disturb
v_disturb = disturb_decay * (cos_prec(eight * PI * x / xlx) &
+ cos_prec(four * PI * x / xlx) / eight &
+ cos_prec(two * PI * x / xlx) / sixteen)
v_disturb = disturb_decay * (cos(eight * PI * x / xlx) &
+ cos(four * PI * x / xlx) / eight &
+ cos(two * PI * x / xlx) / sixteen)

ux1(i, j, k) = ux1(i, j, k) + u_disturb
uy1(i, j, k) = uy1(i, j, k) + v_disturb
Expand Down
11 changes: 5 additions & 6 deletions src/Case-TBL.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,6 @@ subroutine blasius()
use decomp_2d_io
use MPI
use param, only : zero, zptwo, zpeight, one, nine
use dbg_schemes, only: exp_prec, sqrt_prec

implicit none

Expand Down Expand Up @@ -213,7 +212,7 @@ subroutine blasius()
-17.018112827391400_mytype*eta_bl**4 +0.819582894357566_mytype*eta_bl**3 &
-0.0601348202321954_mytype*eta_bl**2 +2.989739912704050_mytype*eta_bl**1

f_bl=f_bl+(1-exp_prec(-delta_eta/delta_int))*eps_eta
f_bl=f_bl+(1-exp(-delta_eta/delta_int))*eps_eta



Expand All @@ -232,14 +231,14 @@ subroutine blasius()
+130.109496961069_mytype*eta_bl**4 -7.64507811014497000_mytype*eta_bl**3 &
+6.94303207046209_mytype*eta_bl**2 -0.00209716712558639_mytype*eta_bl**1 ! &

g_bl=g_bl+(1-exp_prec(-delta_eta/delta_int))*eps_eta
g_bl=g_bl+(1-exp(-delta_eta/delta_int))*eps_eta



x_bl=one/(4.91_mytype**2*xnu)

bxx1(j,k)=f_bl/1.0002014996204402_mytype/1.0000000359138641_mytype !To assure 1.0 in infinity
bxy1(j,k)=g_bl*sqrt_prec(xnu/x_bl)/1.000546554_mytype
bxy1(j,k)=g_bl*sqrt(xnu/x_bl)/1.000546554_mytype
bxz1(j,k)=zero
enddo
enddo
Expand Down Expand Up @@ -269,7 +268,7 @@ subroutine blasius()
-0.0601348202321954_mytype*eta_bl**2 +2.989739912704050_mytype*eta_bl**1


f_bl_inf=f_bl_inf+(1-exp_prec(-delta_eta/delta_int))*eps_eta
f_bl_inf=f_bl_inf+(1-exp(-delta_eta/delta_int))*eps_eta
f_bl_inf=f_bl_inf/1.0002014996204402_mytype/1.0000000359138641_mytype !To assure 1.0 in infinity

#ifdef DEBG
Expand All @@ -292,7 +291,7 @@ subroutine blasius()
+6.94303207046209_mytype*eta_bl**2 -0.00209716712558639_mytype*eta_bl**1


g_bl_inf=g_bl_inf+(1-exp_prec(-delta_eta/delta_int))*eps_eta
g_bl_inf=g_bl_inf+(1-exp(-delta_eta/delta_int))*eps_eta
g_bl_inf=g_bl_inf/1.000546554_mytype
#ifdef DEBG
if (nrank == 0) write(*,*)'g_bl_inf ', g_bl_inf
Expand Down
Loading