! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE DYNAMIC_MODEL (X, I_BC, NODES, N_RES_EQ, COEF_C_EQ, & NDEX_C_EQ, MASSES, DAMPENER, CC_IN) !b NDEX_C_EQ, CC_IN) ! ****************************************************** ! * MODULAR PROGRAMS FOR FINITE ELEMENT ANALYSES * ! * VERSION 1, DYNAMIC ANALYSIS SOLUTION * ! * COPYRIGHT J. E. AKIN, 2004 * ! ****************************************************** ! Fortran 90 version, uses double precision and "automatic arrays" ! ! A SET OF BUILDING BLOCK PROGRAMS ! BY ! DR. J. E. AKIN, P.E. ! DEPT. OF MECHANICAL ENGINEERING & MATERIALS SCIENCE ! RICE UNIVERSITY ! HOUSTON, TEXAS 77251-1892 ! ! email: akin@rice.edu ! Use System_Constants ! for DD_DOT (N_D_FRE), DD_2DOT (N_D_FRE) etc Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_FREE), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_FREE, LT_QP), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! EL_M (LT_FREE, LT_FREE), DIAG_M (LT_FREE), ! D (LT_FREE), ELEM_NODES (LT_N) Use Sys_Properties_Data Use GEOMETRIC_PROPERTIES IMPLICIT NONE LOGICAL :: FACT, BACK CHARACTER (LEN=3) :: T_STRING ! for file names CHARACTER (LEN=21) :: N_NAME ! file name CHARACTER (LEN=20) :: D_NAME ! file name ! Given Arrays REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP), INTENT(IN) :: COEF_C_EQ (MAX_ACT, N_CEQ) !b REAL(DP), INTENT(IN) :: MASSES (N_D_FRE) !b ! point masses REAL(DP), INTENT(IN) :: DAMPENER (N_D_FRE) !b ! point damping REAL(DP), INTENT(IN) :: CC_IN (IN_RHS) !b forcing amplitude INTEGER, INTENT(IN) :: I_BC (MAX_NP) INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER, INTENT(IN) :: N_RES_EQ (MAX_TYP) INTEGER, INTENT(IN) :: NDEX_C_EQ (MAX_ACT, N_CEQ) !b ! Allocatable Arrays INTEGER, ALLOCATABLE :: L_TO_L_NEIGH (:, :) ! neighbors INTEGER, ALLOCATABLE :: L_TO_N_NEIGH (:, :) ! neighbors ! Allocatable Arrays REAL(DP), ALLOCATABLE :: M_SKY (:) ! symmetric mass REAL(DP), ALLOCATABLE :: M_SKY_LOWER (:) ! symmetric damping !b REAL(DP), ALLOCATABLE :: S_SKY (:) ! symmetric part REAL(DP), ALLOCATABLE :: S_SKY_LOWER (:) ! unsymmetric part INTEGER, ALLOCATABLE :: I_DIAG (:), I_DOF_HI (:) ! sky parts ! Automatic Arrays !b REAL(DP) :: C_MD (N_D_FRE) !b C_MD == M_SKY * D_OLD REAL(DP) :: CC (N_D_FRE) !b, CC_OLD (N_D_FRE) !b REAL(DP) :: DD (N_D_FRE), DD_OLD (N_D_FRE) !b REAL(DP) :: SCP_AVERAGES (MAX_NP, SCP_FIT) !b REAL(DP) :: EL_ENORM (N_ELEMS), ERR_MAX !b REAL(DP) :: EL_REFINE (N_ELEMS) !b INTEGER :: NULL_BC (N_CEQ), N_BC ! currently only null EBC INTEGER :: HIST_INDEX (N_G_DOF) !b REAL(DP) :: WILSON_THETA = 1.25d0, HUGHES_ALPHA = -0.05d0 !b REAL(DP) :: NEWMARK_BETA = 0.25d0, NEWMARK_GAMMA = 0.5d0 !b REAL(DP) :: MASS_DAMPING = 0.d0, STIF_DAMPING = 0.d0 REAL(DP) :: W_SQ_MAX, W_MAX, DT INTEGER :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER :: L_TO_L_SUM (N_ELEMS) INTEGER :: L_TO_N_SUM (MAX_NP) INTEGER, SAVE :: BAR_DOF = 1 ! limit bar chart loop INTEGER :: ITER, IBAR, IPARM, IE, OLD_DEL = -1 INTEGER :: IO_1, I, I_WARN, J, K, METHOD, N_STEPS, I_PRINT INTEGER :: T_G, T_S ! TIME GROUP AND TIME STEP NUMBERS !b REAL(DP) :: R_TEST, RATIO, OLD_NORM, NEW_NORM, DIFF_NORM ! SEE NOTATION.F FOR VARIABLE AND ARRAY NOTATIONS ! METHOD: 1-Wilson, 2-Newmark, 3-Hughes, 4-Trapezoidal ! N_BC = NUMBER OF (NULL) ESSENTUAL BOUNDARY CONDITIONS ! NULL_BC = LIST OF DOF NUMBERS HAVING NULL EBC METHOD = TIME_METHOD ! Global control DT = TIME_STEP ! Global control N_STEPS = T_STEPS ! Global control SELECT CASE (METHOD) ! titles done below CASE (1) ; WILSON_METHOD = .TRUE. CASE (2) ; NEWMARK_METHOD = .TRUE. CASE (3) ; HUGHES_METHOD = .TRUE. CASE DEFAULT ! trapezoidal WILSON_METHOD = .FALSE. NEWMARK_METHOD = .FALSE. HUGHES_METHOD = .FALSE. END SELECT IF ( I_BUG > 0 ) THEN PRINT *,'=== USING DYNAMIC_MODEL ===', METHOD !b PRINT *,'MAX_XYZ size ', SIZE (MAX_XYZ) END IF CALL ALLOCATE_GEOM_PROP ALLOCATE ( I_DOF_HI (N_D_FRE), I_DIAG (N_D_FRE) ) IF ( I_BUG > 0 ) PRINT *, 'ALLOCATED SKYLINE INTEGER ARRAYS' I_DIAG = 0 ; I_DOF_HI = 0 ! GENERATE & LIST MESH NEIGHBORS (PATCH) INFORMATION IF ( SCP_NEIGH_FACE ) CALL SET_L_SHAPE_FACE_NEEDS ! for NEEDS !s I_WARN = 1 CALL FIRST_LAST_L_AT_PT (NODES, L_FIRST, L_LAST, I_WARN) !b NEEDS = 1 CALL COUNT_ELEMS_AT_ELEM (N_ELEMS, NOD_PER_EL, MAX_NP, L_FIRST, & L_LAST, NODES, NEEDS, L_TO_L_SUM, N_WARN) NEIGH_L = MAXVAL (L_TO_L_SUM) NEIGH_P = NEIGH_L ! IF ELEMENT BASED PATCHES !b ALLOCATE ( L_TO_L_NEIGH (NEIGH_L, N_ELEMS) ) CALL FORM_ELEMS_AT_EL (N_ELEMS, NOD_PER_EL, MAX_NP, & L_FIRST, L_LAST, NODES, N_SPACE, & L_TO_L_SUM, L_TO_L_NEIGH, & NEIGH_L, NEEDS, ON_BOUNDARY, U_ON_B) !b !b L_TO_L_SUM, L_TO_L_NEIGH, NEIGH_L) MAX_FACES = 2 !b to avoid warning from compiler IF ( LIST_EL_TO_EL ) CALL LIST_ELEMENTS_AT_EL (L_TO_L_NEIGH) CALL COUNT_L_AT_NODE (N_ELEMS, NOD_PER_EL, MAX_NP, NODES, & L_TO_N_SUM) IF ( I_BUG > 0 .OR. PT_EL_SUMS .OR. PT_EL_LIST ) THEN CALL LIST_L_AT_NODE_INFO (L_TO_N_SUM, L_FIRST, L_LAST) IF ( PT_EL_LIST .OR. PT_EL_LIST ) THEN NEIGH_N = MAXVAL (L_TO_N_SUM) NEIGH_P = NEIGH_N ! SINCE NODAL BASED PATCHES !b ALLOCATE ( L_TO_N_NEIGH (NEIGH_N, MAX_NP ) ) CALL FORM_L_ADJACENT_NODES (N_ELEMS, NOD_PER_EL, MAX_NP, & NODES, NEIGH_N, L_TO_N_NEIGH) IF ( PT_EL_LIST ) CALL LIST_ELEMENTS_AT_PT (L_TO_N_NEIGH) DEALLOCATE ( L_TO_N_NEIGH ) !s ? END IF END IF ! BUILD NEIGHBORS, TO DEBUG ! BUILD MESH GEOMETRIC PROPERTIES (INERTIA MATRIX) CALL ASSEMBLE_MESH_GEOM (NODES, X) !--> COMPUTE SYSTEM MATRIX STORAGE SIZE ! ! DETERMINE SYSTEM SKYLINE STORAGE CALL MESH_SKYLINE (NODES, I_DOF_HI, I_DIAG) IF ( I_BUG > 0 ) PRINT *, 'max I_DIAG ', I_DIAG (N_D_FRE) ! CHANGE COLUMN HEIGHTS FOR LINEAR CONSTRAINTS IF ( MAX_ACT > 1 ) CALL CEQ_SKYLINE (N_RES_EQ, NDEX_C_EQ, & I_DOF_HI, I_DIAG) N_COEFF = I_DIAG (N_D_FRE) NOT_SYM = 0 ; N_STORE = N_COEFF IF ( .NOT. SYMMETRIC ) THEN NOT_SYM = N_COEFF ; N_STORE = N_COEFF*2 END IF !b WRITE (6, 5) N_D_FRE !b 5 FORMAT (/, 'TOTAL NUMBER OF SYSTEM EQUATIONS ....',I10 ) ! ** STOP IF DATA CHECK ONLY RUN ** IF ( INPUT_ONLY ) THEN IF ( N_WARN == 0 ) THEN PRINT *, 'END OF MODEL_F90 DATA CHECK (NO WARNINGS)' ELSE PRINT *, 'END MODEL_F90 DATA CHECK, WITH ', N_WARN,' WARNINGS' END IF STOP 'REQUESTED input_only CHECK COMPLETED. EXITING' END IF ! Data check only ! ALLOCATE AND ZERO THE SYSTEM MATRICES M_COEFF = N_COEFF ALLOCATE ( M_SKY(M_COEFF) ) ! add status M_SKY = 0.d0 ! zero mass ALLOCATE ( S_SKY(N_COEFF) ) ! add status S_SKY = 0.d0 ! zero stiffness !b IF ( .NOT. SYMMETRIC ) THEN IF ( SYMMETRIC ) THEN ALLOCATE ( S_SKY_LOWER (1) ) ! Fake lower stiffness array ALLOCATE ( M_SKY_LOWER (N_COEFF) ) ! Fake damping array ELSE ! non-symmetric PRINT *,'WARNING: UN_SYMMETRIC MATRIX IN DYNAMIC_MODEL' N_WARN = N_WARN + 1 !b STOP 'DYNAMIC MODEL CURRENTLY INVALID FOR UN_SYMMETRIC MATRIX' ALLOCATE ( S_SKY_LOWER (NOT_SYM) ) ! add status ALLOCATE ( M_SKY_LOWER (N_COEFF) ) ! add status END IF ! additional arrays S_SKY_LOWER = 0.d0 ! zero stiffness M_SKY_LOWER = 0.d0 IF ( I_BUG > 0 ) PRINT *, 'ALLOCATED SKYLINE REAL ARRAYS' ! INITIALIZE THE SYSTEM SOURCE MATRIX DD = 0.d0 ; CC = 0.d0 ! zero load ! *** BEGIN TIME INTEGRATION GROUPS *** !b CC_OLD = 0.d0 !b IF ( RESTART_STEP > 0 ) THEN ! it governs IC CALL GET_STEP_RESULTS (DD_OLD) ! for RESTART_STEP THIS_STEP = RESTART_STEP IF ( SAVE_1248 ) NEXT_SAVE = 2 * RESTART_STEP ELSE IF ( INITIAL_FUNCTION ) THEN ! spatial function PRINT *,'NOTE: USING SET_INITIAL_CONDITION ROUTINE' CALL SET_INITIAL_CONDITION (X, DD_OLD) ELSE ! constant Initial Conditions PRINT *,'NOTE: DYNAMIC_MODEL INITIAL_VALUE =',INITIAL_VALUE DD_OLD = INITIAL_VALUE END IF ! constant value THIS_STEP = 0 CALL SAVE_STEP_RESULTS (DD_OLD) ! to node_results_000.tmp END IF ! restart step > 0 DD = DD_OLD ! RECOVER TIME INTEGRATION DATA FOR EACH GROUP START_TIME = START_TIME !b GROUPS: DO T_G = 1, T_SETS !b T_STEPS = T_STEPS !b TIME_STEP = TIME_STEP !b TS_B = TS_B ! INITIALIZE THIS TIME GROUP PRINT *, ' ' PRINT *, '*** BEGINNING TIME GROUP NUMBER ', T_G, ' ***' PRINT *, 'STARTING TIME ..............= ', START_TIME PRINT *, 'RESTART FROM STEP NUMBER ...= ', RESTART_STEP PRINT *, 'NUMBER OF TIME STEPS........= ', T_STEPS PRINT *, 'TIME STEP SIZE .............= ', TIME_STEP PRINT *, ' ' ! Open a time-history graph file IF ( HISTORY_DOF > 0 ) THEN T_STRING = TO3STR (RESTART_STEP+1) ! integer to string D_NAME = 'dof_history_' // T_STRING // '.tmp' OPEN (U_DOF, FILE = D_NAME, ACTION='WRITE', STATUS='REPLACE', & IOSTAT = IO_1) IF ( IO_1 > 0 ) STOP 'NOTE: FILE OPEN U_DOF FAILED' ! Insert initial values WRITE (U_DOF, '(I6, 2(1PE13.5))') HISTORY_DOF, START_TIME, & DD (HISTORY_DOF) END IF IF ( HISTORY_NODE > 0 ) THEN T_STRING = TO3STR (RESTART_STEP+1) ! integer to string N_NAME = 'node_history_' // T_STRING // '.tmp' OPEN (U_NODE, FILE = N_NAME, ACTION='WRITE', STATUS='REPLACE', & IOSTAT = IO_1) IF ( IO_1 > 0 ) STOP 'NOTE: FILE OPEN U_NODE FAILED' ! Insert initial values HIST_INDEX = GET_INDEX_AT_PT (HISTORY_NODE) WRITE (U_NODE, '(I6, (8(1PE13.5)) )') HISTORY_NODE, START_TIME, & DD (HIST_INDEX) END IF !--> *** CALCULATE AND ASSEMBLE ELEMENT MATRICES *** !--> *** GENERATE POST SOLUTION MATRICES & STORE *** !b print *,'DD_OLD ' ; call rprint (DD_OLD, LT_FREE, 1, 1) CC = 0.0_DP ! zero load CALL ASSEMBLE_SKYLINE_EQS (I_DIAG, NODES, S_SKY, CC, X, & DD_OLD, S_SKY_LOWER, ITER, M_SKY, M_SKY_LOWER) IF ( IN_RHS > 0 ) CC = CC + CC_IN ! CONSTANT NODAL SOURCES IF ( DEBUG_ASSEMBLY ) CALL SUM_RHS (CC) !b CC_OLD = CC !b IF ( PT_MASSES ) THEN ! Add diagonal mass terms DO J = 1, N_D_FRE M_SKY ( I_DIAG (J) ) = M_SKY ( I_DIAG (J) ) + MASSES (J) END DO END IF ! point masses were input IF ( PT_DAMPING ) THEN ! Add diagonal dampener terms DO J = 1, N_D_FRE M_SKY_LOWER (I_DIAG (J)) = M_SKY_LOWER (I_DIAG (J))+DAMPENER (J) END DO END IF ! point masses were input ! *** ASSEMBLY COMPLETED, CHECK SOURCES *** IF ( DEBUG_ASSEMBLY .AND. N_D_FRE < M_SHOW .AND. T_G == 1 ) THEN PRINT *,'DEBUG ASSEMBLY IN FULL MATRIX FORM: STIFFNESS' CALL SHOW_SKY_FUL (I_DIAG) ! show symbols CALL SKY_FUL_LOCATIONS (I_DIAG) ! show locations CALL SKY_FUL_GEN (S_SKY, I_DIAG, CC, S_SKY_LOWER) ! show values CALL SUM_RHS (CC) PRINT *,'MASS MATRIX' CALL SKY_FUL_SYM (M_SKY, I_DIAG, CC) ! show values PRINT *,'DAMPING MATRIX' CALL SKY_FUL_SYM (M_SKY_LOWER, I_DIAG, CC) ! show values END IF CALL ESTIMATE_SYSTEM_MAX_FREQ (M_SKY, S_SKY, I_DIAG, W_SQ_MAX) W_MAX = DSQRT (W_SQ_MAX) PRINT *,'MAX SYSTEM FREQUENCY = ', W_MAX PRINT *,'STABLE EXPLICIT TIME STEP > ', 2 / W_MAX PRINT *,'CURRENT IMPLICIT TIME STEP = ', DT PRINT *,'RATIO (5-20 COMMON) = ', DT * 0.5d0 * W_MAX ! Allow only null ebc at this time N_BC = SIZE ( NDEX_C_EQ (1, :) ) NULL_BC (1:N_BC) = NDEX_C_EQ (1, 1:N_BC) I_PRINT = 1 !b XXX control this IF ( WILSON_METHOD ) THEN CALL WILSON_METHOD_SYM (I_DIAG, M_SKY, M_SKY_LOWER, S_SKY, & CC_IN, DD, DD_DOT, DT, N_STEPS, I_PRINT, & N_BC, NULL_BC) ELSE ! Newmark, or Hughes, or trapezoidal CALL NEWMARK_METHOD_SYM (I_DIAG, M_SKY, M_SKY_LOWER, S_SKY, & CC_IN, DD, DD_DOT, DT, N_STEPS, I_PRINT,& N_BC, NULL_BC) END IF END DO GROUPS ! over time groups !b XXX deallocate arrays here IF ( N_WARN > 0 ) PRINT *,'DYNAMIC_MODEL GAVE ', N_WARN, ' WARNINGS' END SUBROUTINE DYNAMIC_MODEL SUBROUTINE WILSON_METHOD_SYM (I_DIAG, A_SKY, B_SKY, C_SKY, P_IN, R, & DR, DT, N_STEPS, I_PRINT, N_BC, NULL_BC) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! WILSON THETA UNCONDITIONALLY STABLE METHOD OF ! STEP BY STEP INTEGRATION OF MATRIX EQUATIONS: ! A*D2R(T)/DT2 + B*DR(T)/DT + C*R(T) = P(T) ! INITIAL VALUES OF R AND DR ARE PASSED THRU ARGUMENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_COEFF, N_D_FRE IMPLICIT NONE INTEGER, INTENT (IN) :: N_STEPS, I_PRINT, N_BC INTEGER, INTENT (IN) :: I_DIAG (N_D_FRE) INTEGER, INTENT (IN) :: NULL_BC (N_BC) REAL(DP), INTENT (IN) :: DT REAL(DP), INTENT (INOUT) :: A_SKY (N_COEFF), C_SKY (N_COEFF) REAL(DP), INTENT (INOUT) :: B_SKY (N_COEFF) REAL(DP), INTENT (IN) :: P_IN (N_D_FRE) ! source amplitude REAL(DP), INTENT (INOUT) :: R (N_D_FRE), DR (N_D_FRE) REAL(DP), PARAMETER :: ZERO = 0.d0 INTEGER :: I, K, I_STEP, I_COUNT REAL(DP) :: THETA, DELT, TAU, T, T_PLUS, SCALE_1, SCALE_2 !b REAL(DP) :: WILSON_THETA = 1.25d0, HUGHES_ALPHA = -0.05d0 !b REAL(DP) :: NEWMARK_BETA = 0.25d0, NEWMARK_GAMMA = 0.5d0 !b REAL(DP) :: MASS_DAMPING = 0.d0, STIF_DAMPING = 0.d0 ! Automatic Arrays REAL(DP) :: DR_PLUS (N_D_FRE), D2R (N_D_FRE), P (N_D_FRE) REAL(DP) :: WORK_1 (N_D_FRE), WORK_2 (N_D_FRE) ! N_COEFF = TOTAL NUMBER OF TERMS IN SKYLINE ! N_BC = NUMBER OF D.O.F. WITH SPECIFIED VALUES OF ZERO ! NULL_BC = ARRAY CONTAINING THE N_BC DOF NUMBERS WITH ZERO BC ! R,DR,D2R = 0,1,2 ORDER DERIV. OF R W.R.T. T AT TIME=T ! DR_PLUS = VALUE OF DR AT TIME = T + DELT ! N_STEPS = NO. OF INTEGRATION STEPS ! I_PRINT = NO. OF INTEGRATION STEPS BETWEEN PRINTING ! ** INITIAL CALCULATIONS ** WRITE (6, '( "WILSON STEP BY STEP INTEGRATION",//)') THETA = WILSON_THETA DELT = (THETA - 1.d0) * DT TAU = THETA * DT SCALE_1 = 1.d0 - 1.d0 / THETA SCALE_2 = 1.d0 - 2.d0 / THETA IF ( N_BC < 1 ) PRINT *,'NOTE: NO CONSTRAINTS IN WILSON_METHOD' P = 0.d0 ; DR_PLUS = 0.d0 ; D2R = 0.d0 WORK_1 = 0.d0 ; WORK_2 = 0.d0 IF ( I_BUG > 0 ) THEN PRINT *,'DEBUG: A, B, C' CALL SKY_FUL_SYM (A_SKY, I_DIAG, P_IN) CALL SKY_FUL_SYM (B_SKY, I_DIAG, P_IN) CALL SKY_FUL_SYM (C_SKY, I_DIAG, P_IN) END IF IF ( .NOT. EL_DAMPING ) THEN !b B_SKY = 0.d0 IF ( MASS_DAMPING > 0.d0 ) B_SKY = B_SKY + MASS_DAMPING * A_SKY IF ( STIF_DAMPING > 0.d0 ) B_SKY = B_SKY + STIF_DAMPING * C_SKY END IF ! damping defined at element level IF ( FROM_REST ) THEN ! vel = acc = 0 DR = 0.d0 ; D2R = 0.d0 ELSE ! approx initial acc from diagonal scaled mass WORK_1 = A_SKY (I_DIAG) ! extract diagonal of mass matrix D2R = (SUM (A_SKY) / SUM (WORK_1)) * WORK_1 ! scaled M CALL FORCE_AT_TIME (ZERO, P_IN, P) ! initial force CALL VECTOR_MULT_SKY_SYM (B_SKY, DR, I_DIAG, WORK_1) CALL VECTOR_MULT_SKY_SYM (C_SKY, R, I_DIAG, WORK_2) D2R = (P - WORK_1 - WORK_2) / D2R ! approx initial acc END IF ! need initial acc ! Print and/or save initial conditions I_STEP = 0 ; WRITE (6, 5020) I_STEP, ZERO 5020 FORMAT ( /,' STEP NUMBER = ',I5,5X,'TIME = ',E14.8,/, & & ' I R(I) DR/DT ', & & 'D2R/DT2') WRITE (6, 5030) (K, R (K), DR (K), D2R (K), K = 1, N_D_FRE) 5030 FORMAT ( I10, 2X, E14.8, 2X, E14.8, 2X, E14.8 ) write(u_node,'(I3, 4E14.4)') 3, zero, r(3), dr(3), d2r(3) ! Form system work array to factor for each time step group B_SKY = B_SKY + 2 * A_SKY / TAU + TAU * C_SKY / 3 ! ** APPLY BOUNDARY CONDITIONS (ZERO B_SKY ROWS, COLS) ** DO I = 1, N_BC CALL SKY_TYPE_1_SYM (NULL_BC (I), ZERO, B_SKY, P, I_DIAG) END DO ! initial test of bc ! ** TRIANGULARIZE B_SKY, in B_SKY * R = P ** CALL SKY_SOLVE_SYM (B_SKY, P, R, I_DIAG, .TRUE., .FALSE.) ! ** END OF INITIAL CALCULATIONS ** ! ** CALCULATE SOLUTION AT TIME T ** I_COUNT = I_PRINT - 1 DO I_STEP = 1, N_STEPS I_COUNT = I_COUNT + 1 T = DT * I_STEP T_PLUS = T + DELT IF ( I_COUNT == I_PRINT ) WRITE (6, 5020) I_STEP, T ! FORM MODIFIED FORCING FUNCTION AT T + DELT CALL FORCE_AT_TIME (T_PLUS, P_IN, P) CALL VECTOR_MULT_SKY_SYM (A_SKY, (2*DR / TAU + D2R), & I_DIAG, WORK_1) CALL VECTOR_MULT_SKY_SYM (C_SKY, & (R + 2.d0 * TAU * DR / 3 & + TAU * TAU * D2R / 6), & I_DIAG, WORK_2) P = P + WORK_1 - WORK_2 ! ** APPLY NULL BOUNDARY CONDITIONS ( TO P ) ** P ( NULL_BC ) = 0.d0 ! vector subscripts ! SOLVE FOR DR_PLUS AT TIME T+DELT CALL SKY_SOLVE_SYM (B_SKY, P, DR_PLUS, I_DIAG, .FALSE., .TRUE.) ! USING DATA AT T-DT AND T+DELT CALCULATE VALUES AT T WORK_1 = SCALE_1 * DR + DR_PLUS / THETA WORK_2 = SCALE_2 * D2R + 2 * (WORK_1 - DR) / THETA / DT R = R + DT * (2.d0 * DR + WORK_1) / 3 + DT * DT * D2R / 6 DR = WORK_1 ! overwrite previous DR D2R = WORK_2 ! overwrite previous D2R ! OUTPUT RESULTS FOR TIME T write(u_node,'(I3, 4E14.4)') 3, t, r(3), dr(3), d2r(3) IF ( I_COUNT /= I_PRINT ) CYCLE ! to next time step WRITE (6, 5030) (K, R (K), DR (K), D2R (K), K = 1, N_D_FRE) I_COUNT = 0 END DO END SUBROUTINE WILSON_METHOD_SYM SUBROUTINE NEWMARK_METHOD_SYM (I_DIAG, A_SKY, B_SKY, C_SKY, P_IN, R, & DR, DT, N_STEPS, I_PRINT, N_BC, NULL_BC) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! NEWMARK (BETA, GAMMA) CONDITIONALLY STABLE METHOD OF ! STEP BY STEP INTEGRATION OF MATRIX EQUATIONS: ! A*D2R(T)/DT2 + B*DR(T)/DT + C*R(T) = P(T) ! ALSO INCLUDES HUGHES AND TRAPEZOIDAL METHODS ! INITIAL VALUES OF R AND DR ARE PASSED THRU ARGUMENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_COEFF, N_D_FRE IMPLICIT NONE INTEGER, INTENT (IN) :: N_STEPS, I_PRINT, N_BC INTEGER, INTENT (IN) :: I_DIAG (N_D_FRE) INTEGER, INTENT (IN) :: NULL_BC (N_BC) REAL(DP), INTENT (IN) :: DT REAL(DP), INTENT (INOUT) :: A_SKY (N_COEFF), C_SKY (N_COEFF) REAL(DP), INTENT (INOUT) :: B_SKY (N_COEFF) REAL(DP), INTENT (IN) :: P_IN (N_D_FRE) REAL(DP), INTENT (INOUT) :: R (N_D_FRE), DR (N_D_FRE) REAL(DP), PARAMETER :: ZERO = 0.d0 INTEGER :: I, K, I_STEP, I_COUNT REAL(DP) :: T !b REAL(DP) :: WILSON_THETA = 1.25d0, HUGHES_ALPHA = -0.05d0 !b REAL(DP) :: NEWMARK_BETA = 0.25d0, NEWMARK_GAMMA = 0.5d0 !b REAL(DP) :: MASS_DAMPING = 0.d0, STIF_DAMPING = 0.d0 ! average acceleration (1/4, 1/2) unconditional O(dt^2) ! linear acceleration (1/6. 1/2) conditional O(dt^2) ! Fox-Goodwin (1/12, 1/2) conditional O(dt^4) ! Hilber-Hughes-Taylor ([1-alpha]^2 / 4, [1 - 2 alpha]/2), ! -1/3 <= alpha <= 0, unconditional O(dt^2) REAL(DP) :: a1, a2, a3, a4, a5, a6, BETA, GAMMA ! Automatic Arrays REAL(DP) :: R_PLUS (N_D_FRE), D2R (N_D_FRE), P (N_D_FRE) REAL(DP) :: WORK_1 (N_D_FRE), WORK_2 (N_D_FRE) ! N_COEFF = TOTAL NUMBER OF TERMS IN SKYLINE ! N_BC = NUMBER OF D.O.F. WITH SPECIFIED VALUES OF ZERO ! NULL_BC = ARRAY CONTAINING THE N_BC DOF NUMBERS WITH ZERO BC ! R,DR,D2R = 0,1,2 ORDER DERIV. OF R W.R.T. T AT TIME=T ! R_PLUS = VALUE OF R AT TIME = T + DELT ! N_STEPS = NO. OF INTEGRATION STEPS ! I_PRINT = NO. OF INTEGRATION STEPS BETWEEN PRINTING ! ** INITIAL CALCULATIONS ** IF ( NEWMARK_METHOD ) THEN WRITE (6, '("NEWMARK STEP BY STEP INTEGRATION",/)') BETA = NEWMARK_BETA ; GAMMA = NEWMARK_GAMMA ELSEIF ( HUGHES_METHOD ) THEN WRITE (6, '( "HUGHES STEP BY STEP INTEGRATION",/)') BETA = 0.25d0 * (1.d0 - HUGHES_ALPHA) **2 GAMMA = 0.5d0 * (1.d0 - 2 * HUGHES_ALPHA) ELSE ! default to average acceleration method WRITE (6, '( "TRAPEZOIDAL STEP BY STEP INTEGRATION",/)') BETA = 0.25d0 ; GAMMA = 0.5d0 END IF !b print *, 'beta, gamma ', beta, gamma IF ( I_BUG > 0 ) THEN PRINT *,'DEBUG: A, B, C' CALL SKY_FUL_SYM (A_SKY, I_DIAG, P_IN) CALL SKY_FUL_SYM (B_SKY, I_DIAG, P_IN) CALL SKY_FUL_SYM (C_SKY, I_DIAG, P_IN) END IF IF ( N_BC < 1 ) PRINT *,'NOTE: NO CONSTRAINTS IN NEWMARK_METHOD' P = 0.d0 ; WORK_1 = 0.d0 ; WORK_2 = 0.d0 IF ( .NOT. EL_DAMPING ) THEN !b B_SKY = 0.d0 IF ( MASS_DAMPING > 0.d0 ) B_SKY = B_SKY + MASS_DAMPING * A_SKY IF ( STIF_DAMPING > 0.d0 ) B_SKY = B_SKY + STIF_DAMPING * C_SKY END IF ! damping defined at element level IF ( FROM_REST ) THEN ! vel = acc = 0 DR = 0.d0 ; D2R = 0.d0 ELSE ! approx initial acc from diagonal scaled mass WORK_1 = A_SKY (I_DIAG) ! extract diagonal of mass matrix WHERE ( WORK_1 <= 0.d0 ) WORK_1 = 1.d0 END WHERE D2R = (SUM (A_SKY) / SUM (WORK_1)) * WORK_1 ! scaled M CALL FORCE_AT_TIME (ZERO, P_IN, P) ! initial force CALL VECTOR_MULT_SKY_SYM (B_SKY, DR, I_DIAG, WORK_1) CALL VECTOR_MULT_SKY_SYM (C_SKY, R, I_DIAG, WORK_2) D2R = (P - WORK_1 - WORK_2) / D2R ! approx initial acc END IF ! need initial acc ! Print and/or save initial conditions I_STEP = 0 ; WRITE (6, 5020) I_STEP, ZERO 5020 FORMAT ( /,' STEP NUMBER = ',I5,5X,'TIME = ',E14.8,/, & & ' I R(I) DR/DT ', & & 'D2R/DT2') WRITE (6, 5030) (K, R (K), DR (K), D2R (K), K = 1, N_D_FRE) 5030 FORMAT ( I10, 2X, E14.8, 2X, E14.8, 2X, E14.8 ) write(u_node,'(I3, 4E14.4)') 3, zero, r(3), dr(3), d2r(3) ! Form integration constants a1 = 1.d0 / (BETA * DT **2) ; a2 = 1.d0 / (BETA * DT) a3 = 0.5d0/BETA - 1.d0 ; a4 = GAMMA * a2 a5 = GAMMA / BETA - 1.d0 ; a6 = 0.5d0 * GAMMA / BETA - 1.d0 !b print *, 'a1-a6 ', a1, a2, a3, a4, a5, a6 ! Form system work array to factor for each time step group C_SKY = C_SKY + a1 * A_SKY + a4 * B_SKY ! ** APPLY BOUNDARY CONDITIONS (ZERO C_SKY ROWS, COLS) ** DO I = 1, N_BC CALL SKY_TYPE_1_SYM (NULL_BC (I), ZERO, C_SKY, P, I_DIAG) END DO ! initial test of bc ! ** TRIANGULARIZE C_SKY, in C_SKY * R = P ** CALL SKY_SOLVE_SYM (C_SKY, P, R, I_DIAG, .TRUE., .FALSE.) ! ** END OF INITIAL CALCULATIONS ** ! ** CALCULATE SOLUTION AT TIME T ** I_COUNT = I_PRINT - 1 DO I_STEP = 1, N_STEPS I_COUNT = I_COUNT + 1 T = DT * I_STEP IF ( I_COUNT == I_PRINT ) WRITE (6, 5020) I_STEP, T ! FORM MODIFIED FORCING FUNCTION AT T + DELT CALL FORCE_AT_TIME (T, P_IN, P) CALL VECTOR_MULT_SKY_SYM (A_SKY, (a1*R + a2*DR + a3*D2R), & I_DIAG, WORK_1) CALL VECTOR_MULT_SKY_SYM (B_SKY, (a4*R + a5*DR + a6*D2R), & I_DIAG, WORK_2) P = P + WORK_1 + WORK_2 ! ** APPLY NULL BOUNDARY CONDITIONS ( TO P ) ** P ( NULL_BC ) = 0.d0 ! vector subscripts ! SOLVE FOR DR_PLUS AT TIME T+DELT CALL SKY_SOLVE_SYM (C_SKY, P, R_PLUS, I_DIAG, .FALSE., .TRUE.) ! UPDATE KINEMATICS OF D2R(+), DR(+), R(+) WORK_1 = a1*(R_PLUS - R - dt*DR) - a3*D2R ! D2R(+) WORK_2 = a4*(R_PLUS - R) - a5*DR - dt*a6*D2R ! DR(+) R = R_PLUS ! R(+) DR = WORK_2 ! DR(+) D2R = WORK_1 ! D2R(+) ! OUTPUT RESULTS FOR TIME T write(u_node,'(I3, 4E14.4)') 3, t, r(3), dr(3), d2r(3) IF ( I_COUNT /= I_PRINT ) CYCLE ! to next time step WRITE (6, 5030) (K, R (K), DR (K), D2R (K), K = 1, N_D_FRE) I_COUNT = 0 END DO END SUBROUTINE NEWMARK_METHOD_SYM SUBROUTINE FORCE_AT_TIME (T, P_IN, FORCE) ! SYSTEM LEVEL TIME DEPENDENT GENERALIZED FORCE ! FORCE = P_IN * F_OF_T Use System_Constants ! for N_COEFF, N_D_FRE, RAMP_DOWN, RAMP_UP REAL(DP), INTENT (IN) :: T REAL(DP), INTENT (IN) :: P_IN (N_D_FRE) REAL(DP), INTENT (OUT) :: FORCE (N_D_FRE) !b REAL(DP), PARAMETER :: RAMP_DOWN = 0.1d0 ! now global !b REAL(DP), PARAMETER :: RAMP_UP = -1.d0 ! now global ! T = TIME ! P_IN = CONSTANT INPUT GENERALIZED FORCES AT NODES ! FORCE = CURRENT GENERALIZED FORCES AT NODES REAL(DP), SAVE :: F_OF_T F_OF_T = 0.d0 IF ( RAMP_DOWN_TIME > 0.d0 .AND. T < RAMP_DOWN_TIME ) THEN F_OF_T = F_OF_T + 1.d0 - T / RAMP_DOWN_TIME END IF IF ( RAMP_UP_TIME > 0.d0 ) THEN IF ( T < RAMP_UP_TIME ) THEN F_OF_T = F_OF_T + T / RAMP_UP_TIME ELSE F_OF_T = F_OF_T + 1.d0 END IF END IF FORCE = P_IN * F_OF_T END SUBROUTINE FORCE_AT_TIME SUBROUTINE ESTIMATE_SYSTEM_MAX_FREQ (M_SKY, K_SKY, I_DIAG, W_SQ_MAX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ESTIMATE MAXIMUM SYSTEM FREQUENCY, GERSCHGORIN BOUND ! FROM SYMMETRIC MASS AND STIFFNESS MATRICES IN SKYLINE MODE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_COEFF, N_D_FRE Use Interface_Header ! for FIRST_COL_AT_SKY_ROW, LAST_COL_AT_SKY_ROW REAL(DP), INTENT(IN) :: M_SKY (N_COEFF), K_SKY (N_COEFF) INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) REAL(DP), INTENT(OUT) :: W_SQ_MAX ! Square of max frequency REAL(DP) :: W_SQ, M_II, ROW_SUM, K_IJ INTEGER :: ROW, COL_FIRST, COL_LAST, COL W_SQ_MAX = 0.d0 DO ROW = 1, N_D_FRE ! loop over all dof ROW_SUM = 0.d0 M_II = M_SKY ( I_DIAG ( ROW) ) IF ( M_II <= 0.d0 ) THEN PRINT *, 'NOTE: FOUND INVALID MASS AT ROW ', ROW CYCLE ! to next row ELSE ! positive mass COL_FIRST = FIRST_COL_AT_SKY_ROW (ROW, I_DIAG) COL_LAST = LAST_COL_AT_SKY_ROW (ROW, I_DIAG) DO COL = COL_FIRST, COL_LAST ! column loop IN_SKY = LOCATION_IN_SKY (ROW, COL, I_DIAG) ! locate IF ( IN_SKY > 0 ) THEN ! non-zero ? K_IJ = K_SKY (IN_SKY) ! upper value ROW_SUM = ROW_SUM + DABS ( K_IJ ) ! update sum ELSE CYCLE ! to next column ! zero term END IF END DO ! over columns in this row END IF W_SQ = ROW_SUM / M_II ! row freq IF ( W_SQ > W_SQ_MAX ) W_SQ_MAX = W_SQ ! max freq END DO END SUBROUTINE ESTIMATE_SYSTEM_MAX_FREQ