!*******************************************
!* !* Solve the advection PDE using Explicit FTCS, !* Explicit Lax, Implicit FTCS, and implicit Crank-Nicolson !* methods for constant and varying speed. !* !* Solve dc/dt = -u dc/dx !* u = t/20 ft/minute !* and !* u constant !* !* Compiler used: gnu 95 (g95) on Cygwin. Gcc 3.4.4 !* Date: June 20 2006 !* !* by Nasser Abbasi !******************************************* PROGRAM advection IMPLICIT NONE REAL :: DT,DX,max_run_time,length,snapshot_delta, & first_limit,second_limit INTEGER :: N,SNAPSHOTS character(10) :: cmd_arg ! to read time step from command line INTEGER :: method ! 1=FTCS, 2=LAX, 3=Implicit FTCS, 4=C-R INTEGER :: mode ! 1=Fixed wind speed, 2=speed function of time REAL :: t_start, t_end, cpu_time_used,end_line(1002) INTEGER :: ALL_DATA_FILE_ID PARAMETER(ALL_DATA_FILE_ID=900) ! Initialize data. All methods will use the same ! parameters to make comparing them easier ! read delta t from command line. CALL getarg(1,cmd_arg) cmd_arg=TRIM(cmd_arg) print *,'= ', cmd_arg read(cmd_arg,*)dt !delta in time, in minutes print *,'Dt=',DT N = 1000 ! number of grid points in space length = 100 ! length of space solution in feet first_limit = 0.25*length second_limit = 0.35*length DX = length/N ! delta in space, in feets max_run_time = 30.0 ! how long to run for in minutes SNAPSHOTS = 200 ! number of snapshots per run. Used for animation snapshot_delta = max_run_time / SNAPSHOTS ! time between each snap shot print *,'DT=',DT,' minutes, DX=',DX,' feets' print *,'taking snapshots every ', snapshot_delta ,' minutes' DO mode=1,2 print*,'=======> processing mode ',mode DO method=1,4 ! No enumeration data types in Fotran 90 CALL CPU_TIME(t_start) ! get current CPU time CALL process(mode,method,N,DT,DX,max_run_time,snapshot_delta,& first_limit,second_limit) CALL CPU_TIME(t_end) ! get current CPU time cpu_time_used = t_end - t_start WRITE(*,FMT='(A,I2,A,F12.5)') 'CPU TIME used for method', method, ' = ', cpu_time_used ! Now record test case parameters in last line end_line=0 end_line(1)=cpu_time_used end_line(2)=DT end_line(3)=DX end_line(4)=mode end_line(5)=method WRITE(UNIT=ALL_DATA_FILE_ID,FMT=*) end_line CLOSE(ALL_DATA_FILE_ID) END DO END DO END PROGRAM advection !************************************ !* !* !************************************ SUBROUTINE process(mode,method,N,DT,DX,max_run_time,snapshot_delta,& first_limit,second_limit) IMPLICIT NONE INTEGER, INTENT(IN) :: mode,method,N REAL, INTENT(IN) :: DT,DX,max_run_time,snapshot_delta,& first_limit,second_limit INTEGER :: I LOGICAL :: snap_shot_at_15_taken INTEGER :: ALL_DATA_FILE_ID PARAMETER(ALL_DATA_FILE_ID=900) REAL :: snap_current_time REAL :: current_time REAL :: C(N) ! current solution REAL :: CNEW(N) ! future solution REAL :: CEXACT(N) ! current exact solution REAL :: current_first_limit REAL :: A(N,N),aa(N),b(2:N),cc(N-1),CTEMP(N) ! for C-R and implicit FTCS REAL :: K,speed REAL :: error,RMS ! root mean square error between current and initial sol. current_time = 0. snap_current_time = 0. CALL initialize_solution(C,N,DX,first_limit,second_limit) CEXACT = C current_first_limit = first_limit CALL pre_loop_initialization(mode,method,current_time,K, & DT,DX,N,C,ALL_DATA_FILE_ID, & A,aa,b,cc ) snap_shot_at_15_taken=.FALSE. DO WHILE(current_time < max_run_time) IF( snap_current_time >= snapshot_delta ) THEN snap_current_time = 0. WRITE( UNIT=ALL_DATA_FILE_ID, FMT=*) current_time, error, C END IF SELECT CASE(method) CASE( 1:2 ) IF(method==1) THEN ! ftcs IF(mode==2)THEN K = speed(mode,current_time)*DT/(2.*DX) ENDIF DO I = 2,N-1 CNEW(I) = C(I) - K * ( C(I+1) - C(I-1) ) END DO ELSE !lax IF(mode == 2) THEN K = speed(mode,current_time)*DT/(DX) ENDIF DO I = 2,N-1 CNEW(I) = C(I) - K/2. * ( C(I+1) - C(I-1) ) + & (K**2.)/2 * ( C(I+1) +C(I-1)-2.*C(I) ) END DO END IF CNEW(1) = C(1) CNEW(N) = C(N) ! Boundary conditions C=CNEW CASE( 3 ) ! implicit ftcs IF( mode == 2) THEN ! only need to update Matrix for varying U K = speed(mode,current_time)*DT/(2.*DX) CALL init_A_matrix(A,K,N) CALL init_diagonal_vectors(N,A,cc,aa,b) END IF CALL solve_thomas_algorithm(N,aa,b,cc,C,CNEW) C = CNEW CASE( 4 ) ! C-R IF(mode == 2) THEN !only need to update A if U changes K = speed(mode,current_time)*DT/(4*DX) ! C-R CALL init_A_matrix(A,K,N) CALL init_diagonal_vectors(N,A,cc,aa,b) END IF CTEMP(1) = C(1) CTEMP(N) = C(N) DO I=2,N-1 CTEMP(I)=C(I)+K*C(I-1)-K*C(I+1) END DO CALL solve_thomas_algorithm(N,aa,b,cc,CTEMP,C) END SELECT IF( current_time>=15.0 .AND. (.NOT. snap_shot_at_15_taken)) THEN snap_shot_at_15_taken = .TRUE. CALL take_one_snap_shot(mode,method,15,N,C,DX) END IF current_time = current_time + DT current_first_limit = current_first_limit + speed(mode,current_time)*DT CALL get_current_exact_solution(CEXACT,N,current_first_limit,DX) error = RMS(CEXACT,C,N) snap_current_time = snap_current_time + DT END DO CALL take_one_snap_shot(mode,method,30,N,C,DX) END SUBROUTINE process !************************************ !* !* !************************************ SUBROUTINE pre_loop_initialization(mode,method,current_time,K,& DT,DX,N,C,ALL_DATA_FILE_ID,& A,aa,b,cc) IMPLICIT NONE INTEGER, INTENT(IN) :: mode,method,N,ALL_DATA_FILE_ID REAL, INTENT(IN) :: C(N),DT,DX,current_time REAL, INTENT(OUT) :: K,A(N,N),aa(N),b(2:N),cc(N-1) REAL :: speed SELECT CASE(method) CASE( 1 ) ! FTCS K = speed(mode,current_time)*DT/(2.*DX) IF(mode==1) THEN OPEN(UNIT=ALL_DATA_FILE_ID, file='expAll.txt') ! all time shots CALL print_to_file(C,'exp0.txt',N,DX) ELSE OPEN(UNIT=ALL_DATA_FILE_ID, file='exp_extraAll.txt') ! all time shots CALL print_to_file(C,'exp_extra0.txt',N,DX) END IF CASE( 2 ) ! Lax K = speed(mode,current_time)*DT/(DX) IF(mode==1) THEN OPEN(UNIT=ALL_DATA_FILE_ID, file='laxAll.txt') ! all time shots CALL print_to_file(C,'lax0.txt',N,DX) ELSE OPEN(UNIT=ALL_DATA_FILE_ID, file='lax_extraAll.txt') ! all time shots CALL print_to_file(C,'lax_extra0.txt',N,DX) END IF CASE( 3 ) ! Implicit FTCS K = speed(mode,current_time)*DT/(2.*DX) CALL init_A_matrix(A,K,N) CALL init_diagonal_vectors(N,A,cc,aa,b) IF(mode==1) THEN OPEN(UNIT=ALL_DATA_FILE_ID, file='impAll.txt') ! all time shots CALL print_to_file(C,'imp0.txt',N,DX) ELSE OPEN(UNIT=ALL_DATA_FILE_ID, file='imp_extraAll.txt') ! all time shots CALL print_to_file(C,'imp_extra0.txt',N,DX) END IF CASE( 4 ) ! C-R K = speed(mode,current_time)*DT/(4*DX) ! C-R CALL init_A_matrix(A,K,N) CALL init_diagonal_vectors(N,A,cc,aa,b) IF(mode==1) THEN OPEN(UNIT=ALL_DATA_FILE_ID, file='crAll.txt') ! all time shots CALL print_to_file(C,'cr0.txt',N,DX) ELSE OPEN(UNIT=ALL_DATA_FILE_ID, file='cr_extraAll.txt') ! all time shots CALL print_to_file(C,'cr_extra0.txt',N,DX) END IF END SELECT WRITE( UNIT=ALL_DATA_FILE_ID, FMT=*) current_time,0, C END SUBROUTINE pre_loop_initialization !************************************ !* !* !************************************ SUBROUTINE init_diagonal_vectors(N,A,cc,aa,b) IMPLICIT NONE INTEGER, INTENT(IN) ::N REAL, INTENT(IN) ::A(N,N) REAL, INTENT(OUT) ::aa(N),b(2:N),cc(N-1) INTEGER ::I,J J=2 DO I=1,N-1 cc(I)=A(I,J) J=J+1 END DO cc(1)=0 DO I=1,N aa(I)=A(I,I) END DO J=1 DO I=2,N b(I)=A(I,J) J=J+1 END DO END SUBROUTINE init_diagonal_vectors !************************************ !* !* !************************************ SUBROUTINE initialize_solution(C,N,DX,first_limit,second_limit) IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: DX,first_limit,second_limit REAL, INTENT(INOUT) :: C(0:N-1) INTEGER :: I REAL :: x, PI,av,R PARAMETER( PI = ACOS(-1.) ) x = 0 av = (second_limit+first_limit)/2.0 R = av - first_limit C = 0.0 DO I=0,N-1 IF( x >= first_limit .AND. x <= second_limit ) THEN C(I) = 1 + COS( PI * (x-av)/R ) END IF x = x + DX END DO END SUBROUTINE initialize_solution !************************************ !* !* !************************************ SUBROUTINE print_to_file(C,file_name,N,DX) IMPLICIT NONE REAL, INTENT(IN) :: C(N),DX INTEGER, INTENT(IN) :: N CHARACTER* (*), INTENT(IN) :: file_name INTEGER :: I INTEGER :: FILE_ID PARAMETER(FILE_ID=999) REAL :: current_position OPEN(UNIT=FILE_ID, file=file_name) current_position = 0; DO I=1,N WRITE( UNIT=FILE_ID, FMT=* ) current_position ,'\t', C(I) current_position = current_position + DX END DO CLOSE(FILE_ID) END SUBROUTINE print_to_file !************************************ !* !* !************************************ SUBROUTINE init_A_matrix(A,K,N) IMPLICIT NONE INTEGER, INTENT(IN) ::N REAL, INTENT(IN) ::K REAL, INTENT(OUT) ::A(N,N) INTEGER ::I DO I = 2,N-1 A(I,I-1) = -K A(I,I) = 1 A(I,I+1) = K END DO A(1,1) = 1 A(N,N) = 1 END SUBROUTINE init_A_matrix !************************************ !* !* !************************************ SUBROUTINE solve_thomas_algorithm(N,aa,b,c,old_c,new_c) IMPLICIT NONE REAL, INTENT(IN) :: aa(N),b(2:N),c(N-1),old_c(N) INTEGER, INTENT(IN) :: N REAL, INTENT(INOUT) :: new_c(N) INTEGER :: I REAL :: alpha(N),beta(2:N),g(N) alpha(1) = aa(1) DO I=2,N beta(I)=b(I)/alpha(I-1) alpha(I)=aa(I)-beta(I)*c(I-1) END DO g(1)=old_c(1) DO I=2,N g(I)=old_c(I)-beta(I)*g(I-1) END DO new_c(N)=g(N)/alpha(N) DO I=N-1,1,-1 new_c(I)=(g(I)-c(I)*new_c(I+1))/alpha(I) END DO END SUBROUTINE solve_thomas_algorithm !************************************ !* !* !************************************ REAL FUNCTION speed(MODE,time) IMPLICIT NONE INTEGER, INTENT(IN) :: MODE REAL, INTENT(IN) :: time IF( MODE == 1 ) THEN speed=2.0 ELSE speed=time/20.0 END IF END FUNCTION speed !************************************ !* !* !************************************ SUBROUTINE take_one_snap_shot(mode,method,TIME,N,C,DX) IMPLICIT NONE INTEGER, INTENT(IN) ::TIME,mode,method,N REAL, INTENT(IN) ::C(N),DX IF(TIME==15) THEN SELECT CASE(method) CASE(1) IF(mode==1) THEN CALL print_to_file(C,'exp15.txt',N,DX) ELSE CALL print_to_file(C,'exp_extra15.txt',N,DX) END IF CASE(2) IF(mode==1) THEN CALL print_to_file(C,'lax15.txt',N,DX) ELSE CALL print_to_file(C,'lax_extra15.txt',N,DX) ENDIF CASE(3) IF(mode==1) THEN CALL print_to_file(C,'imp15.txt',N,DX) ELSE CALL print_to_file(C,'imp_extra15.txt',N,DX) END IF CASE(4) IF(mode==1) THEN CALL print_to_file(C,'cr15.txt',N,DX) ELSE CALL print_to_file(C,'cr_extra15.txt',N,DX) END IF END SELECT ELSE SELECT CASE(method) CASE(1) IF(mode==1) THEN CALL print_to_file(C,'exp30.txt',N,DX) ELSE CALL print_to_file(C,'exp_extra30.txt',N,DX) END IF CASE(2) IF(mode==1) THEN CALL print_to_file(C,'lax30.txt',N,DX) ELSE CALL print_to_file(C,'lax_extra30.txt',N,DX) ENDIF CASE(3) IF(mode==1) THEN CALL print_to_file(C,'imp30.txt',N,DX) ELSE CALL print_to_file(C,'imp_extra30.txt',N,DX) END IF CASE(4) IF(mode==1) THEN CALL print_to_file(C,'cr30.txt',N,DX) ELSE CALL print_to_file(C,'cr_extra30.txt',N,DX) END IF END SELECT END IF END SUBROUTINE take_one_snap_shot !************************************ !* !* !************************************ REAL FUNCTION RMS(CEXACT,C,N) IMPLICIT NONE REAL, INTENT(IN) :: CEXACT(N),C(N) INTEGER, INTENT(IN) :: N INTEGER :: I RMS=0. DO I=1,N RMS = RMS+(CEXACT(I)-C(I))**2 END DO RMS = RMS/N RMS = SQRT(RMS) END FUNCTION RMS !************************************ !* !* !************************************ SUBROUTINE get_current_exact_solution(CEXACT,N,current_first_limit,DX) IMPLICIT NONE REAL, INTENT(IN) :: current_first_limit,DX REAL, INTENT(OUT) :: CEXACT(0:N-1) INTEGER, INTENT(IN) :: N INTEGER :: I REAL :: first_limit REAL :: second_limit REAL :: av,R,shift,x,PI PARAMETER( PI = ACOS(-1.) ) first_limit = 25.0 second_limit = 35.0 shift = current_first_limit - FIRST_LIMIT first_limit = current_first_limit second_limit = second_limit + shift av = (second_limit+first_limit)/2.0 R = av - first_limit CEXACT = 0. x = 0. DO I = 0,N-1 IF( x >= first_limit .AND. x <= second_limit ) THEN CEXACT(I) = 1 + COS( PI * (x -av)/R ) END IF x = x + DX END DO END SUBROUTINE get_current_exact_solution