! ----------------------------------------------------------- ! ! THERMODYNAMIC FIELD ESTIMATOR MODULE ! ! ----------------------------------------------------------- ! MODULE ESTIMATOR ! ======================================================================= USE CG_MINIMIZER IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER, PRIVATE :: DIM,NP,TNB INTEGER, PRIVATE, ALLOCATABLE, DIMENSION (:,:) :: ORDER,KB INTEGER, PRIVATE, ALLOCATABLE, DIMENSION (:) :: ORDER_MAX,NB LOGICAL, PRIVATE :: ASSUME_NO_DRIFT, PRINT_INFO INTEGER MAXITER ! ----------------------------------------------------------------------- ! ORDER(VARIABLE,DIMENSION) IS THE HIGHEST ORDER OF SUB-BASIS FUNCTIONS. ! ORDER_MAX IS THE MAXIMUM ORDER OF SUB-BASIS FUNCTIONS WITH X,Y,Z. ! NB IS THE NUMBER OF BASIS FUNCTIONS FOR FIELD VARIABLES (T,VX,VY,VZ). ! KB(VARIABLE,1) AND KB(VARIABLE,2) ARE POINTERS IN BASIS. ! ----------------------------------------------------------------------- DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION (:,:) :: BASIS, SUB_BASIS DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION (:,:) :: BOUNDARY, MU DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION (:,:) :: PV, AA, BT DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION (:) :: BETA, PIVOT DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION (:) :: COEF, XU, XBASIS DOUBLE PRECISION, PRIVATE :: VALUE_LAST, LAST_VALUE ! ----------------------------------------------------------------------------- ! BASIS(PARTICLE,INDEX) ARE BASIS FUNCTION VALUES AT PARTICLE POSITIONS, ! INDEX ARE SECTIONED BY KB(VARIABLE,1) AND KB(VARIABLE,2) FOR EACH VARIABLE. ! BOUNDARY(DIMENSION,1|2) ARE THE MINIMUM|MAXIMUM COORDINATE OF THAT DIMENSION. ! SUB_BASIS(ORDER, DIMENSION) ARE SUB-BASIS FUNCTION VALUES AT CERTAIN ! PARTICLE POSITION (X, Y OR Z), ORDER IS FROM 1 TO ORDER_MAX(DIMENSION). ! MU(NP,DIMENSION) ARE REDUCED PARTICLE COORDINATES RANGING -1 TO 1. ! PV(NP,1|2|3) ARE PARTICLE VELOCITIES, PV(NP,0) IS VELOCITY (MINUS DRIFT) SQUARED. ! BETA(NP) ARE THE INVERSE TEMPERATURE EVALUATED AT PARTICLE POSITIONS. ! COEF(TNB) ARE THE BASIS FUNCTION EXPANSION COEFFICIENTS. ! ------------------------------------------------------------------------------ CONTAINS SUBROUTINE ESTIMATOR_INITIALIZE ( & PARTICLE_NUMBER, DIMENSION, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, & ORDER_T_X, ORDER_T_Y, ORDER_T_Z, ORDER_VX_X, ORDER_VX_Y, ORDER_VX_Z, & ORDER_VY_X, ORDER_VY_Y, ORDER_VY_Z, ORDER_VZ_X, ORDER_VZ_Y, ORDER_VZ_Z, P) ! ------------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H, O-Z) INTEGER PARTICLE_NUMBER INTEGER, OPTIONAL :: DIMENSION, & ORDER_T_X, ORDER_T_Y, ORDER_T_Z, ORDER_VX_X, ORDER_VX_Y, ORDER_VX_Z, & ORDER_VY_X, ORDER_VY_Y, ORDER_VY_Z, ORDER_VZ_X, ORDER_VZ_Y, ORDER_VZ_Z DOUBLE PRECISION, OPTIONAL :: XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX LOGICAL, OPTIONAL :: P DIM = 3 IF (PRESENT(DIMENSION)) DIM = DIMENSION ! THE FOUR FIELD VARIABLES (T,VX,VY,VZ) ARE NUMBERED 0 TO 3. ALLOCATE (ORDER(0:DIM,DIM)) ! DEFAULT IS A CONSTANT IN SPACE ORDER = 1 IF (PRESENT(ORDER_T_X)) ORDER(0,1) = ORDER_T_X IF (PRESENT(ORDER_T_Y)) ORDER(0,2) = ORDER_T_Y IF (PRESENT(ORDER_T_Z)) ORDER(0,3) = ORDER_T_Z IF (PRESENT(ORDER_VX_X)) ORDER(1,1) = ORDER_VX_X IF (PRESENT(ORDER_VX_Y)) ORDER(1,2) = ORDER_VX_Y IF (PRESENT(ORDER_VX_Z)) ORDER(1,3) = ORDER_VX_Z IF (PRESENT(ORDER_VY_X)) ORDER(2,1) = ORDER_VY_X IF (PRESENT(ORDER_VY_Y)) ORDER(2,2) = ORDER_VY_Y IF (PRESENT(ORDER_VY_Z)) ORDER(2,3) = ORDER_VY_Z IF (PRESENT(ORDER_VZ_X)) ORDER(3,1) = ORDER_VZ_X IF (PRESENT(ORDER_VZ_Y)) ORDER(3,2) = ORDER_VZ_Y IF (PRESENT(ORDER_VZ_Z)) ORDER(3,3) = ORDER_VZ_Z ALLOCATE (ORDER_MAX(DIM)) ! FIND OUT THE HIGHEST ORDER OF SUB-BASIS FUNCTIONS IN EACH DIMENSION: ORDER_MAX = MAXVAL(ORDER, DIM=1) ALLOCATE (SUB_BASIS(1:MAXVAL(ORDER_MAX), DIM)) ! A BUFFER THAT STORES THE TABLE OF SUB-BASIS FUNCTION VALUES & ! OF ALL DIMENSIONS AT A CERTAIN SPATIAL LOCATION. ALLOCATE (BOUNDARY(DIM,2)) ! DEFAULT IS 0. TO 1. BOUNDARY(:,1) = 0.D0 BOUNDARY(:,2) = 1.D0 IF (PRESENT(XMIN)) BOUNDARY(1,1) = XMIN IF (PRESENT(XMAX)) BOUNDARY(1,2) = XMAX IF (PRESENT(YMIN)) BOUNDARY(2,1) = YMIN IF (PRESENT(YMAX)) BOUNDARY(2,2) = YMAX IF (PRESENT(ZMIN)) BOUNDARY(3,1) = ZMIN IF (PRESENT(ZMAX)) BOUNDARY(3,2) = ZMAX NP = PARTICLE_NUMBER ALLOCATE (MU(NP,DIM)) ALLOCATE (PV(NP,0:DIM)) ALLOCATE (BETA(NP)) ! THE BASIS FUNCTIONS ARE DIRECT PRODUCTS OF ELEMENTARY ! SUB-BASIS FUNCTIONS IN ALL DIMENSIONS ALLOCATE (NB(0:DIM)) NB = PRODUCT(ORDER, DIM=2) TNB = SUM(NB) ! NB(VARIABLE) IS THE NUMBER OF BASIS FUNCTIONS FOR EACH VARIABLE. ! TNB IS THE TOTAL NUMBER OF BASIS FUNCTIONS FOR ALL VARIABLES COMBINED. ALLOCATE (BASIS(NP,TNB)) ALLOCATE (KB(0:DIM,2)) DO I = 0,DIM ! LOWER AND UPPER POINTERS IN BASIS FOR EACH FIELD VARIABLE. KB(I,1) = SUM(NB(0:I))-NB(I)+1 KB(I,2) = SUM(NB(0:I)) ENDDO ALLOCATE (COEF(TNB)) ! ALLOCATE BUFFER FOR VELOCITY FIELD COEFFICIENTS RELAXATION. MAXV = MAXVAL(NB(1:DIM)) ALLOCATE (BT(MAXV,NP)) ALLOCATE (AA(MAXV,MAXV)) ALLOCATE (PIVOT(MAXV)) ALLOCATE (XU(DIM)) ALLOCATE (XBASIS(TNB)) ! INITIALIZE CONJUGATE GRADIENT MINIMIZER: ONLY TEMPERATURE FIELD ! VARIABLES (BECAUSE IT'S NONLINEAR) NEED TO BE APPARANT. CALL CG_INITIALIZE(NB(0)) PRINT_INFO = .TRUE. IF (PRESENT(P)) PRINT_INFO = P END SUBROUTINE ESTIMATOR_INITIALIZE ! ------------------------------------------------------------------------- SUBROUTINE ESTIMATOR_FINALIZE ! ------------------------------------------------------------------------- DEALLOCATE (ORDER) DEALLOCATE (ORDER_MAX) DEALLOCATE (SUB_BASIS) DEALLOCATE (BOUNDARY) DEALLOCATE (MU) DEALLOCATE (PV) DEALLOCATE (BETA) DEALLOCATE (NB) DEALLOCATE (BASIS) DEALLOCATE (KB) DEALLOCATE (COEF) DEALLOCATE (BT) DEALLOCATE (AA) DEALLOCATE (PIVOT) DEALLOCATE (XU) DEALLOCATE (XBASIS) CALL CG_FINALIZE END SUBROUTINE ESTIMATOR_FINALIZE ! ------------------------------------------------------------------------- SUBROUTINE ESTIMATE_FIELD (PXX,PX,PY,PZ,PVV,PVX,PVY,PVZ,PV2,TOL,MAXI) ! ------------------------------------------------------------------------- ! INPUT PARTICLE POSITIONS: USE SEPERATE ARRAYS (PX,PY,PZ) OR COMBINED PXX. ! INPUT PARTICLE VELOCITIES: USE SEPERATE (PVX,PVY,PVZ) OR COMBINED PVV. ! IF INPUT VELOCITY SQUARED PV2, THAT IMPLIES YOU ONLY FIT TEMPERATURE. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION, OPTIONAL, DIMENSION (:), INTENT(IN) :: & PXX,PX,PY,PZ,PVV,PVX,PVY,PVZ,PV2 ! PARTICLE DATA DOUBLE PRECISION, OPTIONAL, INTENT(IN) :: TOL INTEGER, OPTIONAL, INTENT(IN) :: MAXI DOUBLE PRECISION :: ACCURACY MU = 0.D0 IF (PRESENT(PXX)) THEN DO K = 1,DIM MU(:,K) = & -1+2*(PXX((K-1)*NP+1:K*NP)-BOUNDARY(K,1))/(BOUNDARY(K,2)-BOUNDARY(K,1)) ENDDO ELSE IF (PRESENT(PX)) & MU(:,1) = -1+2*(PX(1:NP)-BOUNDARY(1,1))/(BOUNDARY(1,2)-BOUNDARY(1,1)) IF (PRESENT(PY)) & MU(:,2) = -1+2*(PY(1:NP)-BOUNDARY(2,1))/(BOUNDARY(2,2)-BOUNDARY(2,1)) IF (PRESENT(PZ)) & MU(:,3) = -1+2*(PZ(1:NP)-BOUNDARY(3,1))/(BOUNDARY(3,2)-BOUNDARY(3,1)) ENDIF DO N = 1,NP CALL EVALUATE_BASIS(MU(N,:),BASIS(N,:)) ENDDO PV = 0.D0 COEF = 0.D0 IF (PRESENT(PV2)) THEN PV(:,0) = PV2(1:NP) ASSUME_NO_DRIFT = .TRUE. ELSE IF (PRESENT(PVV)) THEN DO K = 1,DIM PV(:,K) = PVV((K-1)*NP+1:K*NP) ENDDO ELSE IF (PRESENT(PVX)) PV(:,1) = PVX(1:NP) IF (PRESENT(PVY)) PV(:,2) = PVY(1:NP) IF (PRESENT(PVZ)) PV(:,3) = PVZ(1:NP) ENDIF ! STORE THE AVERAGE VELOCITY IN THE LOWEST BASIS COEFFICIENT DO K = 1,DIM COEF(KB(K,1)) = SUM(PV(:,K))/NP PV(:,0) = PV(:,0) + (PV(:,K)-COEF(KB(K,1)))**2 ENDDO ASSUME_NO_DRIFT = .FALSE. ENDIF ACCURACY = 1D-6 IF (PRESENT(TOL)) ACCURACY = TOL MAXITER = 1000 IF (PRESENT(MAXI)) MAXITER = MAXI ! ASSIGN INITIAL VALUES TO THE TEMPERATURE EXPANSION COEFFICIENTS: ! HERE JUST TAKE THE AVERAGE ENERGY AND PUT INTO THE LOWEST ORDER. COEF(KB(0,1)) = NP*DIM/SUM(PV(:,0))/(SUM(BASIS(:,KB(0,1)))/NP) ! INVOKE CONJUGATE GRADIENT MINIMIZER CALL CG_MINIMIZE(COEF(1:NB(0)),NB(0),ACCURACY,ITER,FMIN,VALUE,FORCE) IF (PRINT_INFO) THEN PRINT *,' ' PRINT *,'MINIMIZATION CONVERGED AFTER',ITER,' FORCE EVALUATIONS.' PRINT *,' ' ENDIF END SUBROUTINE ESTIMATE_FIELD ! ------------------------------------------------------------------------- SUBROUTINE EVALUATE_BASIS (X,U) ! ------------------------------------------------------------------------- ! EVALUATE BASIS FUNCTIONS AT A CERTAIN POINT: INPUT IS POSITION X(DIM) ! STORED IN REDUCED COORDINATES RANGING FROM -1 TO 1. OUTPUT IS U(TNB). DOUBLE PRECISION X(:),U(:) INTEGER INDICES(0:3),K,L ! BUILD UP SUB-BASIS FUNCTION TABLE: HERE WE USE CHEBYSHEV POLYNOMIALS. DO K = 1,DIM CALL CHEBYSHEV (X(K),ORDER_MAX(K),SUB_BASIS(:,K)) ENDDO DO L = 1, TNB U(L) = 1.D0 CALL ONE_TO_MANY (L,INDICES) ! PRINT *, 'L = ', L, ' INDICES = ', INDICES(0:DIM) DO K = 1,DIM U(L) = U(L)*SUB_BASIS(INDICES(K),K) ENDDO ENDDO END SUBROUTINE EVALUATE_BASIS ! ------------------------------------------------------------------------- SUBROUTINE ONE_TO_MANY (L,INDICES) ! ------------------------------------------------------------------------- ! TRANSFORM BASIS FUNCTION INDEX L TO SUB-BASIS FUNCTION INDICES ! INDICES(1..DIM). INDICES(0) IS THE VARIABLE NUMBER RANGING 0 TO 3. INTEGER INDICES(0:3), SUM, K, PRO INDICES(0) = 0 DO WHILE (L.GT.KB(INDICES(0),2)) INDICES(0) = INDICES(0)+1 ENDDO SUM = L-KB(INDICES(0),1) DO K = 1,DIM PRO = PRODUCT(ORDER(INDICES(0),K+1:DIM)) INDICES(K) = SUM/PRO SUM = SUM-INDICES(K)*PRO INDICES(K) = INDICES(K)+1 ENDDO END SUBROUTINE ONE_TO_MANY ! ------------------------------------------------------------------------- SUBROUTINE CHEBYSHEV (X,MAXORDER,T) ! ------------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION X, T(:) T(1) = 1.D0 IF (MAXORDER.GT.1) THEN T(2) = X DO L = 3,MAXORDER T(L) = 2*X*T(L-1)-T(L-2) ENDDO ENDIF END SUBROUTINE CHEBYSHEV ! ------------------------------------------------------------------------- FUNCTION VALUE (COEF_T) ! ------------------------------------------------------------------------- ! THE OBJECTIVE FUNCTION OF MINIMIZATION: ONLY NEED TO PASS TEMPERATURE ! BASIS COEFFICIENTS. THE VELOCITY BASIS COEFFICIENTS ARE AUTOMATICALLY ! RELAXED BY SOLVING A MATRIX EQUATION AND ARE STORED IN ! COEF(KB(1,1):KB(DIM,2)). IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (TINY=1D-20, HUGE=1D+20) ! DOUBLE PRECISION, DIMENSION(NB(0)) :: COEF_T DOUBLE PRECISION, DIMENSION(*) :: COEF_T DOUBLE PRECISION VALUE,SUMM,PRO,BETA_AVG INTEGER ITOT SAVE ITOT DATA ITOT /0/ IF (ITOT.GE.MAXITER) THEN VALUE = VALUE_LAST RETURN ENDIF BETA(1:NP) = MATMUL(BASIS(:,KB(0,1):KB(0,2)),COEF_T(1:NB(0))) ! TO PREVENT NEGATIVE TEMPERATURE ! BETA = ABS(BETA) IF (.NOT.ASSUME_NO_DRIFT) THEN DO K = 1,DIM BT(1:NB(K),:) = TRANSPOSE(BASIS(:,KB(K,1):KB(K,2))) DO N = 1,NP BT(1:NB(K),N) = BT(1:NB(K),N)*BETA(N) ENDDO AA(1:NB(K),1:NB(K)) = MATMUL(BT(1:NB(K),:),BASIS(:,KB(K,1):KB(K,2))) ! CHOLESKY SOLVER OF LINEAR EQUATIONS FOR SYMMETRIC AND ! POSITIVE DEFINITE MATRIX: AA*W=B CALL CHOLDC (AA,NB(K),PIVOT) CALL CHOLSL (AA,NB(K),PIVOT, & MATMUL(BT(1:NB(K),:),PV(:,K)),COEF(KB(K,1):KB(K,2))) ENDDO PV(:,0) = 0.D0 DO K = 1, DIM PV(:,0) = PV(:,0) + & (PV(:,K)-MATMUL(BASIS(:,KB(K,1):KB(K,2)),COEF(KB(K,1):KB(K,2))))**2 ENDDO ENDIF BETA_AVG = SUM(BETA)/NP SUMM = DOT_PRODUCT(BETA,PV(:,0)) - DIM*NP*LOG(BETA_AVG) BETA = BETA / BETA_AVG PRO = 1.D0 DO N =1,NP PRO = BETA(N)*PRO IF (PRO.LE.TINY.OR.PRO.GT.HUGE) THEN SUMM = SUMM - DIM*LOG(PRO) PRO = 1.D0 ENDIF ENDDO SUMM = SUMM - DIM*LOG(PRO) VALUE = SUMM/NP VALUE_LAST = VALUE ITOT = ITOT + 1 IF (PRINT_INFO) & PRINT *,'At',ITOT,'th function evaluation, the error =',VALUE END FUNCTION VALUE ! ------------------------------------------------------------------------- SUBROUTINE CHOLDC (A,N,P) ! ------------------------------------------------------------------------- integer n,i,j,k double precision a(:,:),p(:),summ do 13 i=1,n do 12 j=i,n summ=a(i,j) do 11 k=i-1,1,-1 summ=summ-a(i,k)*a(j,k) 11 continue if(i.eq.j)then if(summ.le.0.d0) stop 'choldc failed' p(i)=dsqrt(summ) else a(j,i)=summ/p(i) endif 12 continue 13 continue END SUBROUTINE CHOLDC ! ------------------------------------------------------------------------- SUBROUTINE CHOLSL (A,N,P,B,X) ! ------------------------------------------------------------------------- integer n,i,k double precision a(:,:),b(:),p(:),x(:),summ do 12 i=1,n summ=b(i) do 11 k=i-1,1,-1 summ=summ-a(i,k)*x(k) 11 continue x(i)=summ/p(i) 12 continue do 14 i=n,1,-1 summ=x(i) do 13 k=i+1,n summ=summ-a(k,i)*x(k) 13 continue x(i)=summ/p(i) 14 continue END SUBROUTINE CHOLSL ! ------------------------------------------------------------------------- SUBROUTINE FORCE (COEF_T,F) ! ------------------------------------------------------------------------- ! NUMERICAL DIFFERENTIATION OF THE OBJECTIVE FUNCTION TO GET THE FORCE: IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION, PARAMETER :: ALPHA=0.001 ! DOUBLE PRECISION COEF_T(NB(0)),F(NB(0)),LAST_VALUE DOUBLE PRECISION COEF_T(*),F(*) IF (PRINT_INFO) PRINT *, 'FORCE EVALUATION STARTS...' ! IN CG_MINIMIZE, BEFORE CALLING FORCE(), IT CALLS VALUE(), SO ! VALUE_LAST SHOULD BE EVALUATED AT THE POINT. LAST_VALUE = VALUE_LAST ! DUE TO A COMPILER BUG ON SGI WITH -O3 OPTION, THIS LINE MUST BE KEPT. VALUE_LAST = LOG(EXP(VALUE_LAST)) ! DUE TO A COMPILER BUG ON SGI WITH -O3 OPTION, THIS LINE MUST BE KEPT. DO I = 1,NB(0) ! if (abs(coef_T(i)).lt.alpha*coef_T(1)) then DV = ALPHA*COEF_T(1) ! else ! dv = alpha*coef_T(i) ! endif COEF_T(I) = COEF_T(I) + DV F(I) = -(VALUE(COEF_T)-LAST_VALUE)/DV COEF_T(I) = COEF_T(I) - DV ENDDO IF (PRINT_INFO) PRINT *, 'FORCE EVALUATION FINISHES.' END SUBROUTINE FORCE ! ------------------------------------------------------------------------- SUBROUTINE GIVE_ESTIMATE (X,Y,Z,XX,T,VX,VY,VZ,VV) ! ------------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION, OPTIONAL, INTENT(IN) :: X,Y,Z,XX(:) DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: T,VX,VY,VZ,VV(:) XU = 0.D0 IF (PRESENT(XX)) THEN XU = -1+2*(XX(1:DIM)-BOUNDARY(:,1))/(BOUNDARY(:,2)-BOUNDARY(:,1)) ELSE IF (PRESENT(X)) & XU(1) = -1+2*(X-BOUNDARY(1,1))/(BOUNDARY(1,2)-BOUNDARY(1,1)) IF (PRESENT(Y)) & XU(2) = -1+2*(Y-BOUNDARY(2,1))/(BOUNDARY(2,2)-BOUNDARY(2,1)) IF (PRESENT(Z)) & XU(3) = -1+2*(Z-BOUNDARY(3,1))/(BOUNDARY(3,2)-BOUNDARY(3,1)) ENDIF CALL EVALUATE_BASIS(XU,XBASIS) IF (PRESENT(T)) & T = 1/DOT_PRODUCT(XBASIS(KB(0,1):KB(0,2)),COEF(KB(0,1):KB(0,2))) IF (PRESENT(VX)) & VX = DOT_PRODUCT(XBASIS(KB(1,1):KB(1,2)),COEF(KB(1,1):KB(1,2))) IF (PRESENT(VY)) & VY = DOT_PRODUCT(XBASIS(KB(2,1):KB(2,2)),COEF(KB(2,1):KB(2,2))) IF (PRESENT(VZ)) & VZ = DOT_PRODUCT(XBASIS(KB(3,1):KB(3,2)),COEF(KB(3,1):KB(3,2))) IF (PRESENT(VV)) THEN DO K = 1, DIM VV(K) = DOT_PRODUCT(XBASIS(KB(K,1):KB(K,2)),COEF(KB(K,1):KB(K,2))) ENDDO ENDIF END SUBROUTINE GIVE_ESTIMATE ! ------------------------------------------------------------------------- END MODULE ESTIMATOR ! =============================================================================