Skip to content

Commit

Permalink
Merge pull request #53 from ucns3d-team/ens
Browse files Browse the repository at this point in the history 
enstrophy-stat
  • Loading branch information
@TakisCFD
TakisCFD authored Jun 2, 2023
2 parents bec3e1f + 1a76a4a commit 49d6fe6
Show file tree
Hide file tree
Showing 7 changed files with 242 additions and 27 deletions.
2 changes: 1 addition & 1 deletion src/declarations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
118 changes: 118 additions & 0 deletions src/flow_operations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
10 changes: 10 additions & 0 deletions src/gradients.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/grid_p.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
4 changes: 2 additions & 2 deletions src/parameters.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
39 changes: 23 additions & 16 deletions src/reconstruct_p.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -701,7 +701,7 @@ SUBROUTINE WENOWEIGHTS(N)

if (ees.eq.5)then

power=2
power=4



Expand Down Expand Up @@ -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

Expand All @@ -1497,7 +1497,7 @@ SUBROUTINE WENOWEIGHTS(N)

if (ees.eq.5)then

power=2
power=4



Expand Down Expand Up @@ -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

Expand All @@ -2455,7 +2455,7 @@ SUBROUTINE WENOWEIGHTS2d(N)

if (ees.eq.5)then

power=2
power=4


else
Expand Down Expand Up @@ -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

Expand All @@ -3225,7 +3225,7 @@ SUBROUTINE WENOWEIGHTS2d(N)

if (ees.eq.5)then

power=2
power=4



Expand Down Expand Up @@ -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
!
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 49d6fe6

Please sign in to comment.