c-------------------------------------------------------------: program biharmonic_movie c-------------------------------------------------------------: c---- c---- u + 2*u + u + alpha*(u +u ) + beta*u = f(x,y) c---- xxxx xxyy yyyy xx yy c---- c-------------------------------------------------------------: implicit none external f1,f2,f3,gettime,approxit integer m,n,wn,idf,nmax,mmax parameter (mmax=999) parameter (nmax=999) parameter (idf=mmax+2) parameter (wn=3*mmax+4*nmax+4*nmax+2*mmax+0.5*(nmax+1)**2+19+100) integer iflag,itcg,lw,i,t,j,problem,nstop,step,count double precision a,b,c,d,alpha,beta,tol double precision bda(nmax) double precision bdb(nmax) double precision bdc(mmax) double precision bdd(mmax) double precision f(mmax+2,nmax+2) double precision psi(5*mmax+2,5*nmax+2) double precision w(wn) double precision x,y,dt,runtime,pi,q0,epsilon,tstop,z double precision mean_psi,f_tau,p,dx,dy,del_psi,e_i double precision f1,f2,f3,gettime,approxit pi = 3.141592653589793 100 format('T = ',F6.1) 101 format(I6,' frames.') c-------------------------------------------------------------: c---- Initialize problem c-------------------------------------------------------------: c---- Problem & constants c-------------------------------------------------------------: problem = 3 WRITE(*,*)'Solving problem' c-------------------------------------------------------------: c---- Problem size & time c-------------------------------------------------------------: m = 101 n = 101 tstop =3000 dt = 0.1 step = 50 count = 0 nstop = tstop/dt c-------------------------------------------------------------: c---- Domain c-------------------------------------------------------------: a = 0 b = 16*pi c = 0 d = 4*pi dx = (b-a)/((m+1)) dy = (d-c)/((n+1)) c-------------------------------------------------------------: c---- Write outut file c-------------------------------------------------------------: write(*,*)'Writing ouptufile.' open(unit=10,file='bih_movie_output') open(unit=15,file='bih_epsilon_output') open(unit=20,file='bih_f_psi_output') c open(unit=30,file='bla1') c open(unit=40,file='bla2') write(10,*)a write(10,*)b write(10,*)c write(10,*)d write(10,*)m+2 write(10,*)n+2 write(10,*)nstop/step c-------------------------------------------------------------: c---- Constant of the equation c-------------------------------------------------------------: DO e_i=1,25 epsilon = e_i*10d0**(-2) c epsilon = 0.25 q0 = 1 alpha = 2*q0 beta = q0-epsilon+1/dt lw = wn c-------------------------------------------------------------: c---- Option iflag c-------------------------------------------------------------: iflag = 4 c-------------------------------------------------------------: c---- Function values, f1 c-------------------------------------------------------------: do i = 2,m+1 do j = 2,n+1 x =a+(i-1)*(b-a)/(m+1) y =c+(j-1)*(d-c)/(n+1) f(i,j) = f1(problem,x,y,alpha,beta,dt) enddo enddo c-------------------------------------------------------------: c---- Boundary values, f2 c-------------------------------------------------------------: do i = 1,m+2 x = a+(i-1)*(b-a)/(m+1) y = c f(i,1) = f2(problem,x,y,dt) x = a+(i-1)*(b-a)/(m+1) y = d f(i,n+2) = f2(problem,x,y,dt) end do do j = 1,n+2 x = a y = c+(j-1)*(d-c)/(n+1) f(1,j) = f2(problem,x,y,dt) y = b x = c+(j-1)*(d-c)/(n+1) f(m+2,j) = f2(problem,x,y,dt) end do c-------------------------------------------------------------: c---- Exterior derivatives values, f3 c-------------------------------------------------------------: do j = 1,n x = a y = c+j*(d-c)/(n+1) bda(j) = f3(problem,x,y,4,dt) x = b y = c+j*(d-c)/(n+1) bdb(j) = f3(problem,x,y,2,dt) enddo do i = 1,m x = a+i*(b-a)/(m+1) y = c bdc(i) = f3(problem,x,y,1,dt) x = a+i*(b-a)/(m+1) y = d bdd(i) = f3(problem,x,y,3,dt) enddo c-------------------------------------------------------------: c---- x- and y- values c-------------------------------------------------------------: c write(40,*)m+2 DO i = 0,n+1 WRITE(10,*)(b-a)*i/(n+1) END DO DO i = 0,m+1 WRITE(10,*)(d-c)*i/(n+1) c WRITE(40,*)(b-a)*i/(n+1) END DO c write(30,*)5*(m+1)+1 c DO i = 0,5*(n+1) c WRITE(30,*)(b-a)*i/(5*(n+1)) c END DO c-------------------------------------------------------------: c---- Write first timestep c-------------------------------------------------------------: c do i = 1,m+2 c do j = 1, n+2 cc Hack to avoid small numbers without 'E' c c z = f(i,j) c if (z .gt. -1d-99 .and. z .lt. 1d-99) z =0 c write(10,*)z c end do c end do c-------------------------------------------------------------: c---- Loop c-------------------------------------------------------------: do t = 0,nstop IF (MOD(t,step) .EQ. 0) THEN write(*,100)(t*dt) end if c-------------------------------------------------------------: c---- Solve problem c-------------------------------------------------------------: runtime = gettime() call dbihar(a,b,m,bda,bdb,bdc,bdd,c,d,n,f,idf,alpha,beta $ ,iflag,tol,itcg,w,lw) runtime = gettime()-runtime IF (MOD(t,step) .EQ. 0 .and. t .ne. 0) count = count +1 do j = 1, n+2 do i = 1,m+2 c if(t .eq. 0) write(40,*)f(i,j) IF (MOD(t,step) .EQ. 0 .and. t .ne. 0) THEN c-------------------------------------------------------------: c---- Write f() to outputfile c-------------------------------------------------------------: c Hack to avoid small numbers without 'E' c z = f(i,j) if (z .gt. -1d-99 .and. z .lt. 1d-99) z =0 c write(10,*)z end if c-------------------------------------------------------------: c---- Update right hand side c-------------------------------------------------------------: psi(i,j) = f(i,j) f(i,j) = f1(-problem,f(i,j),f(i,j),alpha,beta,dt) end do end do c-----------------------------------------------------------------: c End Tau Loop, only for epsilon c-----------------------------------------------------------------: end do c-----------------------------------------------------------------: c Write mean_psi, F and the solution to file for every step c number of timesteps c-----------------------------------------------------------------: c c IF (MOD(t,step) .EQ. 0 .and. t .ne. 0) THEN write(*,100)(t*dt) c DO i=1,5*(m+1)+1 c DO j=1,5*(n+1)+1 c psi(i,j) = approxit(f,mmax,nmax,m,n,a,b,c,d,a+(i-1)*dx,c+(j c -1)*dy) cc IF (psi(i,j) .GT. -1d-99 .AND. psi(i,j) .LT. 1d-99) THEN cc psi(i,j) = 0 cc IF (MOD(l,10) .EQ. 0) THEN cc WRITE(10,*)psi(i,j) cc END IF cc END IF c if(t .eq. 0) then c write(30,*)psi(i,j) c endif c ENDDO c ENDDO mean_psi = 0 f_tau = 0 DO i=1,(m+1)+1 DO j=1,(n+1)+1 c-----------------------------------------------------------------: c Mean value for psi^2 c-----------------------------------------------------------------: IF(j .LT. (n+1)+1 .AND. i .LT. (m+1)+1) THEN p = (psi(i,j)+psi(i+1,j)+psi(i+1,j+1)+psi(i,j+1) $ )/4 mean_psi = mean_psi + p*p*dx*dy c-----------------------------------------------------------------: c Energy-function c The Laplace operator is approximated using a second order c finite difference scheme based on Taylor series. c-----------------------------------------------------------------: c IF(j .NE. 1 .AND. i .NE. 1) THEN c del_psi=(psi(i+1,j)-2*psi(i,j)+psi(i-1,j)) c $ /(dx**2)+(psi(i,j+1)-2*psi(i,j)+psi(i,j c $ -1))/(dy**2) c c f_tau= f_tau+0.5*(0.5*psi(i,j)**4-epsilon c $ *psi(i,j)**2+(del_psi+psi(i,j))**2)*dx c $ *dy c END IF END IF ENDDO ENDDO mean_psi = mean_psi/((b-a)*(d-c)) c WRITE(20,*)t*dt,mean_psi,f_tau c END IF c-------------------------------------------------------------: c---- End loop c-------------------------------------------------------------: WRITE(*,*)epsilon,mean_psi WRITE(15,*)epsilon,mean_psi END DO c-------------------------------------------------------------: c---- Info output of bihar c-------------------------------------------------------------: close(10) close(15) close(20) c close(30) c close(40) c-------------------------------------------------------------: c---- Done c-------------------------------------------------------------: write(*,101)count write(*,*)'Done.' stop end