PROGRAM BISPEC1 ! ! ! 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, & petteru,petterfrhs,error ! 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 INTEGER i,ii,flag,j,l INTEGER i_save,ii_save,nstop DOUBLE PRECISION alpha,beta,k,x,y,erel,el2 ! 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 :: svar DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: hoyreut DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: svarut double precision, dimension(:,:), allocatable :: fverdier double precision, dimension(:,:), allocatable :: uverdier 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 100 FORMAT(I6,' ',E16.8,' ',E16.8,' ',E16.8,' ',E16.8) 101 FORMAT(I4,' ',I4,' ',E16.8,' ',E16.8) 102 FORMAT(E16.8,' ',E16.8,' ',E16.8,' ',E16.8) problem = 1 method = 1 WRITE(*,'(A)',ADVANCE = 'NO') 'Enter problem size: ' n = 98 ! read(5,*)n WRITE(*,*)'n=',n WRITE(*,*) 'The problem is solved on the area [a,b]x[c,d]' write(*,'(A)',ADVANCE = 'NO') ' Enter a b c d: ' ! read(5,*) au,bu,cu,du au = -1 bu = 1 cu = -1 du = 1 WRITE(*,*)'a,b,c,d',au,bu,cu,du WRITE(*,'(A)',ADVANCE = 'NO') 'Enter konstant: ' k = 30 ! read(5,*)k WRITE(*,*)'Constant k = ',k nstop = 20 WRITE(*,*)'nstop = ',nstop x1 = 4.d0/((du-cu)**2) x2 = 4.d0/((bu-au)**2) alpha = pde_alpha beta = pde_beta !-----------------------------------------------------------------: ! Write output file !-----------------------------------------------------------------: WRITE(*,*)'Writing outputfile.' OPEN(unit=10,file='bis_output') ! WRITE(10,*)au ! WRITE(10,*)bu ! WRITE(10,*)cu ! WRITE(10,*)du ! WRITE(10,*)n*2+2 ! WRITE(10,*)n*2+2 ! WRITE(10,*)nstop !-----------------------------------------------------------------: ! Loop !-----------------------------------------------------------------: DO l=101,150 ! k = l*5 ! IF(l .EQ. 0) k = 1 n = l IF(method.EQ.1.OR.method.EQ.4) THEN legtrans = .TRUE. ELSE legtrans = .FALSE. END IF ! 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 (svar(0:2*N+3,0:2*N+3)) ALLOCATE (hoyreut(0:(2*N+3)**2)) ALLOCATE (svarut(0:(2*N+3)**2)) ALLOCATE (fverdier(0:2*N+3,0:2*N+3)) ALLOCATE (uverdier(0:2*N+3,0:2*N+3)) ! ! ! Compute the appropriate right hand side for ! our linear system and store it in hoyre. ! As part of this we also need to compute the relevant coordinates xjac. ! Also get the solution u (if available), for this problem ! and store it in svar. ! ! call setup(N,xjac,hoyre,svar,problem,legtrans) ! ! Beregner pkt. xjac. ! 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 !-----------------------------------------------------------------: ! Danner hoyreside og svar: !-----------------------------------------------------------------: DO j = 0,2*N+3 DO i = 0,2*N+3 hoyre(i,j) = petterfrhs(xjacx(i),xjacy(j),au,bu,cu,du,alpha,beta,k) svar(i,j) = petteru(xjacx(i),xjacy(j),au,bu,cu,du,k) ! WRITE(10,102)xjacx(i),xjacy(j),loes1(i,j),svar(i,j) END DO END DO flag = 0 loes1 = hoyre method = 1 etime = gettime() CALL leg(loes1,coeff,xjac,alpha,beta,flag,N,x1,x2,tid1,tid2,tid3,tid4) etime = gettime()-etime CALL adjust(N,xjac,loes1,problem) !-----------------------------------------------------------------: ! Uniform coords !-----------------------------------------------------------------: DO j=0,2*N+3 y = cu+(du-cu)/(2*N+4)*j DO i=0,2*N+3 x = au+(bu-au)/(2*N+4)*i fverdier(i,j)=funkval(x,y,coeff,2*N+4) uverdier(i,j)=petteru(x,y,au,bu,cu,du,k) ! IF(i .NE. 0 .AND. j .NE. 0 .AND. i .NE. 2*N+3 .AND. j& ! .NE. 2*N+3 .AND. j .NE. n+2 .AND. i .NE.& ! n+2.AND. j .NE. n+1 .AND. i .NE. n+1) then ! WRITE(10,101)i+1,j+1,fverdier(i,j),uverdier(i,j) ! END if ENDDO ENDDO CALL error(uverdier,fverdier,N,au,bu,cu,du,max1,el2,erel) WRITE(*,100)INT(n*2+2),max1,el2,erel,etime ! WRITE(*,100)INT((n*2+2)**2),max1,el2,erel,etime ! WRITE(10,100)INT(k),max1,el2,erel,etime DEALLOCATE(hoyre,xjac,xjacx,xjacy,loes1,coeff,loes2,svar,hoyreut,svarut,fverdier,uverdier) ! IF(erel .LT. 0.0001) stop END DO CLOSE(10) END PROGRAM bispec1