PROGRAM BISPEC1_MOVIE ! ! ! This program is a test driver for a set of spectral ! solvers for the biharmonic equation. ! ! USE INTERFACE, ONLY:leg,cg_leg,cheb,cg_cheb,cheb_leg, & cg_cheb_leg,setup,adjust,funkval, & petterfrhs ! use Bi_Parameter ! IMPLICIT NONE ! ! Important driver parameters: ! ! problem the test problem number selected ! ! method the method to be used ! ! 1 Direct Legendre Solver ! 2 Direct Chebychev Solver ! 3 Direct Cheychev-Legendre Solver ! 4 Iterative Legendre Solver ! 5 Iterative Chebychev Solver ! 6 Iterative Chevbychev-Legendre Solver ! ! n the problem size. The problem is discretized with ! 2*n+4 points in each coordinate direction. ! INTEGER problem INTEGER method INTEGER n ! ! Note that the parameters alpha and beta in the equation ! can be changed in the module Parameter for this test ! program. For uniqueness purposes these are called ! pde_alpha and pde_beta. ! LOGICAL legtrans,update INTEGER i,ii,flag,j INTEGER i_save,ii_save,nstop,e_i,step,count,t DOUBLE PRECISION alpha,beta,k,x,y,q0,epsilon,dt,z DOUBLE PRECISION mean_psi,tstop,w,dx,dy,f,del_psi ! INTEGER count,count_max,count_rate,ocount DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: hoyre DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: loes1 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: coeff DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: loes2 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: psi DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: hoyreut DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xjac DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xjacx DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xjacy DOUBLE PRECISION max1,max2,diff,gettime,etime DOUBLE PRECISION t1,t2 DOUBLE PRECISION tid1,tid2,tid3,tid4,tid5,tid6,tid7,tid8 DOUBLE PRECISION au,bu,cu,du,x1,x2,xverdi,yverdi,fverdi,x_unif,y_unif !-----------------------------------------------------------------: ! Input parameters !-----------------------------------------------------------------: problem = 1 method = 1 WRITE(*,'(A)',ADVANCE = 'NO') 'Enter problem size: ' n = 12 ! read(5,*)n WRITE(*,*)'n=',n WRITE(*,*) 'The problem is solved on the area [a,b]x[c,d]' ! read(5,*) au,bu,cu,du au = 0 bu = 16*pi cu = 0 du = 4*pi WRITE(*,*)'a,b,c,d',au,bu,cu,du !-----------------------------------------------------------------: ! Time parameters !-----------------------------------------------------------------: x1 = 4.d0/((du-cu)**2) x2 = 4.d0/((bu-au)**2) dt = 0.1 tstop = 3000 nstop = tstop/dt WRITE(*,*)'nstop = ',nstop step=50 count = 0 !-----------------------------------------------------------------: ! Open output file !-----------------------------------------------------------------: WRITE(*,*)'Writing outputfile.' OPEN(unit=10,file='bis_movie_output') OPEN(unit=15,file='bis_epsilon_output') OPEN(unit=20,file='bis_f_psi_output') WRITE(10,*)au WRITE(10,*)bu WRITE(10,*)cu WRITE(10,*)du WRITE(10,*)n*2+4 WRITE(10,*)n*2+4 !-----------------------------------------------------------------: ! Write the number of times the matlab-program shall read the ! solution from file. !-----------------------------------------------------------------: WRITE(10,*)nstop/step !-----------------------------------------------------------------: ! Allocate memory !-----------------------------------------------------------------: ALLOCATE (hoyre(0:2*N+3,0:2*N+3)) ALLOCATE (xjac(0:2*N+3)) ALLOCATE (xjacx(0:2*N+3)) ALLOCATE (xjacy(0:2*N+3)) ALLOCATE (loes1(0:2*N+3,0:2*N+3)) ALLOCATE (coeff(0:2*N+3,0:2*N+3)) ALLOCATE (loes2(0:2*N+3,0:2*N+3)) ALLOCATE (psi(0:10*N+3,0:10*N+3)) ALLOCATE (hoyreut(0:(2*N+3)**2)) !-----------------------------------------------------------------: ! Legendre !-----------------------------------------------------------------: IF(method.EQ.1.OR.method.EQ.4) THEN legtrans = .TRUE. ELSE legtrans = .FALSE. END IF CALL legendrel(2*N+3,xjac) !-----------------------------------------------------------------: ! Transformerer til [a,b] i x-retning og [c,d] i y-retn. !-----------------------------------------------------------------: DO i = 0,2*N+3 xjacx(i) = 0.5d0*((bu+au)+(bu-au)*xjac(i)) xjacy(i) = 0.5d0*((du+cu)+(du-cu)*xjac(i)) END DO !-----------------------------------------------------------------: ! Write Non-uniform coords !-----------------------------------------------------------------: DO i = 0,2*N+3 WRITE(10,*) xjacx(i) END DO DO i = 0,2*N+3 WRITE(10,*) xjacy(i) END DO !-----------------------------------------------------------------: ! The integral-error is of order dx^2 !-----------------------------------------------------------------: dx=(bu-au)/(10*N+3) dy=(du-cu)/(10*N+3) !-----------------------------------------------------------------: ! Epsilon loop !-----------------------------------------------------------------: DO e_i=1,25 flag = 0 q0 = 1 epsilon = e_i*10d0**(-2) pde_alpha = -2*q0 pde_beta = q0**2 -epsilon+1/dt alpha = pde_alpha beta = pde_beta !-----------------------------------------------------------------: ! Initial state !-----------------------------------------------------------------: DO i=0,2*n+3 DO j=0,2*n+3 loes1(i,j) = SIN(xjacy(j)) END DO END DO update = .false. !-----------------------------------------------------------------: ! Tau Loop !-----------------------------------------------------------------: DO t=0,nstop !-------------------------------------------------------------: !---- Update right hand side !-------------------------------------------------------------: ! IF (MOD(t,step) .EQ. 0) THEN ! WRITE(*,*)(t*dt) ! end if IF(update) THEN DO j = 0,2*N+3 DO i = 0,2*N+3 loes1(i,j) = loes1(i,j)/dt-loes1(i,j)**3 END DO END DO END IF update = .TRUE. method = 1 etime = gettime() CALL leg(loes1,coeff,xjac,alpha,beta,flag,N,x1,x2,tid1,tid2,tid3,tid4) etime = gettime()-etime flag=1 !-----------------------------------------------------------------: ! End Tau Loop, only for epsilon !-----------------------------------------------------------------: ENDDO !-----------------------------------------------------------------: ! Uniform coords on [-1 1]x[-1 1] ! (Funkval can only take values from the transformed region) !-----------------------------------------------------------------: !-----------------------------------------------------------------: ! Write mean_psi, F and the solution to file for every teller ! number of timesteps !-----------------------------------------------------------------: ! IF (MOD(t,step) .EQ. 0 .and. t .ne. 0) count = count +1 DO j=0,10*N+3 DO i=0,10*N+3 x_unif = -1+2.d0/(10*N+3)*i y_unif = -1+2.d0/(10*N+3)*j !-----------------------------------------------------------------: ! Can not use funkval on the boundary !-----------------------------------------------------------------: IF (x_unif .EQ. -1 .OR. x_unif .EQ. 1 .OR. y_unif& .EQ. -1 .OR. y_unif .EQ. 1 ) THEN psi(i,j) = 0 ELSE psi(i,j) = funkval(x_unif,y_unif,coeff,2*N+4) END IF ENDDO ENDDO mean_psi = 0 f = 0 DO j=0,10*N+3 DO i=0,10*N+3 !-----------------------------------------------------------------: ! Mean value for psi^2 !-----------------------------------------------------------------: IF(j .LT. 10*N+3 .AND. i .LT. 10*N+3) THEN w = (psi(i,j)+psi(i+1,j)+psi(i+1,j+1)+psi(i,j+1))/4 mean_psi = mean_psi + w*w*dx*dy !-----------------------------------------------------------------: ! Energy-function ! The Laplace operator is approximated using a second order ! finite difference scheme based on Taylor series. !-----------------------------------------------------------------: IF(j .NE. 0 .AND. i .NE. 0) THEN del_psi=(psi(i+1,j)-2*psi(i,j)+psi(i-1,j))/(dx& **2)+(psi(i,j+1)-2*psi(i,j)+psi(i,j-1))/(dy& **2) f= f+0.5*(0.5*psi(i,j)**4-epsilon*psi(i,j)& **2+(del_psi+psi(i,j))**2)*dx*dy END IF END IF ENDDO ENDDO mean_psi = mean_psi/((bu-au)*(du-cu)) WRITE(20,*)t*dt,mean_psi,f DO j = 0,2*N+3 DO i = 0,2*N+3 WRITE(10,*)loes1(i,j) END DO END DO ! END IF WRITE(*,*)epsilon,mean_psi WRITE(15,*)epsilon,mean_psi END DO ! end do DEALLOCATE(hoyre,xjac,xjacx,xjacy,loes1,coeff,loes2,hoyreut,psi) CLOSE(10) CLOSE(15) CLOSE(20) END PROGRAM bispec1_MOVIE