diff --git a/src/declarations.f90 b/src/declarations.f90 index 6f9f535..5626193 100644 --- a/src/declarations.f90 +++ b/src/declarations.f90 @@ -193,7 +193,7 @@ MODULE DECLARATION REAL::X1,ax,ay,az,nx,ny,nz,nnx,nny,nnz !X1 COORDINATE FOR BASIS FUNCTION COMPUTATION OF POLYNOMIALS REAL::Y1 !Y1 COORDINATE FOR BASIS FUNCTION COMPUTATION OF POLYNOMIALS REAL::Z1,x1c,y1c,h1c,z1c !Z1 COORDINATE FOR BASIS FUNCTION COMPUTATION OF POLYNOMIALS -REAL::TAYLOR !ONLY TO BE USED FOR TAYLOR GREEN VORTEX +REAL::TAYLOR,taylor_ens,TAYLOR_ENSX !ONLY TO BE USED FOR TAYLOR GREEN VORTEX REAL::VOLL !TOTAL VOLUME OF THE DOMAIN REAL::UPTURBLIMIT !UPPER TURBULENCE VISCOSITY RATIO REAL::HYBRIDIST !UPPER TURBULENCE VISCOSITY RATIO diff --git a/src/flow_operations.f90 b/src/flow_operations.f90 index 1f0506b..b8c2bc3 100644 --- a/src/flow_operations.f90 +++ b/src/flow_operations.f90 @@ -1156,6 +1156,18 @@ FUNCTION INFLOW2d(INITCOND,POX,POY) +IF (INITCOND.EQ.790)THEN +R=(2.4D0*6**2)/((0.4*6**2)+2) +U=(6*SQRT(1.4))*(70/(2.4*36)) +V=0.0D0 +P=(2.8*36-0.4)/(2.4) + + +END IF + + + + !KINETIC ENERGY FIRST! SKIN=(OO2)*((U**2)+(V**2)) !INTERNAL ENERGY @@ -1949,6 +1961,112 @@ SUBROUTINE VORTEXCALC(N) END SUBROUTINE VORTEXCALC + + +SUBROUTINE ENSTROPHY_CALC(N) +!> @brief +!> This subroutine computes the q-criterion + +IMPLICIT NONE +INTEGER,INTENT(IN)::N +INTEGER::KMAXE,I,IHGT,IHGJ +REAL::SNORM,ONORM +REAL,DIMENSION(3,3)::TVORT,SVORT,OVORT +real,dimension(3,3)::taul,taur,TAU +REAL,DIMENSION(3)::Q,NNN,nall +REAL::UX,UY,UZ,VX,VY,VZ,WX,WY,WZ,RHO12,U12,V12,W12 ,damp,vdamp,TEMPXX + + KMAXE=XMPIELRANK(N) +!$OMP DO SCHEDULE (STATIC) +DO I=1,KMAXE + + VORTET1(1:3,1:3)=ILOCAL_RECON3(I)%GRADS(1:3,1:3) + + DO IHGT=1,3; DO IHGJ=1,3 + TVORT(IHGT,IHGJ)=VORTET1(IHGJ,IHGT) + END DO; END DO + + OVORT=(VORTET1-TVORT) + ONORM=((OVORT(1,1)*OVORT(1,1))+(OVORT(1,2)*OVORT(1,2))+(OVORT(1,3)*OVORT(1,3))+& +(OVORT(2,1)*OVORT(2,1))+(OVORT(2,2)*OVORT(2,2))+(OVORT(2,3)*OVORT(2,3))+(OVORT(3,1)*OVORT(3,1))+& +(OVORT(3,2)*OVORT(3,2))+(OVORT(3,3)*OVORT(3,3))) + + if(boundtype.eq.1)then + + IELEM(N,I)%VORTEX(2)=(0.5D0*(ONORM*u_c(i)%val(1,1))) + + else + + + LEFTV(1:NOF_vARIABLES)=U_C(I)%VAL(1,1:NOF_vARIABLES) + CALL CONS2PRIM(N) + RIGHTV(1:NOF_vARIABLES)=LEFTV(1:NOF_vARIABLES) + CALL SUTHERLAND(N,LEFTV,RIGHTV) + + + + + UX = ILOCAL_RECON3(I)%GRADS(1,1); UY = ILOCAL_RECON3(I)%GRADS(1,2); UZ = ILOCAL_RECON3(I)%GRADS(1,3); + VX = ILOCAL_RECON3(I)%GRADS(2,1); VY = ILOCAL_RECON3(I)%GRADS(2,2); VZ = ILOCAL_RECON3(I)%GRADS(2,3); + WX = ILOCAL_RECON3(I)%GRADS(3,1); WY = ILOCAL_RECON3(I)%GRADS(3,2); WZ = ILOCAL_RECON3(I)%GRADS(3,3); + + + + + + + ! TAU_XX + TAUL(1,1) = (4.0D0/3.0D0)*UX - (2.0D0/3.0D0)*VY - (2.0D0/3.0D0)*WZ + ! TAU_YY + TAUL(2,2) = (4.0D0/3.0D0)*VY - (2.0D0/3.0D0)*UX - (2.0D0/3.0D0)*WZ + ! TAU_ZZ + TAUL(3,3) = (4.0D0/3.0D0)*WZ - (2.0D0/3.0D0)*UX - (2.0D0/3.0D0)*VY + + ! tau_xy + TAUL(1,2) = (UY + VX);TAUL(2,1) = TAUL(1,2) + + ! TAU_XZ + TAUL(1,3) = (WX + UZ);TAUL(3,1) = TAUL(1,3) + + ! TAU_YZ + TAUL(2,3) = (VZ + WY);TAUL(3,2) = TAUL(2,3) + + + IELEM(N,I)%VORTEX(3)=(4.0/3.0)*VISCL(1)*((ux+vy+wz)**2)*ielem(n,i)%TOTVOLUME + + + SNORM=2.0D0*(UX*UX)+2.0D0*(VY*VY)+2.0D0*(WZ*WZ)+(VX+UY)**2+(WZ+VZ)**2+(UZ+WX)**2-((2.0/3.0)*(UX+VY+WZ)**2) + + + ielem(n,i)%vortex(2)=ielem(n,i)%TOTVOLUME*SNORM*viscl(1) +! IELEM(N,I)%VORTEX(2)=VISCL(1)*sNORM!(-ILOCAL_RECON3(I)%GRADS(2,3)+ILOCAL_RECON3(I)%GRADS(3,2)+ILOCAL_RECON3(I)%GRADS(2,1)-ILOCAL_RECON3(I)%GRADS(3,1)+ILOCAL_RECON3(I)%GRADS(1,3)-ILOCAL_RECON3(I)%GRADS(1,2))**2 + + + + + + + end if + +END DO + + +! +! +! DUDX DVDY DWDZ + + + + + + +!$OMP END DO + + + +END SUBROUTINE ENSTROPHY_CALC + + SUBROUTINE VORTEXCALC2D(N) !> @brief !> This subroutine computes the q criterion for 2D diff --git a/src/gradients.f90 b/src/gradients.f90 index 230a545..33991d9 100644 --- a/src/gradients.f90 +++ b/src/gradients.f90 @@ -47,6 +47,10 @@ SUBROUTINE ALLGRADS_INNER(N,ICONSIDERED) CALL COMPUTE_GRADIENTS_MEAN_LSQ(N,ICONSIDERED,NUMBER_OF_DOG,NUMBER_OF_NEI) + IF (INITCOND.EQ.95)THEN + CALL COMPUTE_GRADIENTS_INNER_MEAN_GGS_VISCOUS(N,ICONSIDERED,NUMBER_OF_DOG,NUMBER_OF_NEI) + END IF + CASE(4) CALL COMPUTE_GRADIENTS_MEAN_LSQ(N,ICONSIDERED,NUMBER_OF_DOG,NUMBER_OF_NEI) @@ -188,6 +192,12 @@ SUBROUTINE ALLGRADS_MIX(N,ICONSIDERED) CALL COMPUTE_GRADIENTS_MEAN_LSQ(N,ICONSIDERED,NUMBER_OF_DOG,NUMBER_OF_NEI) + + + IF (INITCOND.EQ.95)THEN + CALL COMPUTE_GRADIENTS_MIX_MEAN_GGS_VISCOUS(N,ICONSIDERED,NUMBER_OF_DOG,NUMBER_OF_NEI) + END IF + CASE(4) CALL COMPUTE_GRADIENTS_MEAN_LSQ(N,ICONSIDERED,NUMBER_OF_DOG,NUMBER_OF_NEI) diff --git a/src/grid_p.f90 b/src/grid_p.f90 index 10a1120..df22f39 100644 --- a/src/grid_p.f90 +++ b/src/grid_p.f90 @@ -1835,7 +1835,7 @@ SUBROUTINE NEIGHBOURSS(N,IELEM,IMAXE,IMAXN,XMPIE,XMPIN,XMPIELRANK,RESTART,INODEr ! if (itestcase.eq.4)then - allocate(ielem(n,i)%vortex(1)); ielem(n,i)%vortex=zero + allocate(ielem(n,i)%vortex(3)); ielem(n,i)%vortex=zero ALLOCATE(IELEM(N,I)%AVARS(NOF_VARIABLES)); ! end if diff --git a/src/parameters.f90 b/src/parameters.f90 index b50d01b..697fc85 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -387,7 +387,7 @@ SUBROUTINE READ_UCNS3D CFLRAMP=0 !CFL RAMPING: |0: DEACTIVATED |1:ACTIVATED emetis=6 !Metis partitioner : 1: Hybrid metis, 2:adaptive weights for hybrid grids, 3: Uniform metis partionioner,4:NODAL,6=PARMETS itold=10000 !TOLERANCE=n_iterations - GRIDAR1=20.0 ! 0 5.0 7.0 LIMIT ASPECT RATIO CELLS, + GRIDAR1=12.0 ! 0 5.0 7.0 LIMIT ASPECT RATIO CELLS, GRIDAR2=20.0 ! LIMIT VOLUME CELLS fastest=0 ! 0 ||Fastest, no coordinate mapping (1: engaged,0:with transformation) lmach_style=0 !0 ||LOW MACH TREATMENT (1 ACTIVATE, 0 DISABLE),lmach_style(0=only normal component,1=all components) @@ -425,7 +425,7 @@ SUBROUTINE READ_UCNS3D CFLRAMP=0 !CFL RAMPING: |0: DEACTIVATED |1:ACTIVATED emetis=6 !Metis partitioner : 1: Hybrid metis, 2:adaptive weights for hybrid grids, 3: Uniform metis partionioner,4:NODAL,6=PARMETS itold=10000 !TOLERANCE=n_iterations - GRIDAR1=8.0 ! 0 5.0 7.0 LIMIT ASPECT RATIO CELLS, + GRIDAR1=20.0 ! 0 5.0 7.0 LIMIT ASPECT RATIO CELLS, GRIDAR2=20.0 ! LIMIT VOLUME CELLS fastest=0 ! 0 ||Fastest, no coordinate mapping (1: engaged,0:with transformation) lmach_style=0 !0 ||LOW MACH TREATMENT (1 ACTIVATE, 0 DISABLE),lmach_style(0=only normal component,1=all components) diff --git a/src/reconstruct_p.f90 b/src/reconstruct_p.f90 index b0f60ba..dfbf5b6 100644 --- a/src/reconstruct_p.f90 +++ b/src/reconstruct_p.f90 @@ -690,7 +690,7 @@ SUBROUTINE WENOWEIGHTS(N) if (poly.eq.4)then - divbyzero=ielem(n,iconsidered)%totvolume**3 + divbyzero=ielem(n,iconsidered)%totvolume**2 else @@ -701,7 +701,7 @@ SUBROUTINE WENOWEIGHTS(N) if (ees.eq.5)then - power=2 + power=4 @@ -1486,7 +1486,7 @@ SUBROUTINE WENOWEIGHTS(N) if (poly.eq.4)then - divbyzero=ielem(n,iconsidered)%totvolume**3 + divbyzero=ielem(n,iconsidered)%totvolume**2 else @@ -1497,7 +1497,7 @@ SUBROUTINE WENOWEIGHTS(N) if (ees.eq.5)then - power=2 + power=4 @@ -2444,7 +2444,7 @@ SUBROUTINE WENOWEIGHTS2d(N) IF ((IELEM(N,I)%FULL.EQ.1).AND.(ielem(n,i)%TROUBLED.EQ.1))THEN if (poly.eq.4)then - divbyzero=ielem(n,i)%totvolume**3 + divbyzero=ielem(n,iconsidered)%totvolume**2 else @@ -2455,7 +2455,7 @@ SUBROUTINE WENOWEIGHTS2d(N) if (ees.eq.5)then - power=2 + power=4 else @@ -3214,7 +3214,7 @@ SUBROUTINE WENOWEIGHTS2d(N) IF((IELEM(N,I)%FULL.EQ.1).AND.(ielem(n,i)%TROUBLED.EQ.1))THEN if (poly.eq.4)then - divbyzero=ielem(n,i)%totvolume**3 + divbyzero=ielem(n,iconsidered)%totvolume**2 else @@ -3225,7 +3225,7 @@ SUBROUTINE WENOWEIGHTS2d(N) if (ees.eq.5)then - power=2 + power=4 @@ -9763,7 +9763,7 @@ SUBROUTINE APPLY_FILTER_DG(n) !end if - if (ielem(n,i)%er.GT.1.1)THen + if (ielem(n,i)%er.GT.(1.15))THen ielem(n,I)%FILTERED=1 ! @@ -10241,16 +10241,20 @@ SUBROUTINE APPLY_ADDA_FILTER(N) IF (ielem(n,i)%er.gt.1.3)THEN - LWCX1=2!INCREASE DISSIPATION + LWCX1=10!INCREASE DISSIPATION end if + ! END IF ! - IF (ielem(n,i)%er.LE.0.9)THEN + IF (ielem(n,i)%er.LE.0.95)THEN + + LWCX1=1000 + + !LWCX1=10**(12-(4*ielem(n,i)%er**0.8)) - LWCX1=10**(2-(1*ielem(n,i)%er**0.8)) END IF @@ -10276,16 +10280,19 @@ SUBROUTINE APPLY_ADDA_FILTER(N) IF (ielem(n,i)%er.gt.1.3)THEN - LWCX1=2!INCREASE DISSIPATION + LWCX1=10!INCREASE DISSIPATION end if ! END IF ! - IF (ielem(n,i)%er.LE.0.9)THEN + IF (ielem(n,i)%er.LE.0.95)THEN + + - LWCX1=10**(5-(2*ielem(n,i)%er**0.8)) + LWCX1=1000 + !LWCX1=100**(12-(4*ielem(n,i)%er**0.8)) END IF @@ -10321,7 +10328,7 @@ SUBROUTINE FIX_DISSIPATION(N) if (ielem(n,i)%full.eq.1)then !1)reduce dissipation - IF (IELEM(N,I)%LWCX2.GT.5)THEN + IF (IELEM(N,I)%LWCX2.GT.10)THEN IF (IELEM(N,I)%WCX(1).GE.0.999)THEN IELEM(N,I)%DISS=max(IELEM(N,I)%DISS-0.1,0.4d0) !reduce dissipation even more Else diff --git a/src/time_p.f90 b/src/time_p.f90 index c7f846d..ffaf272 100644 --- a/src/time_p.f90 +++ b/src/time_p.f90 @@ -2375,6 +2375,10 @@ SUBROUTINE CALL_FLUX_SUBROUTINES_3D CALL CALCULATE_FLUXESHI(N) CASE(3) CALL CALCULATE_FLUXESHI_CONVECTIVE(N) + + + + if ((SOURCE_ACTIVE.EQ.1))then call SOURCES_COMPUTATION_ROT(N) end if @@ -2394,7 +2398,10 @@ SUBROUTINE CALL_FLUX_SUBROUTINES_3D - + IF (INITCOND.EQ.95)THEN + + CALL ENSTROPHY_CALC(N) + END IF IF (statistics.eq.1)THEN @@ -3961,7 +3968,7 @@ SUBROUTINE TIME_MARCHING(N) INTEGER,INTENT(IN)::N real,dimension(1:5)::DUMMYOUT,DUMMYIN INTEGER::I,KMAXE,TTIME -REAL::CPUT1,CPUT2,CPUT3,CPUT4,CPUT5,CPUT6,CPUT8,timec3,TIMEC1,TIMEC4,TIMEC8,TOTV1,TOTV2,DUMEtg1,DUMEtg2,TOTK,TZX1,TZX2,resolx +REAL::CPUT1,CPUT2,CPUT3,CPUT4,CPUT5,CPUT6,CPUT8,timec3,TIMEC1,TIMEC4,TIMEC8,TOTV1,TOTV2,DUMEtg1,DUMEtg2,TOTK,TZX1,TZX2,resolx,totens,totens1,totens2,totensx,totensx1,totensx2 kill=0 T=RES_TIME resolx=0.01 @@ -4011,6 +4018,16 @@ SUBROUTINE TIME_MARCHING(N) !$OMP END MASTER !$OMP BARRIER + + + + if ((it.eq.0).and.(initcond.eq.95))then + call EXCHANGE_HIGHER(N) + call ARBITRARY_ORDER3(N) + call ENSTROPHY_CALC(N) + + + end if DO @@ -4058,13 +4075,21 @@ SUBROUTINE TIME_MARCHING(N) IF (INITCOND.eq.95)THEN - TOTK=0 + TOTK=0;TOTENS=0;totensx=0.0d0 DO I=1,xmpielrank(n) TOTK=TOTK+IELEM(N,I)%TOTVOLUME*U_C(I)%VAL(1,1)*(1.0/2.0)*& (((U_C(I)%VAL(1,2)/U_C(I)%VAL(1,1))**2)+((U_C(I)%VAL(1,3)/U_C(I)%VAL(1,1))**2)+((U_C(I)%VAL(1,4)/U_C(I)%VAL(1,1))**2)) + if (BOUNDTYPE.eq.1)then + + TOTENS=TOTENS+(IELEM(N,I)%TOTVOLUME*U_C(I)%VAL(1,1)*(1.0/2.0)*& + IELEM(N,I)%VORTEX(2)) + else + TOTENS=TOTENS+(IELEM(N,I)%VORTEX(2)) + TOTENSx=TOTENSx+(IELEM(N,I)%VORTEX(3)) + end if END DO ! @@ -4073,10 +4098,24 @@ SUBROUTINE TIME_MARCHING(N) CALL MPI_BARRIER(MPI_COMM_WORLD,IERROR) CALL MPI_ALLREDUCE(DUMEtg1,DUMEtg2,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERROR) TOTK=DUMEtg2 + DUMEtg1=TOTENS + DUMEtg2=0.0 + CALL MPI_BARRIER(MPI_COMM_WORLD,IERROR) + CALL MPI_ALLREDUCE(DUMEtg1,DUMEtg2,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERROR) + TOTENS=DUMEtg2 + DUMEtg1=TOTENSx + DUMEtg2=0.0 + CALL MPI_BARRIER(MPI_COMM_WORLD,IERROR) + CALL MPI_ALLREDUCE(DUMEtg1,DUMEtg2,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERROR) + TOTENSx=DUMEtg2 IF (N.EQ.0)THEN TOTV1=TOTK/((2.0*PI)**3) + TOTENS1=TOTENS/(((2.0*PI)**3)) + TOTENSx1=TOTENSx/((2.0*PI)**3) IF (it.eq.0)THEN TAYLOR=TOTK + TAYLOR_ENS=TOTENS + TAYLOR_ENSx=TOTENSx END IF @@ -4188,14 +4227,29 @@ SUBROUTINE TIME_MARCHING(N) IF (INITCOND.eq.95)THEN - TOTK=0 + TOTK=0; TOTENS=0.0; totensx=0.0d0 DO I=1,xmpielrank(n) TOTK=TOTK+IELEM(N,I)%TOTVOLUME*U_C(I)%VAL(1,1)*(1.0/2.0)*& (((U_C(I)%VAL(1,2)/U_C(I)%VAL(1,1))**2)+((U_C(I)%VAL(1,3)/U_C(I)%VAL(1,1))**2)+((U_C(I)%VAL(1,4)/U_C(I)%VAL(1,1))**2)) - + + + + + + if (BOUNDTYPE.eq.1)then + + TOTENS=TOTENS+(IELEM(N,I)%TOTVOLUME*U_C(I)%VAL(1,1)*(1.0/2.0)*& + IELEM(N,I)%VORTEX(2)) + else + + TOTENS=TOTENS+(IELEM(N,I)%VORTEX(2)) + TOTENSx=TOTENSx+(IELEM(N,I)%VORTEX(3)) + end if + + @@ -4206,10 +4260,29 @@ SUBROUTINE TIME_MARCHING(N) CALL MPI_BARRIER(MPI_COMM_WORLD,IERROR) CALL MPI_ALLREDUCE(DUMEtg1,DUMEtg2,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERROR) TOTK=DUMEtg2 + + DUMEtg1=TOTENS + DUMEtg2=0.0 + CALL MPI_BARRIER(MPI_COMM_WORLD,IERROR) + CALL MPI_ALLREDUCE(DUMEtg1,DUMEtg2,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERROR) + TOTENS=DUMEtg2 + + + + DUMEtg1=TOTENSx + DUMEtg2=0.0 + CALL MPI_BARRIER(MPI_COMM_WORLD,IERROR) + CALL MPI_ALLREDUCE(DUMEtg1,DUMEtg2,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERROR) + TOTENSx=DUMEtg2 + IF (N.EQ.0)THEN TOTV2=TOTK/((2.0*PI)**3) + TOTENS2=TOTENS/((2.0*PI)**3) + TOTENSx2=TOTENSx/((2.0*PI)**3) IF (it.eq.0)THEN TAYLOR=TOTK + TAYLOR_ENS=TOTENS + TAYLOR_ENSx=TOTENSx END IF IF (IT.EQ.0)THEN @@ -4217,8 +4290,15 @@ SUBROUTINE TIME_MARCHING(N) ELSE OPEN(73,FILE='ENERGY.dat',FORM='FORMATTED',STATUS='old',ACTION='WRITE',POSITION='APPEND') END IF - - WRITE(73,*)T,TOTK/TAYLOR,-(TOTV2-TOTV1)/DT + IF (DG.EQ.1)THEN + WRITE(73,'(E14.7,1X,E14.7,1X,E14.7)')T,TOTK/TAYLOR,-(TOTV2-TOTV1)/DT + ELSE + if (boundtype.eq.1)then + WRITE(73,'(E14.7,1X,E14.7,1X,E14.7,1X,E14.7)')T,TOTK/TAYLOR,-(TOTV2-TOTV1)/DT,TOTENS/TAYLOR_ENS + else + WRITE(73,'(E14.7,1X,E14.7,1X,E14.7,1X,E14.7,1X,E14.7)')T,TOTK/TAYLOR,-(TOTV2-TOTV1)/DT,TOTENS,TOTENSx + end if + END IF CLOSE(73) END IF