! copyright 2005, J. E. Akin, all rights reserved. PROGRAM DRIVER_F90 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! DRIVER FOR MODEL_F90 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! NOTE: Any subprogram name that does not contain an underscore ! is the F77 version simply converted to F90. That is, it ! has the same argument lists and types as before. Use System_Constants ! for DOF_PT, DOF_PT_SUM, MIN_XYZ etc Use Elem_Type_Data ! with arrays L_TYPE LT_DATA Use Sys_Properties_Data Use Interface_Header Use Select_Source IMPLICIT NONE REAL(DP), ALLOCATABLE :: X (:,:) ! nodal coordinates REAL(DP), ALLOCATABLE :: CC_IN (:) ! initial external sources REAL(DP), ALLOCATABLE :: MASSES (:) ! initial point masses REAL(DP), ALLOCATABLE :: COEF_C_EQ (:,:) ! constraint coefficients REAL(DP), ALLOCATABLE :: DAMPENER (:) ! initial point dampners !b REAL(DP), ALLOCATABLE :: DD_DOT (:) ! system velocities !b REAL(DP), ALLOCATABLE :: DD_2DOT (:) ! system accelerations INTEGER, ALLOCATABLE :: I_BC (:) ! nodal constraint flags INTEGER, ALLOCATABLE :: NODES (:,:) ! mesh connectivity INTEGER, ALLOCATABLE :: N_RES_EQ (:) ! constraint type counter INTEGER, ALLOCATABLE :: NDEX_C_EQ (:, :)! constraint dof numbers ! below now in system_constants !b INTEGER, ALLOCATABLE :: DOF_PT (:) ! number of dof at each node !b INTEGER, ALLOCATABLE :: DOF_PT_SUM (:) ! total number of equations INTEGER :: IO_TEST, J ! status of test.dat, loop INTEGER :: CASE_EXACT ! dummy check CHARACTER (len = 64) :: FileName ! input data filename CHARACTER (LEN=16) :: STYLE ! style of control input LOGICAL :: PRT_OPEN = .FALSE. ! ** SET DEFAULTS, OPEN INPUT AND LOG FILES ** CALL SET_CONTROL_DEFAULTS ! including N_CRD, N_PRT ! Get input file name !b WRITE (*, *) 'Enter input data file name (default = test.dat):' !b READ (*, '(A64)') FileName !b IF ( FileName == '' ) THEN ! use default name !b FileName = 'test.dat' !b END IF !b ! Put controls in model.ctl INQUIRE (N_PRT, OPENED = PRT_OPEN) !b IF ( PRT_OPEN ) CLOSE (N_PRT) !b OPEN (N_PRT, FILE = 'model.ctl', ACTION = 'WRITE', & FORM = 'FORMATTED', STATUS ='REPLACE', IOSTAT = IO_TEST) !b !b OPEN (N_CRD, file='test.dat', status='old', action='read', & OPEN (N_CRD, file=FileName, status='old', action='read', & iostat = IO_TEST) IF ( IO_TEST > 0 ) THEN !b PRINT *,'ERROR: Required data file test.dat failed to open' PRINT *,'ERROR: Required data ', FileName, ' failed to open' PRINT *, 'Copy your data to test.dat and retry' STOP 'ERROR: Copy your data to test.dat and retry' END IF U_HIST = 0 ; U_LOG = 0 ! no longer needed !b !b IF ( U_HIST > 0 ) OPEN (U_HIST, file='model.his', status='unknown', & !b action='write', position='append') !b IF ( U_LOG > 0 ) OPEN (U_LOG , file='model.log', status='unknown', & !b action='write') ! ** AUTHOR CREDIT ** PRINT *, 'MODular Element Library (MODEL) 4.3.0' ! June 1, 2004 PRINT *, 'J.E. Akin, Rice University, akin@rice.edu' ! ** READ APPLICATION CONTROL DATA ** !b READ (N_CRD, *) STYLE ! style of user control !b SELECT CASE (STYLE) !b CASE ('keywords') ; CALL PARSE_INPUT !b CASE ('title') ; BACKSPACE (N_CRD) !b CALL PARSE_INPUT !b CASE DEFAULT !b PRINT *, "NO keywords DEFAULTING TO old_style INPUT MODE" !b BACKSPACE (N_CRD) !b CALL FIXED_CONTROL_INPUT !b END SELECT CALL PARSE_INPUT !b ! Option to suppress printing, except to *.tmp files INQUIRE (N_PRT, OPENED = PRT_OPEN) IF ( PRT_OPEN ) CLOSE (N_PRT) ! switch from control to general output IF ( NO_PRINTING ) THEN ! via control no_printing OR omit_prints OPEN (N_PRT, STATUS = 'SCRATCH') ! no print output from here on ELSE ! send to a file instead of the screen OPEN (N_PRT, FILE = 'model.out', ACTION = 'WRITE', Form = 'FORMATTED', & STATUS ='REPLACE', IOSTAT = IO_TEST) print *,'NO_PRINTING ', NO_PRINTING END IF IF ( EXAMPLE > 0 ) CALL SELECT_DESCRIBE_EXAMPLE IF ( EXACT_CASE > 0 ) PRINT *,'NOTE: WILL USE EXACT CASE NUMBER ', EXACT_CASE CALL VALIDATE_AND_LIST_CONTROLS IF ( USE_EXACT .OR. USE_EXACT_BC .OR. USE_EXACT_FLUX .OR. & EXACT_CASE > 0 ) CALL SELECT_EXACT_CASE_NUMBER (CASE_EXACT) !b CALL ALLOCATE_SYSTEM_ARRAYS ! currently body force vector CALL ALLOCATE_OUTPUT_NAMES ALLOCATE ( I_BC (MAX_NP) ) ; I_BC = 0 ALLOCATE ( NODES (L_S_TOT, NOD_PER_EL) ) ; NODES = 0 ALLOCATE ( N_RES_EQ (MAX_TYP) ) ALLOCATE ( X (MAX_NP, N_SPACE) ) ALLOCATE ( MAX_XYZ (N_SPACE) ) ! for system constants ALLOCATE ( MIN_XYZ (N_SPACE) ) ! for system constants if ( .not. ALLOCATED (NODES )) stop ' NODES not allocated' if ( .not. ALLOCATED (N_RES_EQ )) stop ' N_RES_EQ not allocated' if ( .not. ALLOCATED (I_BC )) stop ' I_BC not allocated' if ( .not. ALLOCATED (X )) stop ' X not allocated' IF ( I_BUG > 0 ) PRINT *, 'ALLOCATED MESH DATA' ALLOCATE ( L_TYPE (L_S_TOT), LT_DATA (MAX_DAT, N_L_TYPE) ) TYPE_DATA_ALLOC = .TRUE. IF ( I_BUG > 0 ) PRINT *, 'ALLOCATED TYPE DATA' IF ( DYNAMIC ) THEN ! Allocate nodal velocities and accelerations ALLOCATE ( DD_DOT (N_D_FRE) ) ALLOCATE ( DD_2DOT (N_D_FRE) ) IF ( I_BUG > 0 ) PRINT *, 'ALLOCATED VELOCITY & ACCELERATION' ELSE ! debugging static iter ALLOCATE ( DD_DOT (1) ) ALLOCATE ( DD_2DOT (1) ) END IF DD_DOT = 0.d0 ; DD_2DOT = 0.d0 ! initialize !--> *** SET UP, MAYBE READ PROPERTIES *** IF ( IP_TEST > 0 ) THEN CALL ALLOCATE_PROPERTIES IF ( EXACT_INT > 0 .OR. EXACT_REALS > 0 ) & CALL READ_EXACT_SOLUTION_PROPERTIES IF ( PROPERTY_FIRST ) CALL INPUT_PROPERTIES NO_PROPERTIES = .FALSE. ELSE WRITE (6, * ) ' ' WRITE (6, * ) 'WARNING, NO PROPERTY INPUT' N_WARN = N_WARN + 1 ! INCREMENT WARNING NO_PROPERTIES = .TRUE. ALLOCATE ( DELETED ( L_S_TOT)) ; DELETED = .FALSE. PROP_ELEM_ALLOC = .TRUE. END IF ! ALLOW FOR VARIABLE NUMBER OF EQUATIONS PER NODE IF ( DOF_VARY ) THEN ALLOCATE ( DOF_PT (0:MAX_NP) ) ALLOCATE ( DOF_PT_SUM (-1:MAX_NP) ) DOF_PT (0) = 0 DOF_PT (1:MAX_NP) = N_G_DOF ! default to constant DOF_PT_SUM (-1:0) = 0 DO J = 1, MAX_NP DOF_PT_SUM (J) = DOF_PT_SUM (J-1) + DOF_PT (J) END DO !b print *,'DOF_PT ', DOF_PT (0:MAX_NP) !b print *,'DOF_PT_SUM', DOF_PT_SUM (0:MAX_NP) END IF ! INPUT NODES, BC FLAGS, ELEMENTS (AND TYPES) CALL INPUT_MESH (X, I_BC, NODES, N_RES_EQ) MAX_XYZ = MAXVAL (X, DIM = 1) ! Bounding coordinates MIN_XYZ = MINVAL (X, DIM = 1) ! Bounding coordinates IF ( I_BUG > 0 ) THEN PRINT *, '*** BOUNDING COORDINATE VALUES ***' PRINT *, 'SIZES ', SIZE (MIN_XYZ), SIZE (MAX_XYZ) PRINT *, 'MAX: ', MAX_XYZ PRINT *, 'MIN: ', MIN_XYZ END IF ALLOCATE ( COEF_C_EQ (MAX_ACT, N_CEQ) ) ! FROM I_BC DATA ALLOCATE ( NDEX_C_EQ (MAX_ACT, N_CEQ) ) ! FROM I_BC DATA !--> READ MIXED_BC SEGMENT NODES IF ( N_MIXED > 0 ) THEN CALL INPUT_MIXED_BC_SEG (L_TYPE, NODES) !b CALL SAVE_MIXED_BC_SEG (L_TYPE, NODES) !b to be added END IF !--> READ FLUX SEGMENT NODES AND FLUX COMPONENTS IF ( N_F_SEG > 0 ) THEN IF ( USE_EXACT_FLUX ) THEN CALL INPUT_EXACT_FLUX_SEG (L_TYPE, NODES, X) ELSE CALL INPUT_FLUX_SEG (L_TYPE, NODES) END IF !b CALL SAVE_FLUX_SEG (L_TYPE, NODES) !b to be added END IF ! SAVE ELEMENT AND SEGMENT TYPES AND CONNECTIVITY CALL SAVE_MESH (X, I_BC, L_TYPE, NODES) !b IF ( I_BUG > 0 ) THEN ! LIST THE ONLY ARRAYS KNOWN AT THIS POINT PRINT *, 'X '; CALL RPRINT (X , 1 , MAX_NP ,1) PRINT *, 'I_BC '; CALL IPRINT (I_BC , 1 , MAX_NP) PRINT *, 'L_TYPE '; CALL IPRINT (L_TYPE , 1 , L_S_TOT) PRINT *, 'NODES '; CALL IPRINT (NODES , L_S_TOT, NOD_PER_EL) PRINT *, 'N_RES_EQ '; CALL IPRINT (N_RES_EQ, 1 , MAX_TYP) END IF ! OPEN SEQUENTIAL FILES IF ( N_REACT > 0) OPEN (N_REACT, STATUS = 'SCRATCH', & ACTION = 'READWRITE', FORM = 'UNFORMATTED') !b FILE = 'N_REACT_TMP.BIN', & IF ( L_REACT > 0) OPEN (L_REACT, STATUS = 'SCRATCH', & ACTION = 'READWRITE', FORM = 'UNFORMATTED') !b FILE = 'L_REACT_TMP.BIN', & IF ( U_FLUX > 0) OPEN (U_FLUX, FILE = 'FLUX_SAVE.BIN ', & STATUS = 'REPLACE', & ACTION = 'READWRITE', FORM = 'UNFORMATTED') IF ( U_INTGR > 0) OPEN (U_INTGR, FILE = 'INTEGRAL_SAVE.BIN ', & STATUS = 'REPLACE', & ACTION = 'READWRITE', FORM = 'UNFORMATTED') IF ( N_FILE1 > 0) OPEN (N_FILE1, FILE = 'N_FILE1.BIN', & ACTION = 'READWRITE', STATUS = 'REPLACE', & FORM = 'UNFORMATTED') IF ( N_FILE2 > 0) OPEN (N_FILE2, FILE = 'N_FILE2.BIN', & ACTION = 'READWRITE', STATUS = 'REPLACE', & FORM = 'UNFORMATTED') IF ( N_FILE3 > 0) OPEN (N_FILE3, FILE = 'N_FILE3.DAT', &!b ACTION = 'READWRITE', STATUS = 'REPLACE', &!b FORM = 'FORMATTED')!b IF ( N_FILE4 > 0) OPEN (N_FILE4, FILE = 'N_FILE4.DAT', &!b ACTION = 'READWRITE', STATUS = 'REPLACE', &!b FORM = 'FORMATTED')!b IF ( N_FILE5 > 0) OPEN (N_FILE5, FILE = 'N_FILE5.DAT', &!b ACTION = 'READWRITE', STATUS = 'REPLACE', &!b FORM = 'FORMATTED')!b !--> *** READ NODAL PARAMETER CONSTRAINT EQUATIONS *** IF ( I_BUG > 0 ) THEN ! Print new array PRINT *, 'LT_DATA' ; CALL IPRINT (LT_DATA, MAX_DAT, N_L_TYPE) END IF IF ( USE_EXACT_BC ) THEN CALL INPUT_EXACT_BC_AND_MPC (X, N_RES_EQ, COEF_C_EQ, NDEX_C_EQ) ELSE CALL INPUT_CONSTRAINT_EQS (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ) END IF ! EXACT ESSENTIAL BC !--> *** READ PROPERTIES (DEFAULT CASE) *** IF ( .NOT. PROPERTY_FIRST ) CALL INPUT_PROPERTIES !--> ** INPUT INITIAL FORCING VECTOR ** IF ( IN_RHS > 0 ) THEN IN_RHS = N_D_FRE ALLOCATE ( CC_IN (IN_RHS) ) CALL INPUT_RHS (CC_IN) END IF !--> ** INPUT INITIAL POINT MASSES ** IF ( PT_MASSES ) THEN IN_RHS = N_D_FRE ALLOCATE ( MASSES (N_D_FRE) ) CALL INPUT_MASSES (MASSES) ELSE IF ( JACOBI .OR. DYNAMIC ) PRINT *,'NOTE: NO POINT MASSES INPUT.' ALLOCATE ( MASSES (1) ) END IF !--> ** INPUT INITIAL POINT DAMPENER ** IF ( PT_DAMPING ) THEN IN_RHS = N_D_FRE ALLOCATE ( DAMPENER (N_D_FRE) ) CALL INPUT_DAMPENER (DAMPENER) ELSE IF ( DYNAMIC ) PRINT *,'NOTE: NO POINT DAMPENER INPUT.' ALLOCATE ( DAMPENER (1) ) END IF !--> ** INPUT INITIAL DELETED ELEMENT LIST IF ( GET_DELETED_EL ) CALL INPUT_DELETED_ELEMENTS IF ( DEBUG_UNITS ) THEN CALL LIST_SYSTEM_UNIT_NUMBERS CALL LIST_OPEN_IO_UNITS END IF ! -------------------------- ! --- Call the true main program --- ! -------------------------- !b IF ( EIGEN ) CALL EIGEN_MODEL (X, I_BC, NODES, N_RES_EQ, & !b COEF_C_EQ, NDEX_C_EQ, CC_IN) IF ( JACOBI ) CALL JACOBI_MODEL (X, I_BC, NODES, N_RES_EQ, & COEF_C_EQ, NDEX_C_EQ, MASSES) !b COEF_C_EQ, NDEX_C_EQ, CC_IN) IF ( STATIC ) CALL STATIC_MODEL (X, I_BC, NODES, N_RES_EQ, & COEF_C_EQ, NDEX_C_EQ, CC_IN) IF ( TRANSIENT ) CALL TRANSIENT_MODEL (X, I_BC, NODES, N_RES_EQ, & COEF_C_EQ, NDEX_C_EQ, CC_IN) IF ( DYNAMIC ) CALL DYNAMIC_MODEL (X, I_BC, NODES, N_RES_EQ, & COEF_C_EQ, NDEX_C_EQ, & MASSES, DAMPENER, CC_IN) IF ( DYNAMIC ) THEN DEALLOCATE ( DD_2DOT ) !b not really used DEALLOCATE ( DD_DOT ) END IF DEALLOCATE ( LT_DATA, L_TYPE ) TYPE_DATA_ALLOC = .FALSE. IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED TYPE DATA' IF ( ALLOCATED (NDEX_C_EQ)) DEALLOCATE ( NDEX_C_EQ ) IF ( ALLOCATED (COEF_C_EQ)) DEALLOCATE ( COEF_C_EQ ) IF ( ALLOCATED (MIN_XYZ )) DEALLOCATE ( MIN_XYZ ) !system constants IF ( ALLOCATED (MAX_XYZ )) DEALLOCATE ( MAX_XYZ ) !system constants IF ( ALLOCATED (X )) DEALLOCATE ( X ) IF ( ALLOCATED (N_RES_EQ )) DEALLOCATE ( N_RES_EQ ) IF ( ALLOCATED (I_BC )) DEALLOCATE ( I_BC ) IF ( ALLOCATED (NODES )) DEALLOCATE ( NODES ) IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED MESH DATA' IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION IF ( TYPE_NTRP_ALLOC ) CALL DEALLOCATE_TYPE_INTERPOLATIONS IF ( I_BUG > 0 ) CALL LIST_TYPE_ALLOC_STATUS IF ( N_WARN == 0 ) THEN PRINT *, 'NORMAL END OF MODEL_F90 (NO WARNINGS)' !b STOP 'NORMAL END OF MODEL_F90 (NO WARNINGS)' ELSE PRINT *, 'NORMAL END OF MODEL_F90, WITH ', N_WARN,' WARNINGS' !b STOP 'NORMAL END OF MODEL_F90, WITH WARNINGS' END IF END PROGRAM DRIVER_F90 SUBROUTINE STATIC_MODEL (X, I_BC, NODES, N_RES_EQ, COEF_C_EQ, & NDEX_C_EQ, CC_IN) ! ****************************************************** ! * MODULAR PROGRAMS FOR FINITE ELEMENT ANALYSES * ! * VERSION 1, STATIC ITERATIVE OR INCREMETAL SOLUTION * ! * COPYRIGHT J. E. AKIN, 1998 * ! ****************************************************** ! 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 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 ! for ON_BOUNDARY IMPLICIT NONE LOGICAL :: FACT, BACK ! Given Arrays REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP), INTENT(IN) :: COEF_C_EQ (MAX_ACT, N_CEQ) REAL(DP), INTENT(IN) :: CC_IN (IN_RHS) 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 debug itter dynamic !b REAL(DP), INTENT(IN) :: DD_DOT (N_D_FRE), DD_2DOT (N_D_FRE) !b ! Allocatable Arrays REAL(DP), ALLOCATABLE :: M_SKY (:) ! symmetric mass REAL(DP), ALLOCATABLE :: M_SKY_LOWER (:) ! lower unsymmetric mass REAL(DP), ALLOCATABLE :: S_SKY (:) ! symmetric part REAL(DP), ALLOCATABLE :: S_SKY_LOWER (:) ! unsymmetric part INTEGER, ALLOCATABLE :: I_DIAG (:), I_DOF_HI (:) ! sky parts INTEGER, ALLOCATABLE :: L_TO_L_NEIGH (:, :) ! neighbors INTEGER, ALLOCATABLE :: L_TO_N_NEIGH (:, :) ! neighbors ! Automatic Arrays REAL(DP) :: CC (N_D_FRE) REAL(DP) :: DD (N_D_FRE), DD_OLD (N_D_FRE) REAL(DP) :: SCP_AVERAGES (MAX_NP, SCP_FIT) REAL(DP) :: EL_ENORM (N_ELEMS), ERR_MAX REAL(DP) :: EL_REFINE (N_ELEMS) INTEGER :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER :: L_TO_L_SUM (N_ELEMS) INTEGER :: L_TO_N_SUM (MAX_NP) INTEGER :: ITER, IBAR, IPARM, IE, OLD_DEL = -1 INTEGER :: COL_HI_SQ INTEGER :: IO_1, I_WARN REAL(DP) :: R_TEST, RATIO, OLD_NORM, NEW_NORM, DIFF_NORM ! added to static_model in driver_lib.f INTEGER :: IO_R, U_TMP, GET_NEXT_IO_UNIT ! critical speed REAL(DP) :: W_Y, W_Y2, R_CS, R_CM, R_CF, R_CI ! Rayleigh estimate LOGICAL :: CRITICAL_RPM = .TRUE. ! SEE NOTATION.F FOR VARIABLE AND ARRAY NOTATIONS IF ( I_BUG > 0 ) PRINT *,'=== USING STATIC_MODEL ===' CALL ALLOCATE_GEOM_PROP ! includes ON_BOUNDARY !b call at(296) !b print *,'ALLOCATED(ON_BOUNDARY)', ALLOCATED(ON_BOUNDARY) !b print *,'SIZE(ON_BOUNDARY)', SIZE(ON_BOUNDARY) ! GENERATE & LIST MESH NEIGHBORS (PATCH) INFORMATION IF ( SCP_NEIGH_FACE ) CALL SET_L_SHAPE_FACE_NEEDS ! for NEEDS I_WARN = 1 ! dummy CALL FIRST_LAST_L_AT_PT (NODES, L_FIRST, L_LAST, I_WARN) 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) 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 ( PT_EL_SUMS .OR. PT_EL_LIST ) & CALL LIST_L_AT_NODE_INFO (L_TO_N_SUM, L_FIRST, L_LAST) IF ( PT_EL_LIST .OR. SCP_NEIGH_PT ) 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) END IF IF ( I_BUG > 0 ) THEN IF ( NEIGH_L > 0 ) PRINT *, 'L_TO_L Sparseness Value ', & float( count(L_TO_L_NEIGH == 0)) / (NEIGH_L * N_ELEMS) IF ( NEIGH_N > 0 ) PRINT *, 'L_TO_N Sparseness Value ', & float( count(L_TO_N_NEIGH == 0)) / (NEIGH_N * MAX_NP) END IF ! TO DEBUG !!b force proper on_boundary codes !if ( NEEDS == 0 ) then ! OPEN (67, FILE='on_boundary_input', & !b ! ACTION='READ', STATUS='OLD', IOSTAT = IO_1) !b ! IF ( IO_1 > 0 ) THEN ! PRINT *,'OPENED FAILED FOR on_boundary_input' ! ELSE ! ON_BOUNDARY = .FALSE. ! DO ! FOREVER ! READ (67, *, IOSTAT = IO_1) IE ! IF ( IO_1 /= 0 ) EXIT ! THIS LOOP ! ON_BOUNDARY (IE) = .TRUE. ! END DO ! FOR FACE FLAGS ! END IF ! boundary file provided !end if ! NEEDS !!b !--> ** INPUT INITIAL DELETED ELEMENT LIST IF ( GET_DELETED_EL ) CALL INPUT_DELETED_ELEMENTS ! INITIALIZE SYSTEM DOF FOR ITERATIVE SOLUTION DD = 0.d0 ; DD_OLD = 0.d0 !b IF ( N_ITER > 1) THEN IF ( RESTART_STEP > 0 ) THEN THIS_ITER = RESTART_STEP CALL GET_STEP_RESULTS (DD_OLD) ELSE THIS_ITER = 0 CALL START_ITERATION (X, DD_OLD) END IF !b ELSE !b DD_OLD = 0.d0 END IF !b IF ( I_BUG > 0 ) THEN !b print *,'No iteration start debug active' !b END IF 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 ! *** BEGIN ITERATION LOOP *** DO ITER = 1, N_ITER THIS_ITER = RESTART_STEP + ITER ! Set global value IF ( N_ITER > 1 ) THEN PRINT *, '--------------------------------------------------' PRINT *, ' BEGINNING ITERATION NUMBER ', ITER PRINT *, '--------------------------------------------------' END IF !b call at(463) ; !b print *, 'DD_OLD (1:N_EL_FRE) ',DD_OLD (1:N_EL_FRE) !b print *, 'DD (1:N_EL_FRE) ',DD (1:N_EL_FRE) IF ( ITER > 1 ) THEN IF ( N_REACT > 0 ) REWIND (N_REACT) IF ( L_REACT > 0 ) REWIND (L_REACT) IF ( U_FLUX > 0 ) REWIND (U_FLUX ) IF ( N_FILE1 > 0 ) REWIND (N_FILE1) IF ( N_FILE2 > 0 ) REWIND (N_FILE2) IF ( N_FILE3 > 0 ) REWIND (N_FILE3) IF ( N_FILE4 > 0 ) REWIND (N_FILE4) END IF N_L_DEL = COUNT (DELETED) IF ( N_L_DEL > 0 ) THEN print *,' ' print *,'NUMBER OF ELEMENTS OMITTED = ', N_L_DEL print *,'** OMITTED ELEMENT LIST ** ' DO IE = 1, L_S_TOT IF ( DELETED (IE) ) print *, IE END DO END IF IF ( N_L_DEL > OLD_DEL ) THEN ! more elements omitted this iter OLD_DEL = N_L_DEL ! 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) COL_HI_SQ = SUM ( I_DOF_HI **2 ) 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) IF ( N_COEFF < N_D_FRE ) STOP 'STORAGE ERROR IN STATIC_MODEL' 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, N_STORE, COL_HI_SQ !b 5 FORMAT (/, 'TOTAL NUMBER OF SYSTEM EQUATIONS ....',I10,& !b /, 'NUMBER OF STIFFNESS COEFFICIENTS ....',I10,& !b /, 'SUM OF HEIGHTS SQUARED ..............',I10) END IF ! number of active elements reduced ! ** 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 OF 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 = 0 ! NO MASS ALLOCATE ( M_SKY (M_COEFF) ) ! add status IF ( .NOT. SYMMETRIC ) ALLOCATE ( M_SKY_LOWER (M_COEFF) ) ! status ALLOCATE ( S_SKY (N_COEFF) ) ! add status S_SKY = 0.0_DP ! zero stiffness IF ( .NOT. SYMMETRIC ) THEN ALLOCATE ( S_SKY_LOWER(NOT_SYM) ) ! add status S_SKY_LOWER = 0.0_DP ! zero stiffness END IF IF ( I_BUG > 0 ) PRINT *, ' ALLOCATED SKYLINE REAL ARRAYS' ! INITIALIZE THE SYSTEM SOURCE MATRIX DD = 0.0_DP CC = 0.0_DP ! zero load IF ( IN_RHS > 0 ) CC = CC_IN ! initialize load IF ( I_BUG > 0 ) CALL SUM_RHS (CC) !--> *** CALCULATE AND ASSEMBLE ELEMENT MATRICES *** !--> *** GENERATE POST SOLUTION MATRICES & STORE *** !b print *,'ITER 1, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) CALL ASSEMBLE_SKYLINE_EQS (I_DIAG, NODES, S_SKY, CC, X, & DD_OLD, S_SKY_LOWER, ITER, M_SKY, M_SKY_LOWER) !b print *,'ITER 2, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) IF ( I_BUG > 0 ) CALL SUM_RHS (CC) ! *** ASSEMBLY COMPLETED, CHECK SOURCES *** IF ( I_BUG > 0 .AND. N_D_FRE < M_SHOW) THEN PRINT *,'DEBUG PRINT AFTER ASSEMBLY IN FULL MATRIX FORM:' CALL SKY_FUL_LOCATIONS (I_DIAG) ! show locations IF ( I_BUG == 1 ) THEN ! show values CALL SKY_FUL_GEN (S_SKY, I_DIAG, CC, S_SKY_LOWER) END IF ! SHOW VALUES END IF CALL SUM_RHS (CC) ! ** SAVE DATA FOR REACTION RECOVERY ** IF ( SUM (I_BC) > 0 ) THEN ! reactions exist IF ( SYMMETRIC ) THEN CALL SAVE_SKY_ROWS_SYM (I_BC, S_SKY, CC, I_DIAG) ELSE CALL SAVE_SKY_ROWS_GEN (I_BC, S_SKY, CC, I_DIAG, S_SKY_LOWER) END IF END IF ! for reactions ! CRITICAL SPEED ESTIMATE DESIRED ? !b IF ( CRITICAL_RPM ) THEN ! have to save the force vector !b U_TMP = GET_NEXT_IO_UNIT () !b OPEN (U_TMP, FILE='critical_rpm.bin', ACTION='READWRITE', & !b FORM = 'UNFORMATTED', STATUS='REPLACE', IOSTAT = IO_R) !b IF ( IO_R > 0 ) THEN !b PRINT *, 'WARNING: OPEN FAILED, critical_rpm.bin' !b N_WARN = N_WARN + 1 ; RETURN !b ELSE ! save CC that is destroyed by BC & solver !b WRITE (U_TMP) CC ! binary write !b END IF ! OPEN FAILURE !b END IF ! estimate critical speed !b !--> ** APPLY BOUNDARY CONSTRAINTS TO NODAL PARAMETERS ** !b IF ( I_BUG > 0 ) PRINT *, 'Before SKY_BC_GEN' CALL SKY_BC_GEN (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, S_SKY, CC, & I_DIAG, S_SKY_LOWER) !b print *,'ITER 3, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) IF ( I_BUG == 1 .AND. N_D_FRE < M_SHOW ) THEN PRINT *,'DEBUG PRINT AFTER BC APPLICATION' CALL SKY_FUL_GEN (S_SKY, I_DIAG, CC, S_SKY_LOWER) END IF ! ** CHECK OR FIX SQUARE MATRIX ** ! IF ( SYMMETRIC ) THEN ! assume positive definite CALL CHECK_SKY_SYS_SYM (S_SKY, CC, I_DIAG) ! else add checks for general system END IF !--> *** SOLVE FOR UNKNOWN NODAL PARAMETERS *** ! FACT = .TRUE. ; BACK = .TRUE. IF ( SYMMETRIC ) THEN ! assume positive definite CALL SKY_SOLVE_SYM (S_SKY, CC, DD, I_DIAG, FACT, BACK) ELSE CALL SKY_SOLVE_UNSYM (S_SKY, S_SKY_LOWER, CC, DD, I_DIAG, & FACT, BACK) END IF !b print *,'ITER 4, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) !--> *** SOLUTION COMPLETE, GET SYSTEM REACTIONS *** ! IF ( I_BUG > 0 ) THEN print *,'No solution check debug active' END IF IF ( N_REACT > 0 .AND. SUM(I_BC) > 0) CALL GET_REACTIONS (DD) ! *** PRINT RESULTS *** CALL MAX_AND_MIN_RANGE (DD) IF ( NPT_WRT == 1) THEN IF ( USE_EXACT ) THEN CALL WRITE_BY_NODES_WITH_EXACT (X, DD) ELSE CALL WRITE_BY_NODES (X, DD) END IF ! EXACT SOLUTION PROVIDED (IN EXACT_SOLUTION) END IF ! LIST BY NODES IF ( LEM_WRT == 1) CALL WRITE_BY_ELEMENTS (DD, NODES) IF ( U_PLT3 > 0 ) CALL SAVE_RESULTS (DD) IF ( U_PLT8 > 0 .AND. USE_EXACT ) CALL & SAVE_EXACT_SOLUTION_AT_NODES (X) IF ( U_PLT9 > 0 .AND. USE_EXACT ) CALL & SAVE_EXACT_FLUXES_AT_NODES (X) ! CRITICAL SPEED ESTIMATE DESIRED ? !b IF ( CRITICAL_RPM ) THEN !b REWIND (U_TMP) ; READ (U_TMP) CC ! binary read !b CLOSE (U_TMP) !b OPEN (U_TMP, FILE='critical_rpm.tmp', ACTION='WRITE', & !b STATUS='REPLACE', IOSTAT = IO_R) !b IF ( IO_R > 0 ) THEN !b PRINT *, 'WARNING: OPEN FAILED, critical_rpm.tmp' !b N_WARN = N_WARN + 1 ; RETURN !b ELSE ! save results, assume null BC on DD !b W_Y = DOT_PRODUCT ( ABS (CC), ABS (DD) ) !b W_Y2 = DOT_PRODUCT ( ABS (CC), (DD **2) ) + TINY (1.d0) !b R_CS = SQRT (W_Y / W_Y2) * 9.549 !b R_CM = R_CS * 3.13 ; R_CF = R_CS * 5.67 ; R_CI = R_CS * 19.6 !b PRINT *, ' ' !b PRINT *, 'RAYLEIGH CRITICAL SHAFT SPEED ESTIMATE' !b PRINT *, 'RPM = sqrt (G) * ', R_CS !b PRINT *, 'FOR G = 9.81 m/s^2, RPM = ', R_CM !b PRINT *, 'FOR G = 32.2 f/s^2, RPM = ', R_CF !b PRINT *, 'FOR G = 386.4 in/s^2, RPM = ', R_CI !b WRITE (U_TMP, '( 4(1PE11.3) )') R_CS, R_CM, R_CF, R_CI !b PRINT *, 'NOTE: RAYLEIGH RESULT SAVED TO critical_rpm.tmp' !b CLOSE (U_TMP) !b END IF ! OPEN FAILURE !b END IF ! estimate critical speed !b ! RESET THE REACTION DATA AND SQ MATRIX FOR NEXT ITERATION IF ( N_REACT > 0) REWIND (N_REACT) IF ( .NOT. SYMMETRIC ) DEALLOCATE ( S_SKY_LOWER ) DEALLOCATE ( S_SKY ) IF ( .NOT. SYMMETRIC ) DEALLOCATE ( M_SKY_LOWER ) !b DEALLOCATE ( M_SKY ) !b IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED SKYLINE REAL ARRAYS' ! ** OPTIONAL ELEMENT REACTION RECOVERY *** !b IF ( L_REACT > 0 ) CALL LIST_REACTIONS_AT_ELEMS (NODES, DD, CC_IN) IF ( L_REACT > 0 ) CALL LIST_REACTIONS_AT_ELEMS (NODES, DD) !b !--> *** POST SOLUTION CALCULATIONS *** IF ( U_FLUX > 0 .OR. N_FILE1 > 0 .OR. N_FILE2 > 0 ) THEN !sc IF ( U_SCPR > 0 .AND. N_QP > 0 ) THEN ! POSSIBLE SCP RECOVERY, PREPARE RECORDS SCP_PAD = EST_SCP_PAD_SIZE () SCP_RECL = SET_SCP_RECORD_SIZE () CALL INITIALIZE_SCP_RECORDS ! FOR FILLING ELSE ! REMIND USER PRINT *, ' ' PRINT *, 'NOTE: SUPER_CONVERGENT PATCH RECOVERY NOT ACTIVE' PRINT *, ' (UNIT U_SCPR = 0 OR N_QP = 0)' PRINT *, ' ' END IF ! POSSIBLE SCP RECOVERY ! POST-PROCESS GRADIENTS & GENERATE SCP DATA FOR AVERAGING IF ( U_FLUX > 0 ) CALL POST_PROCESS_GRADS (NODES, DD, ITER) IF ( SCP_RECORDS_SAVED ) THEN ! *** GET SCP NODAL FLUX FROM FLUX (DEFAULT) OR DOF *** IF ( SCP_NEIGH_PT ) THEN ! NODE BASED PATCH NEIGHBORS PRINT *, 'NOTE: SCP BASED ON NODE CENTER PATCHES' N_PATCH = MAX_NP IF ( SCP_DOF ) THEN ! ave nodal dof & get gradients !b CALL SURF_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_N_NEIGH, & DD, SCP_AVERAGES) ELSE ! standard scp CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_N_NEIGH, & SCP_AVERAGES) END IF ELSE IF ( SCP_NEIGH_FACE ) THEN ! FACE BASED PATCHES PRINT *, 'NOTE: SCP BASED ON FACING ELEMENT PATCHES' N_PATCH = N_ELEMS ! Only face neighbors in L_TO_L_NEIGH IF ( SCP_DOF ) THEN ! ave nodal dof & get gradients !b CALL SURF_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH, & DD, SCP_AVERAGES) ELSE ! standard scp CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH, & SCP_AVERAGES) END IF ELSE ! ELEMENT BASED NEIGHBORS (DEFAULT) PRINT *, 'NOTE: SCP BASED ON ADJACENT ELEMENT PATCHES' N_PATCH = N_ELEMS ! All element neighbors in L_TO_L_NEIGH IF ( SCP_DOF ) THEN ! ave nodal dof & get gradients !b CALL SURF_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH, & DD, SCP_AVERAGES) ELSE ! standard scp CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH, & SCP_AVERAGES) END IF END IF ! BASIS FOR PATCH ELEMENTS !b ! ** ARE FLUXES FROM FLUX (DEFAULT) OR DOF AVERAGES !b !b IF ( SCP_DOF ) THEN ! ave nodal dof & get gradients !b !b CALL SURF_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH, & !b !b DD, SCP_AVERAGES) !b !b ELSE ! ave flux from quadrature points !b !b CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH, & !b !b SCP_AVERAGES) !b !b END IF ! averaging type desired !b IF ( USE_EXACT_FLUX ) THEN IF ( GRAD_BASE_ERROR ) PRINT *, & 'NOTE: FLUX IS ACTUALLY GRADIENT HERE' CALL LIST_SCP_AVE_AND_EXACT_FLUXES (X, SCP_AVERAGES) CALL SAVE_SCP_NODE_AVERAGES (SCP_AVERAGES) IF ( SCP_2ND_DERIV ) THEN CALL LIST_SCP_AVE_EX_FLUX_GRAD (X, SCP_AVERAGES) CALL SAVE_SCP_AVE_EX_FLUX_GRAD (X, SCP_AVERAGES) END IF ! 2ND DERIVATIVES ELSE IF ( GRAD_BASE_ERROR ) THEN CALL LIST_SCP_AVE_GRADIENTS (X, SCP_AVERAGES) ELSE CALL LIST_SCP_AVE_FLUXES (X, SCP_AVERAGES) ! IF ( SCP_2ND_DERIV ) THEN !b ! CALL LIST_SCP_AVE_FLUX_GRAD (X, SCP_AVERAGES) !b ! CALL SAVE_SCP_AVE_FLUX_GRAD (SCP_AVERAGES) !b ! END IF ! 2ND DERIVATIVES !b END IF ! gradient based only CALL SAVE_SCP_NODE_AVERAGES (SCP_AVERAGES) END IF ! EXACT FLUXES GIVEN ! *** GET ELEMENT ERROR ESTIMATES VIA SCP GRADIENTS *** IF ( .NOT. NO_ERROR_EST ) THEN ! Use the SCP averages CALL SCP_ERROR_ESTIMATES (NODES, X, SCP_AVERAGES, DD, & EL_ENORM, EL_REFINE, ERR_MAX) CALL SAVE_ELEM_ERROR_ESTIMATES (EL_ENORM) CALL SAVE_PT_AVE_ERROR_EST (NODES, EL_ENORM) IF ( BAR_CHART ) THEN IBAR = 1 ; IPARM = 1 PRINT *, '' ; PRINT *,'ELEMENT ERROR ESTIMATE GRAPH' CALL ELEM_BAR_CHART (IBAR, IPARM, IPARM, EL_ENORM) END IF ! Print bar charts too END IF ! no error est requested END IF ! SCP REQUESTED END IF ! POST PROCESS THE ELEMENTS !b print *,'ITER 5, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) ! POST-PROCESS APPLICATION (AND USE SCP AVERAGING) IF ( N_FILE1 > 0 .OR. N_FILE2 > 0 ) & CALL POST_PROCESS_APPLICATION (NODES, DD, DD_OLD, ITER, X) !b CALL POST_PROCESS_APPLICATION (NODES, DD, ITER, X) !b print *,'IN_RHS ', IN_RHS !b if (IN_RHS>0) call rprint(CC_IN,1,IN_RHS,1) ! *** UPDATE VALUES FOR NEXT ITERATION (IF ANY) *** IF ( ITER == 1) THEN ; R_TEST = 1.d0 ; DD_OLD = DD ; END IF IF ( S_T_CONTINUOUS ) CALL SPACE_TIME_EBC_UPDATE & (NDEX_C_EQ, DD, COEF_C_EQ) IF ( ITER > 1) THEN WRITE (6, * ) ' ' WRITE (6, * ) 'RESULTS FOR ITERATION NUMBER = ', ITER CALL CHANGE_IN_SOLUTION (DD, DD_OLD, OLD_NORM, NEW_NORM, & DIFF_NORM, RATIO) IF ( (RATIO / R_TEST) < CUT_OFF) THEN PRINT *, 'SOLUTION HAS CONVERGED, EXIT ITERATIONS' EXIT ! ITERATION LOOP END IF CALL CORRECT_ITERATION (DD, DD_OLD) END IF !b call at(620) ; print *, 'N_REACT ', N_REACT !b END DO ! N_ITER !b print *,'ITER 6, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) ! ** BAR CHART OF SOLUTION ** IF ( BAR_CHART ) THEN PRINT *, '' ; PRINT *, 'NODAL SOLUTION GRAPHS' IBAR = 1 ! DEFAULT CONTROL DO IPARM = 1, N_G_DOF CALL NODAL_BAR_CHART (IBAR, IPARM, X, DD) END DO END IF ! Print bar charts too ! Save latest, or all results IF ( INC_SAVE > 0 ) THEN THIS_STEP = THIS_ITER !b IF ( MOD (THIS_STEP, INC_SAVE) == 0 ) THEN !b IF ( U_PLT3 > 0 ) CALL SAVE_STEP_RESULTS (DD) CALL SAVE_STEP_RESULTS (DD) IF ( U_PLT8 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_NODE_SOL (X) !b IF ( U_PLT9 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_FLUXES_AT_NODES (X) !b END IF ! desired step number END IF ! INC_SAVE > 0 IF ( SAVE_1248 ) THEN THIS_STEP = THIS_ITER !b IF ( MOD (THIS_STEP, NEXT_SAVE ) == 0 ) THEN !b replace 2 here with global NEXT_RATE NEXT_SAVE = 2 * NEXT_SAVE !b IF ( U_PLT3 > 0 ) CALL SAVE_STEP_RESULTS (DD) CALL SAVE_STEP_RESULTS (DD) IF ( U_PLT8 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_NODE_SOL (X) !b IF ( U_PLT9 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_FLUXES_AT_NODES (X) !b END IF ! desired step number END IF ! geometric progression save !b print *,'ITER 7, CC ', ITER !b call rprint(CC,1,N_D_FRE,1) END DO ! N_ITER !b ! *** PROBLEM COMPLETED, DELETE SCRATCH FILES *** ! IF ( L_REACT > 0) CLOSE (L_REACT, STATUS = "DELETE" ) IF ( N_REACT > 0) CLOSE (N_REACT, STATUS = "DELETE" ) IF ( U_FLUX > 0) CLOSE (U_FLUX , STATUS = "DELETE" ) IF ( U_INTGR > 0) CLOSE (U_INTGR, STATUS = "DELETE" ) IF ( N_FILE1 > 0) CLOSE (N_FILE1, STATUS = "DELETE" ) IF ( N_FILE2 > 0) CLOSE (N_FILE2, STATUS = "DELETE" ) !b IF ( N_FILE3 > 0) CLOSE (N_FILE3) !bc , STATUS = "DELETE" ) !b IF ( N_FILE4 > 0) CLOSE (N_FILE4) !bc , STATUS = "DELETE" ) !b IF ( N_FILE5 > 0) CLOSE (N_FILE5) !bc , STATUS = "DELETE" ) ! **** DEALLOCATE ARRAYS **** CALL DEALLOCATE_GEOM_PROP IF ( NO_PROPERTIES ) THEN DEALLOCATE ( DELETED ) ELSE CALL DEALLOCATE_PROPERTIES END IF IF ( ALLOCATED (I_DIAG) ) DEALLOCATE ( I_DIAG ) IF ( ALLOCATED (I_DOF_HI) ) DEALLOCATE ( I_DOF_HI ) IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED SKYLINE INTEGER ARRAYS' !b DEALLOCATE ( MISC_FIX, NP_FIX, SP_FIX, LP_FIX , MX_FIX) !b DEALLOCATE ( FLT_MISC, FLT_NP, FLT_SP, FLT_EL , FLT_MX) IF ( N_WARN > 0 ) PRINT *,'STATIC_MODEL GAVE ', N_WARN, ' WARNINGS' END SUBROUTINE STATIC_MODEL SUBROUTINE TRANSIENT_MODEL (X, I_BC, NODES, N_RES_EQ, COEF_C_EQ, & NDEX_C_EQ, CC_IN) ! ****************************************************** ! MODULAR PROGRAMS FOR FINITE ELEMENT ANALYSES ! TIME INTERGRATION OF M * D_DOT + K * D = C ! ****************************************************** ! DR. J. E. AKIN, P.E. ! email: akin@rice.edu ! Use System_Constants 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) REAL(DP), INTENT(IN) :: CC_IN (IN_RHS) !b 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) ! Allocatable Arrays REAL(DP), ALLOCATABLE :: M_SKY (:) ! symmetric mass REAL(DP), ALLOCATABLE :: M_SKY_LOWER (:) ! lower unsymmetric mass 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 REAL(DP) :: C_MD (N_D_FRE) !b C_MD == M_SKY * D_OLD REAL(DP) :: CC (N_D_FRE), CC_OLD (N_D_FRE) !b REAL(DP) :: DD (N_D_FRE), DD_OLD (N_D_FRE) INTEGER :: HIST_INDEX (N_G_DOF) INTEGER :: ITER = 1, IE, OLD_DEL = -1, IBAR, IPARM, IO_1 INTEGER :: COL_HI_SQ INTEGER :: J, T_G, T_S ! TIME GROUP AND TIME STEP NUMBERS INTEGER, SAVE :: METHOD = 1 ! time step rule !b LOGICAL :: INITIAL = .TRUE. !b ! SEE NOTATION.F FOR VARIABLE AND ARRAY NOTATIONS IF ( I_BUG > 0 ) THEN PRINT *,'=== USING TRANSIENT_MODEL ===' !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 ! 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) COL_HI_SQ = SUM ( I_DOF_HI **2 ) 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, N_STORE, COL_HI_SQ !b 5 FORMAT (/, 'TOTAL NUMBER OF SYSTEM EQUATIONS .....',I10,& !b /, 'NUMBER OF STIFFNESS COEFFICIENTS .....',I10,& !b /, 'SUM OF HEIGHTS SQUARED ...............',I10) M_COEFF = N_COEFF !b IF ( DIAGONAL_SYS_MASS ) M_COEFF = N_D_FRE !b !b WRITE (6, 6) M_COEFF !b 6 FORMAT ( 'TOTAL NUMBER OF MASS COEFFICIENTS ...',I10) ! ALLOCATE AND ZERO THE SYSTEM MATRIICES ALLOCATE ( M_SKY(M_COEFF) ) ! add status M_SKY = 0.0_DP ! zero mass ALLOCATE ( S_SKY(N_COEFF) ) ! add status S_SKY = 0.0_DP ! zero stiffness IF ( .NOT. SYMMETRIC ) THEN ALLOCATE ( S_SKY_LOWER(NOT_SYM) ) ! add status S_SKY_LOWER = 0.0_DP ! zero stiffness ALLOCATE ( M_SKY_LOWER (M_COEFF) ) ! add status M_SKY_LOWER = 0.d0 END IF IF ( I_BUG > 0 ) PRINT *, ' ALLOCATED SKYLINE REAL ARRAYS' ! *** BEGIN TIME INTEGRATION GROUPS *** CC_OLD = 0.d0 !b DD_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: TRANSIENT_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 !b T_STRING = TO3STR (T_G) ! integer to string 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 !b T_STRING = TO3STR (T_G) ! integer to string 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) CC_OLD = CC !b ! *** ASSEMBLY COMPLETED, CHECK SOURCES *** IF ( DEBUG_ASSEMBLY .AND. N_D_FRE < M_SHOW .AND. T_G == 1 ) THEN PRINT *,'DEBUG PRINT AFTER ASSEMBLY IN FULL MATRIX FORM:' 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 PRINT *,'M_SKY at 1104' CALL SKY_FUL_GEN (M_SKY, I_DIAG, CC, M_SKY_LOWER) ! show values CALL SUM_RHS (CC) END IF !--> ** APPLY BOUNDARY CONSTRAINTS TO NODAL PARAMETERS ** IF ( I_BUG > 0 ) PRINT *, 'Before SKY_BC_GEN' ! KNOWN VALUES TO RHS CALL SKY_BC_GEN (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, S_SKY, CC, & I_DIAG, S_SKY_LOWER) IF ( DEBUG_ASSEMBLY .AND. N_D_FRE < M_SHOW ) THEN ! PRINT FULL SIZE PRINT *,'DEBUG PRINT AFTER BC APPLICATION' CALL SKY_FUL_GEN (S_SKY, I_DIAG, CC, S_SKY_LOWER) END IF ! DEBUG PRINT !--> *** FACTOR THE LHS SQUARE MATRIX, ONCE PER TIME GROUP *** FACT = .TRUE. ; BACK = .FALSE. IF ( SYMMETRIC ) THEN ! assume positive definite CALL SKY_SOLVE_SYM (S_SKY, CC, DD, I_DIAG, FACT, BACK) ELSE CALL SKY_SOLVE_UNSYM (S_SKY, S_SKY_LOWER, CC, DD, I_DIAG, & FACT, BACK) END IF ! *** TIME STEPPING WITHIN EACH TIME GROUP *** IF ( DEBUG_UNITS ) THEN CALL LIST_SYSTEM_UNIT_NUMBERS CALL LIST_OPEN_IO_UNITS END IF ! Now THIS_STEP = 0 or RESTART_STEP STEPS: DO T_S = 1, T_STEPS THIS_STEP = THIS_STEP + 1 TIME = START_TIME + T_S * TIME_STEP !b PRINT *,' TIME STEP ', T_S, ' AT TIME ', TIME PRINT *,' TIME STEP ', THIS_STEP, ' AT TIME ', TIME !--> *** UPDATE SOURCE VECTOR FROM PREVIOUS STEP *** ! Form C_MD == M_SKY * D_OLD, THEN UPDATE BC IF ( SYMMETRIC ) THEN CALL VECTOR_MULT_SKY_SYM (M_SKY, DD, I_DIAG, C_MD) ELSE CALL VECTOR_MULT_SKY_GEN (M_SKY, DD, I_DIAG, M_SKY_LOWER, C_MD) END IF ! Zero rows at essential boundary conditions CALL TYPE_1_COL_ZERO (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, C_MD) IF ( DEBUG_ASSEMBLY .AND. N_D_FRE < M_SHOW ) THEN PRINT *,'DEBUG PRINT AFTER BC APPLICATION TO C_MD & CC' CALL RPRINT (C_MD, 1, N_D_FRE, 1) CALL RPRINT (CC, 1, N_D_FRE, 1) END IF ! DEBUG PRINT !--> *** SOLVE FOR UNKNOWN NODAL PARAMETERS *** FACT = .FALSE. ; BACK = .TRUE. IF ( SYMMETRIC ) THEN ! assume positive definite CALL SKY_SOLVE_SYM (S_SKY, (C_MD + CC), DD, I_DIAG, & FACT, BACK) ELSE CALL SKY_SOLVE_UNSYM (S_SKY, S_SKY_LOWER, (C_MD + CC),& DD, I_DIAG, FACT, BACK) END IF !--> *** SOLUTION COMPLETE, GET SYSTEM REACTIONS *** ! GET REACTION ROWS FILE IF TRANSIENT REACTIONS DESIRED !b like SAVE_SKY_ROWS_GEN for use in GET_REACTIONS !b IF ( N_REACT > 0 .AND. SUM(I_BC) > 0) THEN !b CALL GET_REACTIONS (DD) !b IF ( SAVE_REACTIONS ) CALL SAVE_STEP__REACTIONS (DD) !b END IF ! output/save reactions ! *** PRINT RESULTS *** CALL MAX_AND_MIN_RANGE (DD) IF ( NPT_WRT == 1) THEN IF ( USE_EXACT ) THEN CALL WRITE_BY_NODES_WITH_EXACT (X, DD) ELSE CALL WRITE_BY_NODES (X, DD) END IF ! EXACT SOLUTION PROVIDED (IN EXACT_SOLUTION) END IF ! LIST BY NODES IF ( HISTORY_DOF > 0 ) THEN ! update history plot WRITE (U_DOF, '(I6, 2(1PE13.5))') HISTORY_NODE, TIME, & DD (HISTORY_DOF) END IF IF ( HISTORY_NODE > 0 ) THEN WRITE (U_NODE, '( I6, (8(1PE13.5)) )') HISTORY_NODE, TIME, & DD (HIST_INDEX) END IF IF ( BAR_CHART ) THEN PRINT *, '' ; PRINT *, 'NODAL SOLUTION GRAPHS' IBAR = 1 ! DEFAULT CONTROL DO IPARM = 1, N_G_DOF CALL NODAL_BAR_CHART (IBAR, IPARM, X, DD) END DO END IF ! Print bar charts too DD_OLD = DD ! Save latest, or all results IF ( INC_SAVE > 0 ) THEN IF ( MOD (THIS_STEP, INC_SAVE) == 0 ) THEN !b IF ( U_PLT3 > 0 ) CALL SAVE_STEP_RESULTS (DD) CALL SAVE_STEP_RESULTS (DD) IF ( U_PLT8 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_NODE_SOL (X) !b IF ( U_PLT9 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_FLUXES_AT_NODES (X) !b END IF ! desired step number END IF ! INC_SAVE > 0 IF ( SAVE_1248 ) THEN IF ( MOD (THIS_STEP, NEXT_SAVE ) == 0 ) THEN !b replace 2 below with global NEXT_RATE NEXT_SAVE = 2 * NEXT_SAVE !b IF ( U_PLT3 > 0 ) CALL SAVE_STEP_RESULTS (DD) CALL SAVE_STEP_RESULTS (DD) IF ( U_PLT8 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_NODE_SOL (X) !b IF ( U_PLT9 > 0 .AND. USE_EXACT ) CALL & SAVE_STEP_EXACT_FLUXES_AT_NODES (X) !b END IF ! desired step number END IF ! geometric progression save END DO STEPS ! WITHIN A TIME GROUP END DO GROUPS ! OF TIME INTERVALS IF ( .NOT. SYMMETRIC ) DEALLOCATE ( S_SKY_LOWER ) DEALLOCATE ( S_SKY ) DEALLOCATE ( M_SKY ) IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED SKYLINE REAL ARRAYS' ! *** PROBLEM COMPLETED, DELETE SCRATCH FILES *** IF ( L_REACT > 0) CLOSE (L_REACT, STATUS = "DELETE" ) IF ( N_REACT > 0) CLOSE (N_REACT, STATUS = "DELETE" ) IF ( U_FLUX > 0) CLOSE (U_FLUX , STATUS = "DELETE" ) !b IF ( N_FILE1 > 0) CLOSE (N_FILE1, STATUS = "DELETE" ) IF ( N_FILE2 > 0) CLOSE (N_FILE2, STATUS = "DELETE" ) ! **** DEALLOCATE ARRAYS **** CALL DEALLOCATE_GEOM_PROP CALL DEALLOCATE_PROPERTIES IF ( ALLOCATED (I_DIAG) ) DEALLOCATE ( I_DIAG ) IF ( ALLOCATED (I_DOF_HI) ) DEALLOCATE ( I_DOF_HI ) IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED SKYLINE INTEGER ARRAYS' IF ( N_WARN > 0 ) PRINT *,'TRANSIENT_MODEL GAVE ', N_WARN, ' WARNINGS' END SUBROUTINE TRANSIENT_MODEL