From bd912eebc32126e9452f17717b9e1f6e0b86f54b Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Tue, 25 Jan 2022 15:40:35 +0000 Subject: [PATCH] molecular and hyperviscosity --- src/makefile | 1 - src/parameters.f90 | 27 +++++++++++++--------- src/spectral.f90 | 56 ++++++++++++++++++++++++++++++++++------------ src/strat.f90 | 48 +++++++++++++++++++++++---------------- 4 files changed, 88 insertions(+), 44 deletions(-) diff --git a/src/makefile b/src/makefile index 12cccb6..8cabe10 100644 --- a/src/makefile +++ b/src/makefile @@ -35,7 +35,6 @@ timestep: $(objects) $(sourcedir)/timestep.f90 bminmax: $(objects) $(sourcedir)/bminmax.f90 $(f90) parameters.o constants.o $(sourcedir)/bminmax.f90 -o bminmax $(flags) - $(objects): $(fft_lib) $(sources) $(f90) $(fft_lib) $(sources) -c $(flags) diff --git a/src/parameters.f90 b/src/parameters.f90 index 657f35e..f495803 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -4,20 +4,27 @@ module parameters ! the suite of ps f90 files. !Domain grid dimensions: -integer,parameter:: nx=256,ny=384 +integer,parameter:: nx=N_X,ny=N_Y !Domain width in x and limits in y: -double precision,parameter:: ellx=1000d0 -double precision,parameter:: ymin=0.d0,ymax=1500d0 +double precision,parameter:: ellx=L_X +double precision,parameter:: ymin=Y_MIN,ymax=Y_MAX !Simulation duration and data save interval: -double precision,parameter:: tsim=900.d0,tgsave=10.d0 +double precision,parameter:: tsim=T_SIM,tgsave=T_GSAVE - !Hyperviscosity parameters: -integer,parameter:: nnu=3 -double precision,parameter:: prediss=10.d0 -! Damping rate is prediss*zeta_rms*(k/k_max)^(2*nnu) on wavenumber -! k where k_max is the maximum x or y wavenumber and zeta_rms is -! the rms vorticity. Note: nnu = 3 and prediss = 10 are recommended. + !(Hyper)viscosity parameters: +integer,parameter:: nnu=N_NU +double precision,parameter:: prediss=PRE_DISS +! If nnu = 1, this is the molecular viscosity case. Then, we +! choose the viscosity nu = prediss*((b_max-b_min)/k_{x,max}^3) +! where k_{x_max} is the maximum x wavenumber. +! Note: prediss = 2 is recommended. +! ---------------------------------------------------------------- +! If nnu > 1, this is the hyperviscosity case. Then, the damping +! rate is prediss*zeta_char*(k/k_max)^(2*nnu) on wavenumber k +! where k_max is the maximum x or y wavenumber and zeta_char is +! a characteristic vorticity (see subroutine adapt of strat.f90). +! Note: nnu = 3 and prediss = 10 are recommended. end module diff --git a/src/spectral.f90 b/src/spectral.f90 index 3847309..e46dcdf 100644 --- a/src/spectral.f90 +++ b/src/spectral.f90 @@ -17,7 +17,7 @@ module spectral integer:: kmag(0:nxm1,0:ny),kmax !====================================================================! -! From main code: call init_spectral to initialise ! +! From main code: call init_spectral(bbdif) to initialise ! ! then call main_invert(zz,uu,vv) to perform inversion ! !====================================================================! @@ -25,11 +25,14 @@ module spectral !===================================================================== -subroutine init_spectral +subroutine init_spectral(bbdif) !Declarations: implicit none + !Passed variable (bbdif = max(b) - min(b) at t = 0): +double precision:: bbdif + !Local variables: !Maximum x wavenumber: @@ -154,20 +157,45 @@ subroutine init_spectral enddo !--------------------------------------------------------------------- - !Define dissipation operator: -visc=prediss/max(rkxmax,rkymax)**(2*nnu) -hdis(0,0)=zero -do kx=1,nxm1 - hdis(kx,0)=visc*rkx(kx)**(2*nnu) -enddo -do ky=1,ny - hdis(0,ky)=visc*rky(ky)**(2*nnu) -enddo -do ky=1,ny +! Damping, viscous or hyperviscous: +if (nnu .eq. 1) then + !Define viscosity: + visc=prediss*sqrt(bbdif/rkxmax**3) + write(*,'(a,1p,e14.7)') ' Viscosity nu = ',visc + + !Define spectral dissipation operator: + hdis(0,0)=zero do kx=1,nxm1 - hdis(kx,ky)=visc*(rkx(kx)**2+rky(ky)**2)**nnu + hdis(kx,0)=visc*rkx(kx)**2 enddo -enddo + do ky=1,ny + hdis(0,ky)=visc*rky(ky)**2 + enddo + do ky=1,ny + do kx=1,nxm1 + hdis(kx,ky)=visc*(rkx(kx)**2+rky(ky)**2) + enddo + enddo + +else + !Define hyperviscosity: + visc=prediss/max(rkxmax,rkymax)**(2*nnu) + write(*,'(a,1p,e14.7)') ' Hyperviscosity nu = ',visc + + !Define dissipation operator: + hdis(0,0)=zero + do kx=1,nxm1 + hdis(kx,0)=visc*rkx(kx)**(2*nnu) + enddo + do ky=1,ny + hdis(0,ky)=visc*rky(ky)**(2*nnu) + enddo + do ky=1,ny + do kx=1,nxm1 + hdis(kx,ky)=visc*(rkx(kx)**2+rky(ky)**2)**nnu + enddo + enddo +endif return end subroutine init_spectral diff --git a/src/strat.f90 b/src/strat.f90 index 1fa0f71..bde6521 100644 --- a/src/strat.f90 +++ b/src/strat.f90 @@ -93,31 +93,33 @@ subroutine initialise implicit none - !Local variable: + !Local variables: double precision:: qq(0:ny,0:nxm1) +double precision:: bbdif !---------------------------------------------------------------------- - !Initialise inversion constants and arrays: -call init_spectral - -!---------------------------------------------------------------------- - !Read vorticity and convert to spectral space as zs: -open(11,file='zz_init.r8',form='unformatted', & + !Read buoyancy and find min/max value for use in determining viscosity: +open(11,file='bb_init.r8',form='unformatted', & & access='direct',status='old',recl=2*nbytes) read(11,rec=1) t,qq close(11) -call ptospc_fc(nx,ny,qq,zs,xfactors,yfactors,xtrig,ytrig) +bbdif=maxval(qq)-minval(qq) + + !Initialise inversion constants and arrays: +call init_spectral(bbdif) + + !Convert buoyancy (in qq) to spectral space as bs: +call ptospc_fc(nx,ny,qq,bs,xfactors,yfactors,xtrig,ytrig) !qq is overwritten here. -!---------------------------------------------------------------------- - !Read buoyancy and convert to spectral space as bs: -open(11,file='bb_init.r8',form='unformatted', & + !Read vorticity and convert to spectral space as zs: +open(11,file='zz_init.r8',form='unformatted', & & access='direct',status='old',recl=2*nbytes) read(11,rec=1) t,qq close(11) -call ptospc_fc(nx,ny,qq,bs,xfactors,yfactors,xtrig,ytrig) +call ptospc_fc(nx,ny,qq,zs,xfactors,yfactors,xtrig,ytrig) !qq is overwritten here. !---------------------------------------------------------------------- @@ -383,13 +385,21 @@ subroutine adapt dt4=dt/four !--------------------------------------------------------------------- - !Update hyperdiffusion operator used in time stepping: -dfac=zzch*dt/two - !zzch is the characteristic vorticity defined above. -diss=two/(one+dfac*hdis) - !hdis = C*(K/K_max)^{2p} where K^2 = k_x^2+k_y^2, p is the order, - !K_max is the maximum x or y wavenumber and C is a dimensionless - !prefactor (see spectral.f90 and parameters.f90 where C = prediss). +if (nnu .eq. 1) then + !Update diffusion operator used in time stepping: + dfac=dt/two + diss=two/(one+dfac*hdis) + !hdis = nu*(k_x^2+k_y^2) where nu is the viscosity coefficient + !(see spectral.f90 and parameters.f90). +else + !Update hyperdiffusion operator used in time stepping: + dfac=zzch*dt/two + !zzch is the characteristic vorticity defined above. + diss=two/(one+dfac*hdis) + !hdis = C*(K/K_max)^{2p} where K^2 = k_x^2+k_y^2, p is the order, + !K_max is the maximum x or y wavenumber and C is a dimensionless + !prefactor (see spectral.f90 and parameters.f90 where C = prediss). +endif !--------------------------------------------------------------------- !Save |u|_max, N_max and gamma_max to monitor.asc: