PROGRAM Analyze_Velocity_Evolution !by Peter Bird, UCLA, 2003.09.27 !This utility program is useful for gaining insight into those (rare) cases in ! in which program NeoKinema does not converge to a solution as it iterates. !The user can re-run NeoKinema with the final switch, "dump_all_solutions" ! (in the input parameter file p_[token].nki) set = .TRUE., so that ! velocity solutions from all iterations will be written to v_log.nko. !Then, this utility program can be run to analyze that giant v_log.nko file. IMPLICIT NONE INTEGER :: bytes_used, & & components_wanted, & & first, & & i, i_from_ic, i_from_jc, ic, ios, & & j, j_from_ic, j_from_jc, jc, & & k, last, & & MB_used, & & ndof, num_nod, number, & & p, passes, pcap, & & v REAL*8 :: ic_noise, jc_noise, lever, middle, range, slope_b, sum, sum_x, sum_xy, sum_x2, sum_y REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: v_mps ! input: all velocity solutions REAL*8, DIMENSION(:,:), ALLOCATABLE :: v_mean, v_trend, v_mode ! outputs : computed velocity-like solution (mean, trend, mode, ...) REAL*8, DIMENSION(:,:), ALLOCATABLE :: covariance, eigenvectors, recipes REAL*8, DIMENSION(:), ALLOCATABLE :: eigenvalues WRITE (*, "(' PROGRAM Analyze_Velocity_Evolution')") WRITE (*, "(' by Peter Bird, UCLA, 2003.09.27')") WRITE (*, "(' ')") WRITE (*, "(' This utility program is useful for gaining insight into those (rare) cases in')") WRITE (*, "(' in which program NeoKinema does not converge to a solution as it iterates.')") WRITE (*, "(' The user can re-run NeoKinema with the final switch, ""dump_all_solutions""')") WRITE (*, "(' (in the input parameter file p_[token].nki) set = .TRUE., so that')") WRITE (*, "(' velocity solutions from all iterations will be written to v_log.nko.')") WRITE (*, "(' Then, this utility program can be run to analyze that giant v_log.nko file.')") WRITE (*, "(' ')") CALL Pause() OPEN (UNIT = 1, FILE = "v_log.nko", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' Error: v_log.nko not found in current directory.')") CALL Pause() STOP END IF READ (1, *, IOSTAT = ios) num_nod, passes IF (ios /= 0) THEN WRITE (*, "(' Error: Empty or malformed input file: v_log.nko')") CALL Pause() STOP END IF ALLOCATE ( v_mps(2, num_nod, passes) ) DO i = 1, passes READ (1, *, IOSTAT = ios) ! (pass number = i) IF (ios /= 0) THEN WRITE (*, "(' Error: Incomplete v_log.nko file; adjust #passes in first line')") CALL Pause() STOP END IF DO j = 1, num_nod READ (1, *, IOSTAT = ios) k, v_mps(1, j, i), v_mps(2, j, i) ! k should be == j IF (ios /= 0) THEN WRITE (*, "(' Error: Incomplete v_log.nko file; adjust #passes in first line')") CALL Pause() STOP END IF END DO END DO CLOSE (1) WRITE (*, "(' v_log.nko has been read and is held in memory.')") WRITE (*, "(' It contains ',I6,' velocity solutions.')") passes 10 WRITE (*, "(' ')") WRITE (*, "(' Enter two integers to identify the group of solutions to be analyzed:')") WRITE (*, "(' Enter the lower integer now: '\)") READ (*, *) first IF ((first < 1).OR.(first > passes)) THEN WRITE (*, "(' Error: Enter an integer from 1 to ',I6,'.')") passes CALL Pause() GO TO 10 END IF 20 WRITE (*, "(' Enter the higher integer now: '\)") READ (*, *) last IF ((last < first).OR.(last > passes)) THEN WRITE (*, "(' Error: Enter an integer from ',I6,' to ',I6,'.')") first, passes CALL Pause() GO TO 20 END IF number = 1 + last - first !Compute the MEAN velocity field in this range (first, last). ALLOCATE ( v_mean(2, num_nod) ) v_mean = 0.0D0 ! whole array, initializing sum DO i = 1, 2 DO j = 1, num_nod DO k = first, last v_mean(i, j) = v_mean(i, j) + v_mps(i, j, k) END DO END DO END DO IF (number > 1) THEN DO i = 1, 2 DO j = 1, num_nod v_mean(i, j) = v_mean(i, j) / number END DO END DO END IF OPEN (UNIT = 2, FILE = "v_mean.out") WRITE (2, "('v_mean.out')") WRITE (2, "('for visualization with FiniteMap, to study convergence of NeoKinema')") WRITE (2, "('computed with Analyze_Velocity_Evolution.f90')") DO j = 1, num_nod WRITE (2, "(ES13.6, 1X, ES13.6)") v_mean(1, j), v_mean(2, j) END DO CLOSE (2) WRITE (*, "(' ')") WRITE (*, "(' See v_mean.out for mean velocity; plot with FiniteMap.')") !Compute the TREND (net change from first to last, according to: ! --> simple difference (if number == 2 OR 3); or ! --> trend lines fit by least-squares (if number >= 4) IF (number > 1) THEN ALLOCATE ( v_trend(2, num_nod) ) IF (number < 4) THEN ! use simple difference between first and last values DO i = 1, 2 DO j = 1, num_nod v_trend(i, j) = v_mps(i, j, last) - v_mps(i, j, first) END DO END DO ELSE ! number >= 4; fit least-squares lines range = last - first ! in "x" variable, which is 1 less than "number" middle = (first + last) / 2.0D0 ! mean value of "pass" variable ( = x ) DO i = 1, 2 ! components of velocity DO j = 1, num_nod sum_x = 0.0D0 sum_x2 = 0.0D0 sum_y = 0.0D0 sum_xy = 0.0D0 DO k = first, last sum_x = sum_x + k sum_x2 = sum_x2 + k * k sum_y = sum_y + v_mps(i, j, k) sum_xy = sum_xy + k * v_mps(i, j, k) END DO slope_b = ( number * sum_xy - sum_x * sum_y ) / & & ( number * sum_x2 - sum_x**2 ) v_trend(i, j) = slope_b * range END DO ! j = 1, num_nod END DO ! i = 1, 2; components of velocity END IF ! number < 4 or >= 4 (two methods) OPEN (UNIT = 3, FILE = "v_trend.out") WRITE (3, "('v_trend.out')") WRITE (3, "('for visualization with FiniteMap, to study convergence of NeoKinema')") WRITE (3, "('computed with Analyze_Velocity_Evolution.f90')") DO j = 1, num_nod WRITE (3, "(ES13.6, 1X, ES13.6)") v_trend(1, j), v_trend(2, j) END DO CLOSE (3) WRITE (*, "(' ')") WRITE (*, "(' See v_trend.out for net velocity shift; plot with FiniteMap.')") WRITE (*, "(' Note: Do NOT change reference frame when plotting;')") WRITE (*, "(' this velocity field is already differential and frame-independent,')") WRITE (*, "(' so changing its frame would yield confusing and misleading plots.')") END IF ! number > 1, so that computing a trend is possible IF (number >= 3) THEN ! eigen-vector analysis is possible: WRITE (*, "(' ')") WRITE (*, "(' Begining the task of computing the largest eigenvector of the solution ""noise""')") WRITE (*, "(' (what is left after subtracting the mean and the trend):')") WRITE (*, "(' in some cases this could be the strongest mode of oscillation;')") WRITE (*, "(' but, in the case of exponential behavior, it might just be the nonlinearity.')") WRITE (*, "(' This task is slow and memory-intensive, and might just crash...')") ndof = 2 * num_nod bytes_used = 8 * (2 * ndof**2 + 2 * ndof) MB_used = bytes_used / (1024 * 1024) WRITE (*, "(' Allocating ',I6,' MB of memory...')") MB_used ALLOCATE ( covariance(ndof, ndof) ) ALLOCATE ( eigenvalues(ndof) ) ALLOCATE ( eigenvectors(ndof, ndof) ) ALLOCATE ( recipes(1, ndof) ) ALLOCATE ( v_mode(2, num_nod) ) !Compute the covariance matrix: WRITE (*, "(' ')") DO ic = 1, ndof WRITE (*, "('+computing covariance matrix row ',I6,' out of ',I6)") ic, ndof IF (MOD(ic, 2) == 1) THEN ! ic is odd: i_from_ic = 1 j_from_ic = 1 + (ic - 1) / 2 ELSE ! ic is even i_from_ic = 2 j_from_ic = ic / 2 END IF DO jc = 1, ndof IF (MOD(jc, 2) == 1) THEN ! jc is odd: i_from_jc = 1 j_from_jc = 1 + (jc - 1) / 2 ELSE ! jc is even i_from_jc = 2 j_from_jc = jc / 2 END IF covariance(ic, jc) = 0.0 ! just initializing before sum DO k = first, last lever = k - middle ic_noise = v_mps(i_from_ic, j_from_ic, k) - & & v_mean(i_from_ic, j_from_ic) - & & lever * v_trend(i_from_ic, j_from_ic) / range jc_noise = v_mps(i_from_jc, j_from_jc, k) - & & v_mean(i_from_jc, j_from_jc) - & & lever * v_trend(i_from_jc, j_from_jc) / range covariance(ic, jc) = covariance(ic, jc) + ic_noise * jc_noise END DO covariance(ic, jc) = covariance(ic, jc) / number END DO END DO !Compute largest eigenvalue and corresponding eigenvector by iterative means: components_wanted = 1 WRITE (*, "(' calling Eigen to analyze for largest eigenvalue...')") CALL DEigen (epsiln = 1.0D-5, length = ndof, nda = ndof, ndv = ndof, & ! intent(IN) & a = covariance, number = components_wanted, & ! intent(INOUT) & values = eigenvalues, vectrs = eigenvectors) ! intent(OUT) WRITE (*, "(/' Eigen found ',I3,' eigenvalues and eigenvectors.')") number WRITE (*, "(' Computing and outputing the result...')") pcap = 1 ! compute one of one eigenvalues ! Compute factor Recipes R DO p = 1, pcap DO v = 1, ndof recipes(p, v) = eigenvectors(v, p) * DSQRT(eigenvalues(p)) END DO END DO DO i = 1, num_nod v_mode(1, i) = recipes(1, 2 * i - 1) v_mode(2, i) = recipes(1, 2 * i) END DO OPEN (UNIT = 4, FILE = "v_mode.out") WRITE (4, "('v_mode.out')") WRITE (4, "('for visualization with FiniteMap, to study convergence of NeoKinema')") WRITE (4, "('computed with Analyze_Velocity_Evolution.f90')") DO j = 1, num_nod WRITE (4, "(ES13.6, 1X, ES13.6)") v_mode(1, j), v_mode(2, j) END DO CLOSE (4) WRITE (*, "(' ')") WRITE (*, "(' See v_mode.out for eigenvector of prime oscillation mode; plot with FiniteMap.')") WRITE (*, "(' Note: Do NOT change reference frame when plotting;')") WRITE (*, "(' this velocity field is already differential and frame-independent,')") WRITE (*, "(' so changing its frame would yield confusing and misleading plots.')") WRITE (*, "(' ')") WRITE (*, "(' Computing amplitude of this eigenvector in each iteration:')") DO k = first, last lever = k - middle sum = 0.0D0 DO ic = 1, ndof IF (MOD(ic, 2) == 1) THEN ! ic is odd: i_from_ic = 1 j_from_ic = 1 + (ic - 1) / 2 ELSE ! ic is even i_from_ic = 2 j_from_ic = ic / 2 END IF ic_noise = v_mps(i_from_ic, j_from_ic, k) - & & v_mean(i_from_ic, j_from_ic) - & & lever * v_trend(i_from_ic, j_from_ic) / range sum = sum + recipes(1, ic) * ic_noise END DO sum = sum / eigenvalues(1) WRITE (*, "(' ',I6,F10.4)") k, sum END DO END IF ! number >= 3, so eigenvector analysis is possible WRITE (*, "(' ')") WRITE (*, "(' Job done; plot resulting file(s) with FiniteMap or NeoKineMap.')") CALL Pause() CONTAINS !_____________________________________________________________ SUBROUTINE DEigen (epsiln, length, nda, ndv, & ! intent(IN) & a, number, & ! intent(IN OUT) & values, vectrs) ! intent(OUT) ! COMPUTE EIGENVALUES AND EIGENVECTORS OF A MATRIX WHICH IS: ! -REAL (NOT COMPLEX); ! -SYMMETRIC (A(I,J)=A(J,I)); AND ! -POSITIVE SEMIDEFINATE (NO NEGATIVE EIGENVALUES, SO THAT ! X (AX) >= 0.0 FOR ANY VECTOR X) ! BY THE METHOD OF POWERS AND SWEEPING-OUT OF KNOWN EIGENVALUES. ! BASED ON M.L. JAMES, G.M. SMITH, & J.C. WOLFORD, ! APPLIED NUMERICAL METHODS FOR DIGITAL COMPUTATION, ! HARPER & ROW, PAGES 214-239. ! DETERMINES EIGENVALUES AND EIGENVECTORS SIMULTANEOUSLY, ! FROM THE LARGEST EIGENVALUE DOWNWARD TOWARD ZERO. ! THIS ROUTINE MAY BE USED TO FIND ALL VALUES AND VECTORS, BUT ! ACCURACY DECLINES AFTER THE FIRST 5. ! 'A' IS THE MATRIX OF 'LENGTH' ROWS AND 'LENGTH' COLUMNS. ! IT MUST BE REAL AND SYMMETRIC. ! 'A' IS OVERWRITTEN DURING THE COMPUTATION. ! 'EPSILN' IS THE ACCEPTABLE ERROR IN EACH EIGENVECTOR COMPONENT; ! ITERATION STOPS WHEN THIS PRECISION IS REACHED. ! 'VALUES' IS A VECTOR OF LENGTH 'LENGTH' CONTAINING ORDERED ! EIGENVALUES, LARGEST TO SMALLEST. ! 'VECTRS' IS A MATRIX WHOSE COLUMNS ARE THE NORMALISED EIGENVECTORS ! OF 'A'. THE NUMBER OF THE COLUMN CORRESPONDS TO THE ORDER OF THE ! EIGEN VALUES IN 'VALUES'. ! 'NUMBER' IS THE NUMBER OF EIGENVALUES AND EIGENVECTORS DESIRED ! (ON INPUT) AND THE NUMBER ACTUALLY FOUND (ON OUTPUT). ! NDA AND NDV ARE SIZES OF THE DIMENSIONS OF THE FIRST (ROW) ! SUBSCRIPTS OF 'A' AND 'VECTRS' (RESPECTIVELY) ! IN THE CALLING PROGRAM ( >= LENGTH). ! NOTE: IF BY CHANCE THE LARGEST EIGENVALUE IS ASSOCIATED ! WITH AN EIGENVECTOR EXACTLY PERPENDICULAR TO (1,1,...,1), ! THE SECOND-LARGEST EIGENVALUE WILL BE FOUND BEFORE THE LARGEST. IMPLICIT NONE REAL*8, INTENT(IN) :: epsiln INTEGER, INTENT(IN) :: length, nda, ndv REAL*8, DIMENSION(nda, *), INTENT(IN OUT) :: a INTEGER, INTENT(IN OUT) :: number REAL*8, DIMENSION(*), INTENT(OUT) :: values REAL*8, DIMENSION(ndv, *), INTENT(OUT) :: vectrs INTEGER, PARAMETER :: maxitr = 1000 INTEGER :: ir, j, jc, jv, k LOGICAL :: ok REAL*8 :: lambda, size REAL*8, DIMENSION(:), ALLOCATABLE :: x ALLOCATE ( x(length) ) number = MIN(number, length) ! INITIALIZE TO PROVIDE FOR CASE OF SEVERAL ZERO EIGENVALUES DO jv = 1, number values(jv) = 0.0D0 DO ir = 1, length vectrs(ir, jv) = 0.0D0 END DO END DO ! LOOP, FOR EACH DESIRED EIGENVALUE... DO jv = 1, number ! INITIALIZE POWER METHOD WITH RANDOM UNIT VECTOR DO j = 1, length x(j) = 1.0D0 / DSQRT(1.0D0 * length) END DO ! ITERATE POWER METHOD UP TO MAXITR TIMES DO k = 1, maxitr ! MULTIPLY A TIMES X DO ir = 1, length vectrs(ir, jv) = 0.0 DO jc = 1, length vectrs(ir, jv) = vectrs(ir, jv) + a(ir, jc) * x(jc) END DO END DO ! FIND NEW ESTIMATE OF EIGENVALUE size = 0.0D0 DO ir = 1, length size = size + vectrs(ir, jv)**2 END DO lambda = DSQRT(size) ! NORMALIZE ESTIMATE OF EIGENVECTOR IF (size == 0.0D0) RETURN DO ir = 1, length vectrs(ir, jv) = vectrs(ir, jv) / lambda END DO ! TEST FOR CONVERGENCE ok = .TRUE. DO ir = 1, length ok = ok.AND.(ABS(vectrs(ir, jv) - x(ir)) <= epsiln) END DO IF (ok) GO TO 52 DO ir = 1, length x(ir) = vectrs(ir, jv) END DO END DO ! IF CONVERGENCE HAS FAILED TO OCCUR... number = jv - 1 WRITE (6, 51) jv, maxitr, number 51 FORMAT (' SUBROUTINE EIGEN FAILED TO FIND EIGEN VECTOR #', & & I5/' IN ',I10,' ATTEMPTS. RETURNING NUMBER=',I5) RETURN ! ANOTHER EIGENVALUE HAS BEEN FOUND 52 values(jv) = lambda ! SWEEP OUT KNOWN EIGENVALUE FROM MATRIX A DO ir = 1, length DO jc = 1, length a(ir, jc) = a(ir, jc) - & & lambda * vectrs(ir, jv) * vectrs(jc, jv) END DO END DO END DO DEALLOCATE ( x ) RETURN END SUBROUTINE DEigen SUBROUTINE Pause() !The purpose of this routine is to give the user time to read an error message !before the program stops (and the operating system closes its window). !The programmer should WRITE the error, CALL Pause(), and then STOP. IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause END PROGRAM Analyze_Velocity_Evolution