C*********************************************************************** SUBROUTINE INITIAL(DT,DX,Nx,Q,Delta,beta,Cloak) IMPLICIT NONE INTEGER j, Nx, Nnon0, Nfront, Nback REAL*8 DT, DX, Q(-1:1,-Nx:Nx), Delta, Qdum, front, back, test, / Sdt, beta, Cloak(-Nx:Nx) C Set the initial conditions of the Solution DO j = -Nx, Nx Q(1,j) = 0.0D0 Cloak(j) = Q(1,j) ENDDO Nnon0 = INT(Nx*0.2D0) DO j = -Nnon0, Nnon0 Q(1,j) = 0.04D0*EXP( -(DX*j)**2 / Delta**2 ) Cloak(j) = Q(1,j) ENDDO C Now determine your backwards step using the analytical solution. 400 DO j = -Nx, Nx Q(0,j) = 0.0D0 ENDDO Sdt = ( 1.0D0 - beta )*DT back = -0.2D0 - Sdt Nback = INT(back/DX) test = ABS( NINT(back/DX) - back/DX ) IF ( test .LE. 1.0D-7 ) Nback = NINT(back/DX) front = 0.2D0 - Sdt Nfront = INT(front/DX) test = ABS( NINT(front/DX) - front/DX ) IF ( test .LE. 1.0D-7 ) Nfront = NINT(front/DX) DO j = Nback, Nfront Q(0,j) = 0.04D0*EXP( -(DX*j + Sdt)**2 / Delta**2 ) ENDDO RETURN END C*********************************************************************** SUBROUTINE MATRIX(DT,DX,thetax,thetat,beta,Nbound,Nuse,Mat,Work, / Ord) IMPLICIT NONE INTEGER j, k, Nuse, Nbound, info REAL*8 DT, DX, rho, thetax, thetat, beta, Mat(Nuse,Nuse), / Work(Nuse), rcond, thxa INTEGER Ord(Nuse) C Set the initial conditions of the matrix. DO j = 1, Nuse Work(j) = 0.0D0 DO k = 1, Nuse Mat(j,k) = 0.0D0 ENDDO ENDDO rho = DT / DX thxa = thetax*( 1.0D0 - beta**2 )*rho**2 IF ( Nbound .EQ. -1 ) THEN DO j = 1, Nuse IF ( j .NE. 1 ) / Mat(j,j-1) = 0.5D0*( thetat - thxa + beta*rho ) Mat(j,j) = 1.0D0 - thetat + thxa IF ( j .NE. Nuse ) / Mat(j,j+1) = 0.5D0*( thetat - thxa - beta*rho ) ENDDO ELSEIF ( Nbound .EQ. 0 ) THEN Mat(1,Nuse) = 0.5D0*( thetat - thxa + beta*rho ) Mat(1,1) = 1.0D0 - thetat + thxa Mat(1,2) = 0.5D0*( thetat - thxa - beta*rho ) DO j = 2, Nuse Mat(j,j-1) = 0.5D0*( thetat - thxa + beta*rho ) Mat(j,j) = 1.0D0 - thetat + thxa IF ( j .EQ. Nuse ) THEN Mat(j,1) = 0.5D0*( thetat - thxa - beta*rho ) ELSE Mat(j,j+1) = 0.5D0*( thetat - thxa - beta*rho ) ENDIF ENDDO ELSEIF ( Nbound .EQ. 1 ) THEN IF ( beta .LE. 1.0D0 ) THEN Mat(1,1) = 3.0D0*( rho*( 1.0D0 + beta ) + 1.0D0 ) Mat(1,2) = -4.0D0*rho*( 1.0D0 + beta ) Mat(1,3) = rho*( 1.0D0 + beta ) ELSE Mat(1,Nuse-2) = rho*( 1.0D0 + beta) Mat(1,Nuse-1) = -4.0D0*rho*( 1.0D0 + beta ) Mat(1,Nuse) = 3.0D0*( rho*( 1.0D0 + beta ) - 1.0D0 ) ENDIF DO j = 2, Nuse-1 Mat(j,j-1) = 0.5D0*( thetat - thxa + beta*rho ) Mat(j,j) = 1.0D0 - thetat + thxa Mat(j,j+1) = 0.5D0*( thetat - thxa - beta*rho ) ENDDO IF ( beta .GE. -1.0D0 ) THEN Mat(Nuse,Nuse-2) = rho*( 1.0D0 - beta ) Mat(Nuse,Nuse-1) = -4.0D0*rho*( 1.0D0 - beta ) Mat(Nuse,Nuse) = 3.0D0*( rho*( 1.0D0 - beta ) + 1.0D0 ) ELSE Mat(Nuse,1) = 3.0D0*( rho*( 1.0D0 - beta ) - 1.0D0 ) Mat(Nuse,2) = -4.0D0*rho*( 1.0D0 - beta ) Mat(Nuse,3) = rho*( 1.0D0 - beta ) ENDIF ELSE PRINT *, 'Invalid value of Nbound' STOP ENDIF C CALL DGEFA(Mat,Nuse,Nuse,Ord,INFO) C IF ( INFO .NE. 0 ) THEN C PRINT *, 'The Inversion has failed.' C STOP C ENDIF CALL DGECO(Mat,Nuse,Nuse,Ord,rcond,work) PRINT 15, 1.0D0 / rcond 15 FORMAT(' The condition number = ',G15.7) RETURN END C*********************************************************************** SUBROUTINE EVOL(DT,DX,thetax,thetat,beta,Nbound,Nx,Q,Nuse,Mat, / Vec,Ord,Cloak,Ncloak) IMPLICIT NONE INTEGER j, k, Nx, Nuse, Nbound, Ncloak REAL*8 DT, DX, thetax, thetat, beta, rho, Q(-1:1,-Nx:Nx), / Mat(Nuse,Nuse), Vec(Nuse), Cloak(-Nx:Nx), thxa, thxb, / xpt, Qp, Qs, Qr INTEGER Ord(Nuse) C Evolve to the next timestep. DO j = -Nx, Nx Q(-1,j) = Q(0,j) Q(0,j) = Q(1,j) ENDDO rho = ( DT / DX ) thxa = thetax*( 1.0D0 - beta**2 )*rho**2 thxb = ( 1.0D0 - thetax )*( 1.0D0 - beta**2 )*rho**2 IF ( Nbound .EQ. -1 ) THEN DO j = -Nx+1, Nx-1 k = j + Nx Vec(k) = ( thxb + thetat )*( Q(0,j+1) + Q(0,j-1) ) / + 2.0D0*( 1.0D0 - thetat - thxb )*Q(0,j) / - ( thxa + 1.0D0 - thetat )*Q(-1,j) / + 0.5D0*( thxa - thetat )*( Q(-1,j+1) + Q(-1,j-1) ) / - 0.5D0*beta*rho*( Q(-1,j+1) - Q(-1,j-1) ) ENDDO ELSEIF ( Nbound .EQ. 0 ) THEN Vec(1) = ( thxb + thetat )*( Q(0,-Nx+1) + Q(0,Nx-1) ) / + 2.0D0*( 1.0D0 - thetat - thxb )*Q(0,-Nx) / - ( thxa + 1.0D0 - thetat )*Q(-1,-Nx) / + 0.5D0*( thxa - thetat )*( Q(-1,-Nx+1) + Q(-1,Nx-1) ) / - 0.5D0*beta*rho*( Q(-1,-Nx+1) - Q(-1,Nx-1) ) DO j = -Nx+1, Nx-1 k = j + Nx + 1 Vec(k) = ( thxb + thetat )*( Q(0,j+1) + Q(0,j-1) ) / + 2.0D0*( 1.0D0 - thetat - thxb )*Q(0,j) / - ( thxa + 1.0D0 - thetat )*Q(-1,j) / + 0.5D0*( thxa - thetat )*( Q(-1,j+1) + Q(-1,j-1) ) / - 0.5D0*beta*rho*( Q(-1,j+1) - Q(-1,j-1) ) ENDDO ELSE IF ( beta .LE. 1.0D0 ) THEN Vec(1) = 4.0D0*Q(0,-Nx) - Q(-1,-Nx) ELSE Vec(1) = Q(-1,Nx) - 4.0D0*Q(0,Nx) ENDIF DO j = -Nx+1, Nx-1 k = j + Nx + 1 Vec(k) = ( thxb + thetat )*( Q(0,j+1) + Q(0,j-1) ) / + 2.0D0*( 1.0D0 - thetat - thxb )*Q(0,j) / - ( thxa + 1.0D0 - thetat )*Q(-1,j) / + 0.5D0*( thxa - thetat )*( Q(-1,j+1) + Q(-1,j-1) ) / - 0.5D0*beta*rho*( Q(-1,j+1) - Q(-1,j-1) ) ENDDO IF ( beta .GE. -1.0D0 ) THEN Vec(Nuse) = 4.0D0*Q(0,Nx) - Q(-1,Nx) ELSE Vec(Nuse) = Q(-1,-Nx) - 4.0D0*Q(0,-Nx) ENDIF ENDIF CALL DGESL(Mat,Nuse,Nuse,Ord,Vec,0) IF ( Nbound .EQ. -1 ) THEN Q(1,-Nx) = 0.0D0 DO j = -Nx+1, Nx-1 k = j + Nx Q(1,j) = Vec(k) ENDDO Q(1,Nx) = 0.0D0 ELSEIF ( Nbound .EQ. 0 ) THEN DO j = -Nx, Nx-1 k = j + Nx + 1 Q(1,j) = Vec(k) ENDDO Q(1,Nx) = Q(1,-Nx) ELSE DO j = -Nx, Nx k = j + Nx + 1 Q(1,j) = Vec(k) ENDDO ENDIF IF ( Ncloak .EQ. 0 ) GOTO 500 C Smooth the solution if necessary. IF ( Nbound .EQ. -1 ) THEN Cloak(-Nx) = 0.0D0 Cloak(-Nx+1) = 0.5D0*Vec(1) + 0.25D0*Vec(2) DO j = -Nx+2, Nx-2 k = j + Nx Cloak(j) = 0.25D0*Vec(k-1) + 0.5D0*Vec(k) + 0.25D0*Vec(k+1) ENDDO Cloak(Nx-1) = 0.25D0*Vec(Nuse-1) + 0.5D0*Vec(Nuse) Cloak(Nx) = 0.0D0 ELSEIF ( Nbound .EQ. 0 ) THEN Cloak(-Nx) = 0.25D0*Vec(Nuse) + 0.5D0*Vec(1) + 0.25D0*Vec(2) DO j = -Nx+1, Nx-2 k = j + Nx + 1 Cloak(j) = 0.25D0*Vec(k-1) + 0.5D0*Vec(k) + 0.25D0*Vec(k+1) ENDDO Cloak(Nx-1) = 0.25D0*Vec(Nuse-1) + 0.5D0*Vec(Nuse) / + 0.25D0*Vec(1) Cloak(Nx) = Cloak(-Nx) ELSE Cloak(-Nx) = Q(1,-Nx) Cloak(-Nx+1) = 0.5D0*Vec(2) + 0.5D0*Vec(3) DO j = -Nx+2, Nx-2 k = j + Nx + 1 Cloak(j) = 0.25D0*Vec(k-1) + 0.5D0*Vec(k) + 0.25D0*Vec(k+1) ENDDO Cloak(Nx-1) = 0.6D0*Vec(Nuse-2) + 0.4D0*Vec(Nuse-1) Cloak(Nx) = Q(1,Nx) ENDIF RETURN 500 DO j = -Nx, Nx Cloak(j) = Q(1,j) ENDDO RETURN END