! File: exact_lib.f ! copyright 2005, J. E. Akin, all rights reserved. ! A library of exact analytic solutions used for testing and ! educational purposes. The are listed in order of the allowed ! global integer value of EXACT_CASE. The default value is ! EXACT_CASE = 0, which means the needed code is in various ! include files that were compiled before running MODEL. ! --------------------- EXACT_CASE = 0 ------------------------ SUBROUTINE EXACT_SOLUTION (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE include 'my_exact_inc' IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE EXACT_SOLUTION SUBROUTINE EXACT_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE include 'my_exact_flux_inc' IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE EXACT_FLUX SUBROUTINE EXACT_INTEGRAL (XYZ, INTEGRAL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT INTEGRAL FOR ELEMENT MATRICES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: INTEGRAL ! PT EXACT VALUE include 'my_exact_integral_inc' IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE EXACT_INTEGRAL SUBROUTINE EXACT_NORMAL_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX ! PT EXACT VALUE include 'my_exact_normal_flux_inc' IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE EXACT_NORMAL_FLUX SUBROUTINE EXACT_SOURCE (XYZ, SOURCE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOURCE FOR ELEMENT MATRICES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE include 'my_exact_source_inc' IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE EXACT_SOURCE SUBROUTINE EXACT_FLUX_GRAD (XYZ, FLUX_GRAD) ! 2nd deriv ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET EXACT FLUX GRADIENTS (2ND DERIVATIVES) AT A NODE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for MAX_NP, N_R_B, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT(IN) :: XYZ (N_SPACE) REAL(DP), INTENT(OUT) :: FLUX_GRAD (N_R_B * N_SPACE) !b INTEGER :: HIGHER ! FLUX_GRAD = EXACT FLUX GRAD (2ND DERIV) COMPONENTS AT A POINT ! HIGHER = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS !b HIGHER = N_R_B * N_SPACE include 'my_exact_2nd_der_inc' IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE EXACT_FLUX_GRAD ! ------------------------------------------------------------------ ! ------------------------------------------------------------------ ! TABULATION OF EXACT SOLUTION CASES ! ------------------------------------------------------------------ ! ------------------------------------------------------------------ ! --------------------- EXACT_CASE = 1 ------------------------ SUBROUTINE MYERS_1D_ROD_HEAT (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Myers example conduction and convection along a rod, tip insulated ! K*A*U,XX - h*P*(U-U_ext) = 0, U(0)=U_0, dU/dx(L)=0 ! For EXACT_CASE = 1, EXAMPLE = 101, DATA_SET = 1, or 2 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! For globally constant data the analytic solution is: ! U(x) = U_ext - (U_0-U_ext) * Cosh (m*(L-x)) / Cosh (mL) ! m^2 = h_e*P_e/(K_e*A_e) Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), SAVE :: K_e, A_e, h_e, P_e, L, m, Cosh_mL, U_0, U_ext LOGICAL, SAVE :: FIRST_CALL = .TRUE. IF ( FIRST_CALL ) THEN ! one time calculations FIRST_CALL = .FALSE. IF ( DEBUG_PROPERTY ) PRINT *, 'Entering MYERS_1D_ROD_HEAT' IF ( EXACT_REALS > 0 ) THEN ! (can use_exact_bc) K_e = GET_REAL_EXACT (1) ! thermal conductivity A_e = GET_REAL_EXACT (2) ! area of bar h_e = GET_REAL_EXACT (3) ! convection coefficent on perimeter P_e = GET_REAL_EXACT (4) ! perimeter of area A_e U_ext = GET_REAL_EXACT (5) ! external ref temp L = GET_REAL_EXACT (6) ! for exact U_0 = GET_REAL_EXACT (7) ! essential bc at x = 0 ELSE ! use el and misc prop (but not for bc assignment) K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e U_ext = GET_REAL_MISC (1) ! external ref temp L = GET_REAL_MISC (2) ! for exact U_0 = GET_REAL_MISC (3) ! essential bc at x = 0 END IF ! source of data m = sqrt (h_e*P_e/(K_e*A_e)) ! for exact sol COSH_ML = COSH (m*L) ! for exact IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED MYERS_1D_ROD_HEAT' IF ( N_SPACE /= 1 .OR. N_G_DOF /= 1 ) & STOP 'MYERS_1D_ROD_HEAT: INVALID CONSTANTS' END IF ! first call EXACT_SOL (1) = U_ext + (U_0-U_ext) * Cosh (m*(L-XYZ(1))) / COSH_ML END SUBROUTINE MYERS_1D_ROD_HEAT SUBROUTINE MYERS_1D_ROD_HEAT_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Myers example conduction and convection along a rod, tip insulated ! K*A*U,XX - h*P*(U-U_ext) = 0, U(0)=U_0, dU/dx(L)=0 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE LOGICAL, SAVE :: FIRST_CALL = .TRUE. REAL(DP), SAVE :: K_e, A_e, h_e, P_e, L, m, Cosh_mL, U_0, U_ext IF ( FIRST_CALL ) THEN ! one time calculations FIRST_CALL = .FALSE. IF ( DEBUG_PROPERTY ) PRINT *, 'Entering MYERS_1D_ROD_HEAT_FLUX' IF ( EXACT_REALS > 0 ) THEN ! (can use_exact_bc) K_e = GET_REAL_EXACT (1) ! thermal conductivity A_e = GET_REAL_EXACT (2) ! area of bar h_e = GET_REAL_EXACT (3) ! convection coefficent on perimeter P_e = GET_REAL_EXACT (4) ! perimeter of area A_e U_ext = GET_REAL_EXACT (5) ! external ref temp L = GET_REAL_EXACT (6) ! for exact U_0 = GET_REAL_EXACT (7) ! essential bc at x = 0 ELSE ! use el and misc prop (but not for bc assignment) K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e U_ext = GET_REAL_MISC (1) ! external ref temp L = GET_REAL_MISC (2) ! for exact U_0 = GET_REAL_MISC (3) ! essential bc at x = 0 END IF ! source of data m = sqrt (h_e*P_e/(K_e*A_e)) ! for exact sol COSH_ML = COSH (m*L) ! for exact IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED MYERS_1D_ROD_HEAT_FLUX' IF ( N_SPACE /= 1 .OR. N_G_DOF /= 1 .OR. N_R_B /= 1 ) & STOP 'MYERS_1D_ROD_HEAT_FLUX: INVALID CONSTANTS' END IF ! first call FLUX (1) = -(U_0-U_ext)*m*SINH (m*(L-XYZ(1))) / COSH_ML ! exact grad END SUBROUTINE MYERS_1D_ROD_HEAT_FLUX SUBROUTINE MYERS_1D_ROD_HEAT_LOSS (XYZ, LOSS) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Myers example conduction and convection along a rod, tip insulated ! K*A*U,XX - h*P*(U-U_ext) = 0, U(0)=U_0, dU/dx(L)=0 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: LOSS (N_R_B) ! PT EXACT VALUE LOGICAL, SAVE :: FIRST_CALL = .TRUE. REAL(DP), SAVE :: K_e, A_e, h_e, P_e, L, m, Cosh_mL, U_0, U_ext IF ( FIRST_CALL ) THEN ! one time calculations FIRST_CALL = .FALSE. IF ( DEBUG_PROPERTY ) PRINT *, 'Entering MYERS_1D_ROD_HEAT_LOSS' IF ( EXACT_REALS > 0 ) THEN ! (can use_exact_bc) K_e = GET_REAL_EXACT (1) ! thermal conductivity A_e = GET_REAL_EXACT (2) ! area of bar h_e = GET_REAL_EXACT (3) ! convection coefficent on perimeter P_e = GET_REAL_EXACT (4) ! perimeter of area A_e U_ext = GET_REAL_EXACT (5) ! external ref temp L = GET_REAL_EXACT (6) ! for exact U_0 = GET_REAL_EXACT (7) ! essential bc at x = 0 ELSE ! use el and misc prop (but not for bc assignment) K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e U_ext = GET_REAL_MISC (1) ! external ref temp L = GET_REAL_MISC (2) ! for exact U_0 = GET_REAL_MISC (3) ! essential bc at x = 0 END IF ! source of data m = sqrt (h_e*P_e/(K_e*A_e)) ! for exact sol COSH_ML = COSH (m*L) ! for exact IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED MYERS_1D_ROD_HEAT_LOSS' IF ( N_SPACE /= 1 .OR. N_G_DOF /= 1 .OR. N_R_B /= 1 ) & STOP 'MYERS_1D_ROD_HEAT_LOSS: INVALID CONSTANTS' END IF ! first call LOSS = H_e*P_e*(U_0-U_ext)*Sinh(m*L)/(m*Cosh_mL) IF ( TOUCH ) PRINT *, XYZ(1) ! (touch is never true) END SUBROUTINE MYERS_1D_ROD_HEAT_LOSS ! --------------------- EXACT_CASE = 2 ------------------------ SUBROUTINE LAKHANY_WHITEMAN_TEST (XYZ, EXACT_SOL) ! "Superconvergeny Derivative Operators: Derivative Recovery ! Techniques", A. M. Lakhany, J. R. Whiteman, pp. 196-216, ! Finite Element Methods, M. Krizek, et. al (Eds), Marcel ! Decker, New York, 1998, EXACT_CASE = 2 ! (u,xx+u,yy) = -(2 + pi_sq * (1-x^2)) * Sin (pi y) ! on (-1,-1)X(1,1). Exact u(x, y) = (1 - x^2) * Sin (pi y) ! Symmetric about x=0, anti-symmetric about y=0 ! A second derivative test problem Use System_Constants ! for N_G_DOF, N_SPACE IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE !b REAL(DP), PARAMETER :: PI = 3.141592654d0 !b REAL(DP), PARAMETER :: TWO_PI = (PI + PI), PI_SQ = (PI * PI) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: LAKHANY_WHITEMAN_TEST' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'LAKHANY_WHITEMAN_TEST: INVALID CONSTANTS' END IF EXACT_SOL (1) = (1.d0 - XYZ(1)**2) * SIN (PI * XYZ(2)) END SUBROUTINE LAKHANY_WHITEMAN_TEST SUBROUTINE LAKHANY_WHITEMAN_TEST_SOURCE (XYZ, SOURCE) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT SOURCE !b REAL(DP), PARAMETER :: PI = 3.141592654d0 !b REAL(DP), PARAMETER :: TWO_PI = (PI + PI), PI_SQ = (PI * PI) INTEGER, SAVE :: NOTE = 0, COUNT = 0, LIMIT = 10 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: LAKHANY_WHITEMAN_TEST_SOURCE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'LAKHANY_WHITEMAN_TEST_SOURCE: INVALID CONSTANTS' END IF SOURCE = (2.d0 + PI_SQ * (1.d0 - XYZ(1)**2)) * SIN (PI * XYZ(2)) COUNT = COUNT + 1 IF ( DEBUG_EL_COL .AND. COUNT <= LIMIT ) PRINT *,'PT, SOURCE: ', & XYZ, SOURCE END SUBROUTINE LAKHANY_WHITEMAN_TEST_SOURCE SUBROUTINE LAKHANY_WHITEMAN_TEST_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE !b REAL(DP), PARAMETER :: PI = 3.141592654d0 !b REAL(DP), PARAMETER :: TWO_PI = (PI + PI), PI_SQ = (PI * PI) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: LAKHANY_WHITEMAN_TEST_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'LAKHANY_WHITEMAN_TEST_FLUX: INVALID CONSTANTS' END IF FLUX (1) = -2.d0 * XYZ(1) * SIN (PI * XYZ(2)) FLUX (2) = (1.d0 - XYZ(1)**2) * PI * COS (PI * XYZ(2)) END SUBROUTINE LAKHANY_WHITEMAN_TEST_FLUX SUBROUTINE LAKHANY_WHITEMAN_TEST_DERIV_2 (XYZ, DERIV_2) Use System_Constants ! for N_G_DOF, N_SPACE, N_2_DER, EXACT_CASE ! Exact derivative orders are: U,xx U,xy U,yy ! FEA derivative orders are: U,xx U,yx U,xy U,yy IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: DERIV_2 (N_2_DER) ! PT EXACT VALUE !b REAL(DP), PARAMETER :: PI = 3.141592654d0 !b REAL(DP), PARAMETER :: TWO_PI = (PI + PI), PI_SQ = (PI * PI) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: LAKHANY_WHITEMAN_TEST_DERIV_2' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'LAKHANY_WHITEMAN_TEST_DERIV_2: INVALID CONSTANTS' END IF DERIV_2 (1) = -2.d0 * SIN (PI * XYZ(2)) ! XX DERIV_2 (2) = -TWO_PI * XYZ(1) * COS (PI * XYZ(2)) ! YX = XY DERIV_2 (3) = -PI_SQ * (1.d0 - XYZ(1)**2) * SIN (PI * XYZ(2)) !YY END SUBROUTINE LAKHANY_WHITEMAN_TEST_DERIV_2 SUBROUTINE LAKHANY_WHITEMAN_TEST_DERIV_3 (XYZ, DERIV_3) Use System_Constants ! for N_G_DOF, N_SPACE, N_3_DER, EXACT_CASE ! Exact derivative orders are: U,xxx U,xxy U,xyy U,yyy ! FEA derivative orders are: IMPLICIT NONE !b INTEGER, PARAMETER :: N_3_DER = 4 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: DERIV_3 (N_3_DER) ! PT EXACT VALUE !b REAL(DP), PARAMETER :: PI = 3.141592654d0 !b REAL(DP), PARAMETER :: TWO_PI = (PI + PI), PI_SQ = (PI * PI) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: LAKHANY_WHITEMAN_TEST_DERIV_3', XYZ IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'LAKHANY_WHITEMAN_TEST_DERIV_3: INVALID CONSTANTS' END IF DERIV_3 = 0.d0 STOP 'LAKHANY_WHITEMAN_TEST_DERIV_3: NOT YET CODED' END SUBROUTINE LAKHANY_WHITEMAN_TEST_DERIV_3 ! --------------------- EXACT_CASE = 3 ------------------------ SUBROUTINE POISSON_ZZ_SQ_W_SOURCE (XYZ, EXACT_SOL) ! Isotropic Poisson Eq on a square from (-1, -1) to (+1, +1) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 3 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; PRINT *, 'NOTE: USED POISSON_ZZ_SQ_W_SOURCE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'POISSON_ZZ_SQ_W_SOURCE: INVALID CONSTANTS' END IF EXACT_SOL (1) = (1.d0 - XYZ(1)**2) * (1.d0 - XYZ(2)**2) END SUBROUTINE POISSON_ZZ_SQ_W_SOURCE SUBROUTINE POISSON_ZZ_SQ_W_SOURCE_FLUX (XYZ, FLUX) ! Isotropic Poisson Eq on a square from (-1, -1) to (+1, +1) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 3 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: POISSON_ZZ_SQ_W_SOURCE_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'POISSON_ZZ_SQ_W_SOURCE_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL = (1.d0 - XYZ(1)**2) * (1.d0 - XYZ(2)**2) FLUX (1) = -2.d0 * XYZ(1) * (1.d0 - XYZ(2)**2) FLUX (2) = -(1.d0 - XYZ(1)**2) * 2.d0 * XYZ(2) END SUBROUTINE POISSON_ZZ_SQ_W_SOURCE_FLUX SUBROUTINE POISSON_ZZ_SQ_W_SOURCE_SOURCE (XYZ, SOURCE) ! Isotropic Poisson Eq on a square from (-1, -1) to (+1, +1) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 3 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT SOURCE INTEGER :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; PRINT *, 'NOTE: POISSON_ZZ_SQ_W_SOURCE_SOURCE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'POISSON_ZZ_SQ_W_SOURCE_SOURCE: INVALID CONSTANTS' END IF SOURCE = -2.d0 * (XYZ(1)**2 + XYZ(2)**2 - 2.d0) END SUBROUTINE POISSON_ZZ_SQ_W_SOURCE_SOURCE ! --------------------- EXACT_CASE = 4 ------------------------ SUBROUTINE STRONG_DIAGONAL_SQ (XYZ, EXACT_SOL) ! Solution of the 2-D test case used by Oden and Zienkwicz ! a square from (0,0) to (1,1) with strong diagonal gradients ! The solution is symmetric about y = x ! Comp Meth App Mech Engrg 77 (1989) 79-212 ! I J Num Meth Engrg 33 (1992) 207-224 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 4 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP) :: PSI, X, Y REAL(DP), PARAMETER :: TEMP = 20.d0, PSI_O = 0.8d0 REAL(DP), PARAMETER :: ROOT_2 = 1.4142135623730950488d0 INTEGER :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; PRINT *, 'NOTE: STRONG_DIAGONAL_SQ' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'STRONG_DIAGONAL_SQ: INVALID CONSTANTS' END IF X = XYZ (1) ; Y = XYZ (2) ; PSI = (X + Y) / ROOT_2 EXACT_SOL (1) = X * (1.d0 - X) * Y * (1.d0 - Y) & * ATAN ( TEMP * (PSI - PSI_O) ) END SUBROUTINE STRONG_DIAGONAL_SQ SUBROUTINE STRONG_DIAGONAL_SQ_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 4 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP) :: X, Y INTEGER :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; PRINT *, 'NOTE: STRONG_DIAGONAL_SQ_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'STRONG_DIAGONAL_SQ_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL = X*(1.d0 - X)*Y*(1.d0 - Y)*ATAN(TEMP*(PSI-PSI_O)) ! PSI = (X + Y) / ROOT_2 ! PARAMETER :: TEMP = 20.d0, PSI_O = 0.8d0 ! PARAMETER :: ROOT_2 = 1.4142135623730950488d0 X = XYZ (1) ; Y = XYZ (2) ! Here FLUX = GRADIENT FLUX (1) = (1.d0 - Y)*Y*(((14.14213562373094d0 & - 14.14213562373095d0*X)*X) & / (257.d0 + 200.d0*(-2.262741699796952d0 + X)*X & - 452.5483399593904d0*Y + 400.d0*X*Y + 200.d0*Y**2) & + (-1.d0 + 2.d0*X)*ATAN(16.d0 - 14.14213562373095d0*X & - 14.14213562373095d0*Y)) ! can be simplified FLUX (2) = (1.d0 - X)*X*(((14.14213562373094d0 & - 14.14213562373095d0*Y)*Y) & / (257.d0 + 200.d0*(-2.262741699796952d0 + X)*X & - 452.5483399593904d0*Y + 400.d0*X*Y + 200.d0*Y**2) & + (-1.d0 + 2.d0*Y)*ATAN(16.d0 - 14.14213562373095d0*X & - 14.14213562373095d0*y)) END SUBROUTINE STRONG_DIAGONAL_SQ_FLUX SUBROUTINE STRONG_DIAGONAL_SQ_SOURCE (XYZ, SOURCE) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 4 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT SOURCE REAL(DP) :: X, Y INTEGER :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1 PRINT *, 'NOTE: STRONG_DIAGONAL_SQ_SOURCE' END IF IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'STRONG_DIAGONAL_SQ_SOURCE: INVALID CONSTANTS' X = XYZ (1) ; Y = XYZ (2) SOURCE = (X**4*(5656.854249492382d0 - 11313.70849898476d0*Y) & - 22627.41699796952d0*(X**3)*(-1.631370849898475d0 & + Y)*(-0.5d0 + Y) - 22627.41699796952d0*(X**2) & * (-0.4779882018197526d0 + Y)* (1.855559240944815d0 & + (-2.719068073027961d0 + Y)*Y) - 11313.70849898476d0*X & * (-1.707106781186555d0 + Y)* (-0.2928932188134525d0 & + Y)* (1.284999999999993d0 + (-2.262741699796944d0 + Y) & * Y) + 5656.854249492382d0*(-1.d0 + Y)*Y*(1.285d0 & + (-2.262741699796912d0 + Y)*Y))/(257.d0 + 200.d0 & * (-2.262741699796952d0 + X)*X - 452.5483399593904d0*Y & + 400.d0*X*Y + 200.d0*Y**2)**2 + 2.d0*((-1.d0 + X)*X & + (-1.d0 + Y)*Y)*ATAN(16.d0 - 14.14213562373095d0*X & - 14.14213562373095d0*Y) END SUBROUTINE STRONG_DIAGONAL_SQ_SOURCE ! --------------------- EXACT_CASE = 5 ------------------------ SUBROUTINE POISSON_ANNULUS_ARC (XYZ, EXACT_SOL) ! Solution of the 2-D Poisson Equation on an segment of a circular ! annulus in the first quadrant from r=1 to r=2 ! r = 2 |-----\ ! q_n = 0..| \ .... q_n = 4 X^2 Y^2 ! r = 1 |--\ K,Q \ Q = 4(X^2 + y^2) = 4 R^2 ! T=10-X^2+X^4...\ | K = 2 ! Y | | T(X,Y) = 10 - X^2 Y^2, exact ! o-X ------- ... T = 10 (10 is arbitrary constant) ! Dirchlet BC: constant at Y=0, variable at R=1 ! Normal flux: variable on R = 2, Zero on X = 0 (or on X = Y) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 5 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP) :: X, Y INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED POISSON_ANNULUS_ARC' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'POISSON_ANNULUS_ARC: INVALID CONSTANTS' END IF X = XYZ (1) ; Y = XYZ (2) EXACT_SOL (1) = 10.d0 - X**2 * Y**2 END SUBROUTINE POISSON_ANNULUS_ARC SUBROUTINE POISSON_ANNULUS_ARC_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 5 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP) :: X, Y, K = 2.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED ' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'POISSON_ANNULUS_ARC_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL = 10.d0 - X**2 * Y**2 X = XYZ (1) ; Y = XYZ (2) FLUX (1) = -2.d0 * X * Y**2 ! Gradient FLUX (2) = -2.d0 * X**2 * Y ! Gradient FLUX = FLUX * K END SUBROUTINE POISSON_ANNULUS_ARC_FLUX SUBROUTINE POISSON_ANNULUS_ARC_DERIV_2 (XYZ, DERIV_2) Use System_Constants ! for N_G_DOF, N_SPACE, N_2_DER, EXACT_CASE = 5 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: DERIV_2 (N_2_DER) ! PT EXACT VALUE REAL(DP) :: X, Y, K = 2.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: POISSON_ANNULUS_ARC_DERIV_2' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'POISSON_ANNULUS_ARC_DERIV_2: INVALID CONSTANTS' END IF ! EXACT_SOL = 10.d0 - X**2 * Y**2 X = XYZ (1) ; Y = XYZ (2) DERIV_2 (1) = -2.d0 * Y**2 ! EXACT,XX DERIV_2 (2) = -4.d0 * X * Y ! EXACT,XY == YX DERIV_2 (3) = -2.d0 * X**2 ! EXACT,YY DERIV_2 = DERIV_2 * K END SUBROUTINE POISSON_ANNULUS_ARC_DERIV_2 SUBROUTINE POISSON_ANNULUS_ARC_NORMAL_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 5 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX ! PT EXACT VALUE REAL(DP) :: X, Y, K = 2.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: POISSON_ANNULUS_ARC_NORMAL_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'POISSON_ANNULUS_ARC_NORMAL_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL = 10.d0 - X**2 * Y**2 X = XYZ (1) ; Y = XYZ (2) FLUX = -4.d0 * K * (X**2 * Y**2) / SQRT (X**2 + Y**2) ! on arc ! Note: on a constant radius R the normal flux is ! q_n = - 4.d0 * K * (X**2 * Y**2) / R ! For constant A angle line the normal flux is ! q_n = - 2.d0 * K * X * Y * (X * Sin A - Y * Cos A) END SUBROUTINE POISSON_ANNULUS_ARC_NORMAL_FLUX SUBROUTINE POISSON_ANNULUS_ARC_SOURCE (XYZ, SOURCE) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 5 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT SOURCE REAL(DP) :: X, Y, K = 2.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED POISSON_ANNULUS_ARC_SOURCE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'POISSON_ANNULUS_ARC_SOURCE: INVALID CONSTANTS' END IF X = XYZ (1) ; Y = XYZ (2) SOURCE = 2.d0 * K * (X**2 + Y**2) END SUBROUTINE POISSON_ANNULUS_ARC_SOURCE ! --------------------- EXACT_CASE = 6 ------------------------ SUBROUTINE NAFEMS_CONVECTION_TEST (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 6 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! H = 750 W/mC, 0 C L_AB = 0.6 m, L_BC = 1.0 m ! D----------C AT (0.6, 0.2) T = 18.254 C ! q=0 | | H = 750 W/mC, 0 C ! | k=52W/mC | ! A----------B --> x Point B has very high gradients ! T_AB = 100 C INTEGER :: I INTEGER, PARAMETER :: MAX_ROOT = 40 REAL (dp) :: X, Y, TEMP, TEMP1, TEMP2 REAL (dp), PARAMETER :: A = 0.6d0, B = 1.d0, H = 750.d0, K = 52.d0 REAL (dp), PARAMETER :: T_AB = 100.d0, H_BAR = H/K REAL (dp), PARAMETER :: H_BAR_SQ = H_BAR**2 ! ROOT IS ROOT OF ROOT * TAN(ROOT*A) - H_BAR = 0 REAL (dp), PARAMETER :: ROOT (MAX_ROOT) = (/ & 2.3489248096585258d0, 7.0922953913044458d0, 11.937651219167781d0,& 16.88611595161343d0, 21.914106119234724d0, 26.997701929604876d0,& 32.119360643037886d0, 37.26735387646562d0, 42.433975347748458d0,& 47.614106330819361d0, 52.804274611335842d0, 58.002068114386958d0,& 63.205771294262838d0, 68.414136745272785d0, 73.626238908565369d0,& 78.841378309935209d0, 84.059017498643797d0, 89.278737291769232d0,& 94.500206298463695d0, 99.723159305319967d0, 104.94738168756264d0,& 110.17269799133665d0, 115.39896345131892d0, 120.62605760591623d0,& 125.85387943287243d0, 131.08234360156550d0, 136.31137755558950d0,& 141.54091921974700d0, 146.77091518163468d0, 152.00131923753867d0,& 157.23209122057966d0, 162.46319604943233d0, 167.69460295082322d0,& 172.92628481998585d0, 178.15821769142210d0, 183.39038029845793d0,& 188.62275370473256d0, 193.85532099431336d0, 199.08806700986247d0,& 204.32097813040141d0 /) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED NAFEMS_CONVECTION_TEST' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /= 2 ) & STOP 'NAFEMS_CONVECTION_TEST: INVALID CONSTANTS' END IF X = XYZ(1) ; Y = XYZ(2) ; TEMP = 0d0 DO I = 1, MAX_ROOT TEMP1 = COS(ROOT(I)*X) * (ROOT(I) * COSH(ROOT(I)*(B-Y)) & + H_BAR * SINH(ROOT(I)*(B-Y)) ) TEMP2 = ( (ROOT(I)**2 + H_BAR_SQ) * A + H_BAR ) & * ( ROOT(I)*COSH(ROOT(I)*B) + H_BAR*SINH(ROOT(I)*B) ) & * COS(ROOT(I)*A) TEMP = TEMP + TEMP1/TEMP2 END DO ! TO APPROXIMATE INFINITE SERIES EXACT_SOL = 2.d0 * H_BAR * T_AB * TEMP END SUBROUTINE NAFEMS_CONVECTION_TEST SUBROUTINE NAFEMS_CONVECTION_TEST_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 6 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE ! H = 750 W/mC, 0 C L_AB = 0.6 m, L_BC = 1.0 m ! D----------C AT (0.6, 0.2) T = 18.254 C ! q=0 | | H = 750 W/mC, 0 C ! | k=52W/mC | ! A----------B --> x Point B has very high gradients ! T_AB = 100 C INTEGER :: I INTEGER, PARAMETER :: MAX_ROOT = 40 REAL (dp) :: X, Y, TEMP1, TEMP2 REAL (dp) :: DX_TEMP, DX_TEMP1 REAL (dp) :: DY_TEMP, DY_TEMP1 REAL (dp), PARAMETER :: A = 0.6d0, B = 1.d0, H = 750.d0, K = 52.d0 REAL (dp), PARAMETER :: T_AB = 100.d0, H_BAR = H/K REAL (dp), PARAMETER :: H_BAR_SQ = H_BAR**2 ! ROOT IS ROOT OF ROOT * TAN(ROOT*A) - H_BAR = 0 REAL (dp), PARAMETER :: ROOT (MAX_ROOT) = (/ & 2.3489248096585258d0, 7.0922953913044458d0, 11.937651219167781d0,& 16.88611595161343d0, 21.914106119234724d0, 26.997701929604876d0,& 32.119360643037886d0, 37.26735387646562d0, 42.433975347748458d0,& 47.614106330819361d0, 52.804274611335842d0, 58.002068114386958d0,& 63.205771294262838d0, 68.414136745272785d0, 73.626238908565369d0,& 78.841378309935209d0, 84.059017498643797d0, 89.278737291769232d0,& 94.500206298463695d0, 99.723159305319967d0, 104.94738168756264d0,& 110.17269799133665d0, 115.39896345131892d0, 120.62605760591623d0,& 125.85387943287243d0, 131.08234360156550d0, 136.31137755558950d0,& 141.54091921974700d0, 146.77091518163468d0, 152.00131923753867d0,& 157.23209122057966d0, 162.46319604943233d0, 167.69460295082322d0,& 172.92628481998585d0, 178.15821769142210d0, 183.39038029845793d0,& 188.62275370473256d0, 193.85532099431336d0, 199.08806700986247d0,& 204.32097813040141d0 /) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *,'NOTE: USED NAFEMS_CONVECTION_TEST_FLUX' END IF X = XYZ(1) ; Y = XYZ(2) DX_TEMP = 0.d0 ; DY_TEMP = 0.d0 DO I = 1, MAX_ROOT TEMP1 = COS(ROOT(I)*X) * (ROOT(I) * COSH(ROOT(I)*(B-Y)) & + H_BAR * SINH(ROOT(I)*(B-Y)) ) DX_TEMP1 = -SIN(ROOT(I)*X) * (ROOT(I) * COSH(ROOT(I)*(B-Y)) & + H_BAR * SINH(ROOT(I)*(B-Y)) ) * ROOT(I) DY_TEMP1 = -COS(ROOT(I)*X) * (ROOT(I) * SINH(ROOT(I)*(B-Y)) & + H_BAR * COSH(ROOT(I)*(B-Y)) ) * ROOT(I) TEMP2 = ( (ROOT(I)**2 + H_BAR_SQ) * A + H_BAR ) & * ( ROOT(I)*COSH(ROOT(I)*B) + H_BAR*SINH(ROOT(I)*B) ) & * COS(ROOT(I)*A) DX_TEMP = DX_TEMP + DX_TEMP1/TEMP2 DY_TEMP = DY_TEMP + DY_TEMP1/TEMP2 END DO ! TO APPROXIMATE INFINITE SERIES FLUX (1) = 2.d0 * H_BAR * T_AB * DX_TEMP FLUX (2) = 2.d0 * H_BAR * T_AB * DY_TEMP END SUBROUTINE NAFEMS_CONVECTION_TEST_FLUX SUBROUTINE NAFEMS_CONVECTION_TEST_MIXED (XYZ, S_m, C_m) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 6 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: S_m, C_m ! SQ & COL SOURCES INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: NAFEMS_CONVECTION_TEST_MIXED' END IF ! For a segment with the mixed boundary condition: ! k_x * u,x * n_x + k_y * u,y * n_y + S_m * u + C_m = 0 ! S_m = SQUARE MATRIX SOURCE = convection coeff, H_c ! C_m = COLUMN MATRIX SOURCE ! = - convection coeff * tempature = -H_c*T_c ! H = 750 W/mC, 0 C L_AB = 0.6 m, L_BC = 1.0 m ! D----------C AT (0.6, 0.2) T = 18.254 C ! q=0 | | H = 750 W/mC, 0 C ! | k=52W/mC | ! A----------B --> x Point B has very high gradients ! T_AB = 100 C S_m = 750.d0 ; C_m = 0.d0 ! for edges DC, CB IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE NAFEMS_CONVECTION_TEST_MIXED ! --------------------- EXACT_CASE = 7 ------------------------ SUBROUTINE POTENTIAL_FLOW_CYL_EXACT (XYZ, EXACT_SOL) ! Flow around a cylinder of radius A, located at (X_CL, 0) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 7 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), PARAMETER :: X_CL = 8.d0, UINF = 1.d0 REAL(DP), PARAMETER :: EPS = 0.00001d0 !b , PI = 3.141592654d0 REAL(DP), PARAMETER :: A = 2.d0, A_SQ = A * A REAL(DP) :: ATAN_ARG, R, THETA INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED POTENTIAL_FLOW_CYL_EXACT' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'POTENTIAL_FLOW_CYL_EXACT: INVALID CONSTANTS' END IF R = SQRT( (X_CL - XYZ(1))**2 + XYZ(2)**2 ) IF ( (X_CL - XYZ(1)) <= EPS ) THEN THETA = 90.0 ELSE ATAN_ARG = XYZ(2)/(X_CL - XYZ(1)) THETA = 180.0 - ATAN(ATAN_ARG)*(180.0/PI) ENDIF THETA = (PI/180.0)*THETA EXACT_SOL (1) = UINF*(R+(A_SQ/R))*COS(THETA) END SUBROUTINE POTENTIAL_FLOW_CYL_EXACT SUBROUTINE POTENTIAL_FLOW_CYL_EXACT_VEL (XYZ, FLUX) ! Flow around a cylinder of radius A, located at (X_CL, 0) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 7 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), PARAMETER :: X_CL = 8.d0, UINF = 1.d0 REAL(DP), PARAMETER :: EPS = 0.00001d0 !b , PI = 3.141592654d0 REAL(DP), PARAMETER :: A = 2.d0, A_SQ = A * A REAL(DP) :: ATAN_ARG, R_SQ, THETA, VR, VTHETA INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED ' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'POTENTIAL_FLOW_CYL_EXACT_VEL: INVALID CONSTANTS' END IF R_SQ = (X_CL - XYZ(1))**2 + XYZ(2)**2 IF ( (X_CL - XYZ(1)) <= EPS ) THEN THETA = 90.d0 ELSE ATAN_ARG = XYZ(2)/(X_CL - XYZ(1)) THETA = 180.d0 - ATAN(ATAN_ARG)*(180.d0/PI) ENDIF THETA = (PI/180.d0)*THETA VR = UINF*(1.d0 - (A_SQ/R_SQ))*COS(THETA) VTHETA = -UINF*(1.d0 + (A_SQ/R_SQ))*SIN(THETA) FLUX (1) = VR*COS(THETA) - VTHETA*SIN(THETA) FLUX (2) = VR*SIN(THETA) + VTHETA*COS(THETA) END SUBROUTINE POTENTIAL_FLOW_CYL_EXACT_VEL ! --------------------- EXACT_CASE = 8 ------------------------ SUBROUTINE PATCH_TEST_2D_PE_2_DER (XYZ, EXACT_SOL) ! SOLUTION VALUE FOR ! 2D POISSON EQ WITH CONSTANT SECOND DERIVATIVES Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 8 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_2D_PE_2_DER' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'PATCH_TEST_2D_PE_2_DER: INVALID CONSTANTS' END IF ! EXACT_SOL: u = x + 5y + 2x^2 + 3xy + 4y^2 ! const sec deriv EXACT_SOL (1) = XYZ(1) + 2.d0 * XYZ(1)**2 + 3.d0 * XYZ(1) * XYZ(2) & + 4.d0 * XYZ(2)**2 + 5.d0 * XYZ(2) END SUBROUTINE PATCH_TEST_2D_PE_2_DER SUBROUTINE PATCH_TEST_2D_PE_2_DER_SOURCE (XYZ, SOURCE) ! SOLUTION VALUE FOR ! 2D POISSON EQ WITH CONSTANT SECOND DERIVATIVES Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 8 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: PATCH_TEST_2D_PE_2_DER_SOURCE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'PATCH_TEST_2D_PE_2_DER_SOURCE: INVALID CONSTANTS' END IF ! EXACT_SOL: u = x + 5y + 2x^2 + 3xy + 4y^2 ! const sec deriv !-SOURCE = k_xx u,xx + 2k_xy U,xy + k_yy u,yy = 12 or 18, u,xy /= 0 SOURCE = -12.d0 ! used only if N_FLT_EL < 4 IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE PATCH_TEST_2D_PE_2_DER_SOURCE SUBROUTINE PATCH_TEST_2D_PE_2_DER_FLUX (XYZ, FLUX) ! TWO FLUX COMPONENTS (u,x AND u,y) FOR ! 2D POISSON EQ WITH CONSTANT SECOND DERIVATIVES Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 8 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), PARAMETER :: K_x = 1.d0, K_xy = 0.d0, K_y = 1.d0 ! Add more general K array later REAL(DP) :: U_x, U_y INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *,'NOTE: USED PATCH_TEST_2D_PE_2_DER_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /= 2 ) & STOP 'PATCH_TEST_2D_PE_2_DER_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL: u = x + 5y + 2x^2 + 3xy + 4y^2 ! const sec deriv ! u,x = 1 + 4 x + 3 y, u,y = 5 + 3 x + 8 y U_x = 1.d0 + 4.d0 * XYZ(1) + 3.d0 * XYZ(2) U_y = 5.d0 + 3.d0 * XYZ(1) + 8.d0 * XYZ(2) FLUX (1) = K_x * U_x + K_xy * U_y FLUX (2) = K_xy * U_x + K_y * U_y END SUBROUTINE PATCH_TEST_2D_PE_2_DER_FLUX SUBROUTINE PATCH_TEST_2D_PE_2_DER_FG (XYZ, FLUX_GRAD) ! GRADIENT OF EACH OF TWO FLUX COMPONENTS ! (u,xx, u,xy u,yx AND u,yy) FOR ! 2D POISSON EQ WITH CONSTANT SECOND DERIVATIVES Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 8 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX_GRAD (N_2_DER) ! PT EXACT VALUE REAL(DP), PARAMETER :: K_x = 1.d0, K_xy = 1.d0, K_y = 1.d0 REAL(DP), PARAMETER :: U_xx = 4.d0, U_xy = 3.d0, U_yy = 8.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_2D_PE_2_DER_FG' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /= 2 ) & STOP 'PATCH_TEST_2D_PE_2_DER_FG: INVALID CONSTANTS' END IF ! EXACT_SOL: u = x + 5y + 2x^2 + 3xy + 4y^2 ! const sec deriv ! 1:u,xx = 4 2:u,xy = 3 = u,yx 3:u,yy = 8 FLUX_GRAD (1) = 4.d0 FLUX_GRAD (2) = 3.d0 FLUX_GRAD (3) = 8.d0 IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE PATCH_TEST_2D_PE_2_DER_FG ! --------------------- EXACT_CASE = 9 ------------------------ SUBROUTINE UXX_U_X (XYZ, EXACT_SOL) ! U,XX + U + X = 0, U(0)=0=U(1), U = Sin(x)/Sin(1) - x Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 9 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), PARAMETER :: SIN_1 = 0.841470985d0 ! Sin (1) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_U_X' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 ) & STOP 'UXX_U_X: INVALID CONSTANTS' END IF EXACT_SOL (1) = SIN (XYZ(1)) / SIN_1 - XYZ(1) END SUBROUTINE UXX_U_X SUBROUTINE UXX_U_X_FLUX (XYZ, FLUX) ! U,XX + U + X = 0, U(0)=0=U(1), U = Sin(x)/Sin(1) - x Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 9 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), PARAMETER :: SIN_1 = 0.841470985d0 ! Sin (1) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_U_X_FLUX' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 .OR. N_R_B > 2 ) & STOP 'UXX_U_X_FLUX: INVALID CONSTANTS' END IF FLUX (1) = COS (XYZ(1)) / SIN_1 - 1.d0 END SUBROUTINE UXX_U_X_FLUX ! --------------------- EXACT_CASE = 10 ------------------------ SUBROUTINE UXX_U_X_NBC (XYZ, EXACT_SOL) ! U,XX + U + X = 0, U(0)=0,U'(1)=0, U = Sin(x)/Cos(1) - x Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 10 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), PARAMETER :: COS_1 = 0.540302306d0 ! Cos (1) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_U_X_NBC' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 ) & STOP 'UXX_U_X_NBC: INVALID CONSTANTS' END IF EXACT_SOL (1) = SIN (XYZ(1)) / COS_1 - XYZ(1) END SUBROUTINE UXX_U_X_NBC SUBROUTINE UXX_U_X_FLUX_NBC (XYZ, FLUX) ! U,XX + U + X = 0, U(0)=0,U'(1)=0, U = Sin(x)/Cos(1) - x Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 10 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), PARAMETER :: COS_1 = 0.540302306d0 ! Cos (1) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_U_X_FLUX_NBC' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 .OR. N_R_B > 2 ) & STOP 'UXX_U_X_FLUX_NBC: INVALID CONSTANTS' END IF FLUX (1) = COS (XYZ(1)) / COS_1 - 1.d0 END SUBROUTINE UXX_U_X_FLUX_NBC ! --------------------- EXACT_CASE = 11 ------------------------ SUBROUTINE UXX_X_TO_N (XYZ, EXACT_SOL) ! U,xx + X^N = 0, U(0)=0=U(1), U = (X-X^N)/((N+1)(N+2)) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 11 Use Sys_Properties_Data ! for GET_INTEGER_MISC function IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER :: N ! exponent data in miscellaneous data INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_X_TO_N' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 ) & STOP 'UXX_X_TO_N: INVALID CONSTANTS' END IF N = GET_INTEGER_MISC (1) ! first integer misc property EXACT_SOL (1) = (XYZ(1) - XYZ(1)**(N+2))/((N + 1)*(N + 2)) END SUBROUTINE UXX_X_TO_N SUBROUTINE UXX_X_TO_N_SOURCE (XYZ, SOURCE) ! U,xx + X^N = 0, U(0)=0=U(1), U = (X-X^N)/((N+1)(N+2)) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 11 Use Sys_Properties_Data ! for GET_INTEGER_MISC function IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE INTEGER :: N ! exponent data in miscellaneous data INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_X_TO_N_SOURCE' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 ) & STOP 'UXX_X_TO_N_SOURCE: INVALID CONSTANTS' END IF N = GET_INTEGER_MISC (1) ! first integer misc property SOURCE = XYZ(1)**N END SUBROUTINE UXX_X_TO_N_SOURCE SUBROUTINE UXX_X_TO_N_FLUX (XYZ, FLUX) ! U,xx + X^N = 0, U(0)=0=U(1), U = (X-X^N)/((N+1)(N+2)) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 11 Use Sys_Properties_Data ! for GET_INTEGER_MISC function IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER :: N ! exponent data in miscellaneous data INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_X_TO_N_FLUX' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 .OR. N_R_B > 2 ) & STOP 'UXX_X_TO_N_FLUX: INVALID CONSTANTS' END IF N = GET_INTEGER_MISC (1) ! first integer misc property FLUX (1) = (1.d0 - (N+2)*XYZ(1)**(N+1))/((N + 1)*(N + 2)) END SUBROUTINE UXX_X_TO_N_FLUX ! --------------------- EXACT_CASE = 12 ------------------------ SUBROUTINE PATCH_TEST_2D_STRESS (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 12 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_2D_STRESS' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 2 ) & STOP 'PATCH_TEST_2D_STRESS: INVALID CONSTANTS' END IF ! EXACT_SOL U = V = 1 + 3 * X - 4 * Y ! displacements EXACT_SOL (1) = 1.d0 + 3.d0 * XYZ(1) - 4.d0 * XYZ(2) EXACT_SOL (2) = 1.d0 + 3.d0 * XYZ(1) - 4.d0 * XYZ(2) END SUBROUTINE PATCH_TEST_2D_STRESS SUBROUTINE PATCH_TEST_2D_STRESS_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 12 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP) :: E (N_R_B, N_R_B) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_2D_STRESS_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 2 .OR. N_R_B /=3 ) & STOP 'PATCH_TEST_2D_STRESS_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL U = V = 1 + 3 * X - 4 * Y ! STRAINS : 3, -4, -1 FLUX (1) = 3.d0 ! dU/dX FLUX (2) = -4.d0 ! dV/dY FLUX (3) = -1.d0 ! dU/dY + dV/dX ! STRESSES IF ( PLANE_STRAIN ) THEN CALL E_PLANE_STRAIN (E) ELSE CALL E_PLANE_STRESS (E) END IF FLUX = MATMUL (E, FLUX) IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE PATCH_TEST_2D_STRESS_FLUX ! --------------------- EXACT_CASE = 13 ------------------------ SUBROUTINE PATCH_TEST_AXI_SYM_STRESS (XYZ, EXACT_SOL) ! AXISYMMETRIC SOLID WITH CONSTANT STRAINS (NEED E=1, NU=0) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 13 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_AXI_SYM_STRESS' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 2 ) & STOP 'PATCH_TEST_AXI_SYM_STRESS: INVALID CONSTANTS' END IF ! EXACT: U = R, V = 1 + 3 * R + 2 * Z EXACT_SOL (1) = XYZ(1) EXACT_SOL (2) = 1.d0 + 3.d0 * XYZ(1) + 2.d0 * XYZ(2) END SUBROUTINE PATCH_TEST_AXI_SYM_STRESS SUBROUTINE PATCH_TEST_AXI_SYM_STRESS_FLUX (XYZ, FLUX) ! AXISYMMETRIC SOLID WITH CONSTANT STRAINS (NEED E=1, NU=0) Use System_Constants ! for N_G_DOF, N_SPACEE, EXACT_CASE = 13 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: PATCH_TEST_AXI_SYM_STRESS_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=4 ) & STOP 'PATCH_TEST_AXI_SYM_STRESS_FLUX: INVALID CONSTANTS' END IF ! EXACT: e_rr = U,r = 1, e_zz = V,z = 2, ! e_rz = U,z + V,r = 0 + 3 = 3, e_tt = U/r = 1 FLUX (1) = 1.d0 ! e_rr FLUX (2) = 2.d0 ! e_zz FLUX (3) = 3.d0 ! e_rz FLUX (4) = 1.d0 ! e_tt IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE PATCH_TEST_AXI_SYM_STRESS_FLUX ! --------------------- EXACT_CASE = 14 ------------------------ SUBROUTINE UXX_Q_STEP (XYZ, EXACT_SOL) ! U,xx + Q_STEP = 0, U(0)=0=U(1), Q = 1 for X<=1/2 else = 0 ! U = X(3-4X)/8, X <= 1/2, else U = (1-X)/8 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 14 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_Q_STEP' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 ) & STOP 'UXX_Q_STEP: INVALID CONSTANTS' END IF IF ( XYZ(1) <= 0.5d0 ) THEN EXACT_SOL (1) = XYZ(1) * (3.d0 -4.d0 * XYZ(1)) / 8.d0 ELSE EXACT_SOL (1) = (1.d0 - XYZ(1)) / 8.d0 END IF END SUBROUTINE UXX_Q_STEP SUBROUTINE UXX_Q_STEP_SOURCE (XYZ, SOURCE) ! U,xx + Q_STEP = 0, U(0)=0=U(1), Q = 1 for X<=1/2 else = 0 ! U = X(3-4X)/8, X <= 1/2, else U = (1-X)/8 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 14 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_Q_STEP_SOURCE' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 ) & STOP 'UXX_Q_STEP_SOURCE: INVALID CONSTANTS' END IF IF ( XYZ(1) <= 0.5d0 ) THEN SOURCE = 1.d0 ELSE ; SOURCE = 0.d0 ; END IF END SUBROUTINE UXX_Q_STEP_SOURCE SUBROUTINE UXX_Q_STEP_FLUX (XYZ, FLUX) ! U,xx + Q_STEP = 0, U(0)=0=U(1), Q = 1 for X<=1/2 else = 0 ! U = X(3-4X)/8, X <= 1/2, else U = (1-X)/8 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 14 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_Q_STEP_FLUX' IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 .OR. N_R_B > 2 ) & STOP 'UXX_Q_STEP_FLUX: INVALID CONSTANTS' END IF IF ( XYZ(1) <= 0.5d0 ) THEN FLUX (1) = (3.d0 - XYZ(1)) / 8.d0 ELSE ; FLUX (1) = - 0.125d0 ; END IF END SUBROUTINE UXX_Q_STEP_FLUX ! --------------------- EXACT_CASE = 15 ------------------------ SUBROUTINE UXX_FAUSETT_EBC (XYZ, EXACT_SOL) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 482 ! with Y(0)=2, Y(1)=5/3, Y(X) = X^4/6 - 3X^2/2 + X + 2 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 15 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 LOGICAL, SAVE :: L_SQ ! Least Squares option used ? IF ( NOTE == 0 ) THEN ; NOTE = 1 IF ( N_SPACE <= 2 .AND. N_R_B <= 2 .AND. N_G_DOF == 1 ) THEN L_SQ = .FALSE. ! Galerkin method C_0 element IF ( DEBUG_EXACT ) PRINT *, 'NOTE: UXX_FAUSETT_EBC Galerkin' ELSE IF ( N_SPACE == 1 .AND. N_R_B == 2 .AND. N_G_DOF == 2 ) THEN L_SQ = .TRUE. ! Least Squares method C_1 element IF ( DEBUG_EXACT ) PRINT *, 'NOTE: UXX_FAUSETT_EBC Least Squares' ELSE STOP 'UXX_FAUSETT_EBC: INVALID CONSTANTS' END IF ! type of analysis END IF EXACT_SOL (1) = XYZ(1)**4 / 6.d0 - 1.5d0 * XYZ(1)**2 + XYZ(1) + 2.d0 IF ( L_SQ ) EXACT_SOL (2) = 4.d0 * XYZ(1)**3 / 6.d0 & - 3.d0 * XYZ(1) + 1.d0 END SUBROUTINE UXX_FAUSETT_EBC SUBROUTINE UXX_FAUSETT_EBC_SOURCE (XYZ, SOURCE) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 482 ! with Y(0)=2, Y(1)=5/3, Y(X) = X^4/6 - 3X^2/2 + X + 2 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 15 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_EBC_SOURCE' END IF SOURCE = -XYZ(1)**2 + 1.d0 END SUBROUTINE UXX_FAUSETT_EBC_SOURCE SUBROUTINE UXX_FAUSETT_EBC_FLUX (XYZ, FLUX) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 482 ! with Y(0)=2, Y(1)=5/3, Y(X) = X^4/6 - 3X^2/2 + X + 2 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 15 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 LOGICAL, SAVE :: L_SQ = .FALSE. ! Least Squares option used IF ( NOTE == 0 ) THEN ; NOTE = 1 IF ( N_G_DOF == 2 ) L_SQ = .TRUE. IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_EBC_FLUX' END IF FLUX (1) = 4.d0 * XYZ(1)**3 / 6.d0 - 3.d0 * XYZ(1) + 1.d0 IF ( L_SQ ) FLUX (2) = 2.d0 * XYZ(1)**2 - 3.d0 END SUBROUTINE UXX_FAUSETT_EBC_FLUX ! --------------------- EXACT_CASE = 16 ------------------------ SUBROUTINE UXX_FAUSETT_EBC_RBC (XYZ, EXACT_SOL) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 16 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 LOGICAL, SAVE :: L_SQ = .FALSE. ! Least Squares option used IF ( NOTE == 0 ) THEN ; NOTE = 1 IF ( N_SPACE <= 2 .AND. N_R_B <= 2 .AND. N_G_DOF == 1 ) THEN L_SQ = .FALSE. ! Galerkin method C_0 element IF ( DEBUG_EXACT ) PRINT *, 'NOTE: UXX_FAUSETT_EBC_RBC Galerkin' ELSE IF ( N_SPACE == 1 .AND. N_R_B == 2 .AND. N_G_DOF == 2 ) THEN L_SQ = .TRUE. ! Least Squares method C_1 element IF ( DEBUG_EXACT ) PRINT *, 'NOTE: UXX_FAUSETT_EBC_RBC Least Sq' ELSE STOP 'UXX_FAUSETT_EBC_RBC: INVALID CONSTANTS' END IF ! type of analysis END IF EXACT_SOL (1) = (XYZ(1)**4 - 3.d0 * XYZ(1)**2 - XYZ(1) + 6.d0)/6.d0 IF ( L_SQ ) EXACT_SOL (2) = (4.d0 * XYZ(1)**3 & - 6.d0 * XYZ(1) - 1.d0)/6.d0 END SUBROUTINE UXX_FAUSETT_EBC_RBC SUBROUTINE UXX_FAUSETT_EBC_RBC_SOURCE (XYZ, SOURCE) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 16 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_EBC_RBC_SOURCE' END IF SOURCE = -XYZ(1)**2 + 1.d0 END SUBROUTINE UXX_FAUSETT_EBC_RBC_SOURCE SUBROUTINE UXX_FAUSETT_EBC_RBC_FLUX (XYZ, FLUX) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 16 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 LOGICAL, SAVE :: L_SQ = .FALSE. ! Least Squares option used IF ( NOTE == 0 ) THEN ; NOTE = 1 IF ( N_G_DOF == 2 ) L_SQ = .TRUE. IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_EBC_RBC_FLUX' END IF FLUX (1) = (4.d0 * XYZ(1)**3 - 6.d0 * XYZ(1) - 1.d0)/6.d0 IF ( L_SQ ) FLUX (2) = 2.d0 * XYZ(1)**2 - 1.d0 END SUBROUTINE UXX_FAUSETT_EBC_RBC_FLUX SUBROUTINE UXX_FAUSETT_EBC_RBC_ROBIN (XYZ, ROBIN_1, ROBIN_2) ! Mixed or Robin boundary condition, Standard form: ! K_n * U,n + ROBIN_1_SEG * U + ROBIN_2_SEG = 0 ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 16 !b Use Select_Source ! for ROBIN_1_SEG, ROBIN_2_SEG IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: ROBIN_1, ROBIN_2 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_EBC_RBC_ROBIN' END IF ROBIN_2 = 0.d0 ROBIN_1 = 1.d0 ! since K_n * U,n = U,x at x=x_max IF ( XYZ (1) /= 1.d0 ) THEN PRINT *,'WARNING: UXX_FAUSETT_EBC_RBC_ROBIN, INVALID POINT' N_WARN = N_WARN + 1 END IF END SUBROUTINE UXX_FAUSETT_EBC_RBC_ROBIN ! --------------------- EXACT_CASE = 17 ------------------------ SUBROUTINE UXX_FAUSETT_RBC_RBC (XYZ, EXACT_SOL) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y'(0)+Y(0)=0, Y'(1)-Y(1)=3, Y=X^4/6 +3X^2/2 +X -1 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 17 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 LOGICAL, SAVE :: L_SQ ! Least Squares option used ? IF ( NOTE == 0 ) THEN ; NOTE = 1 IF ( N_SPACE <= 2 .AND. N_R_B <= 2 .AND. N_G_DOF == 1 ) THEN L_SQ = .FALSE. ! Galerkin method C_0 element IF ( DEBUG_EXACT ) PRINT *, 'NOTE: UXX_FAUSETT_RBC_RBC Galerkin' ELSE IF ( N_SPACE == 1 .AND. N_R_B == 2 .AND. N_G_DOF == 2 ) THEN L_SQ = .TRUE. ! Least Squares method C_1 element IF ( DEBUG_EXACT ) PRINT *, 'NOTE: UXX_FAUSETT_RBC_RBC Least Sq' ELSE STOP 'UXX_FAUSETT_EBC: INVALID CONSTANTS' END IF ! type of analysis END IF EXACT_SOL (1) = XYZ(1)**4 /6.d0 + 1.5d0 * XYZ(1)**2 + XYZ(1) - 1.d0 IF ( L_SQ ) EXACT_SOL (2) =4.d0 * XYZ(1)**3 /6.d0 +3.d0 * XYZ(1) +1.d0 END SUBROUTINE UXX_FAUSETT_RBC_RBC SUBROUTINE UXX_FAUSETT_RBC_RBC_SOURCE (XYZ, SOURCE) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y'(0)+Y(0)=0, Y'(1)-Y(1)=3, Y=X^4/6 +3X^2/2 +X -1 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 17 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_RBC_RBC_SOURCE' END IF SOURCE = -XYZ(1)**2 + 1.d0 END SUBROUTINE UXX_FAUSETT_RBC_RBC_SOURCE SUBROUTINE UXX_FAUSETT_RBC_RBC_FLUX (XYZ, FLUX) ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y'(0)+Y(0)=0, Y'(1)-Y(1)=3, Y=X^4/6 +3X^2/2 +X -1 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 17 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 LOGICAL, SAVE :: L_SQ = .FALSE. ! Least Squares option used IF ( NOTE == 0 ) THEN ; NOTE = 1 IF ( N_G_DOF == 2 ) L_SQ = .TRUE. IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_RBC_RBC_FLUX' END IF FLUX (1) = 4.d0 * XYZ(1)**3 /6.d0 + 3.d0 * XYZ(1) + 1.d0 IF ( L_SQ ) FLUX (2) = 2.d0 * XYZ(1)**2 + 3.d0 END SUBROUTINE UXX_FAUSETT_RBC_RBC_FLUX SUBROUTINE UXX_FAUSETT_RBC_RBC_ROBIN (XYZ, ROBIN_1, ROBIN_2) ! Mixed or Robin boundary condition, Standard form: ! K_n * U,n + ROBIN_1_SEG * U + ROBIN_2_SEG = 0 ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 484 ! with Y'(0)+Y(0)=0, Y'(1)-Y(1)=3, Y=X^4/6 +3X^2/2 +X -1 Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 17 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: ROBIN_1, ROBIN_2 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED UXX_FAUSETT_RBC_RBC_ROBIN' END IF IF ( XYZ (1) == 1.d0 ) THEN ROBIN_2 = -3.d0 ROBIN_1 = -1.d0 ! since K_n * U,n = U,x at x=x_max ELSE IF ( XYZ (1) == 0.d0 ) THEN ROBIN_2 = 0.d0 ROBIN_1 = -1.d0 ! since K_n * U,n = -U,x at x=x_min ELSE PRINT *,'WARNING: UXX_FAUSETT_RBC_RBC_ROBIN, INVALID POINT', XYZ N_WARN = N_WARN + 1 END IF END SUBROUTINE UXX_FAUSETT_RBC_RBC_ROBIN ! --------------------- EXACT_CASE = 18 ------------------------ SUBROUTINE SUPG_1D_01 (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 18 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! Pe * u,x + u,xx = 0, u(0) = 1, u(1) = 0 REAL(DP), SAVE :: Pe, E_Pe INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) PRINT *, & 'NOTE: USED SUPG_1D_01' ! Recover global Peclet number from properties IF ( EXACT_FL > 0 ) THEN Pe = GET_REAL_EXACT (1) ! Recover saved Pe E_Pe = EXP (Pe) ELSE IF ( MISC_FL > 0 ) THEN Pe = GET_REAL_MISC (1) ! Recover saved Pe E_Pe = EXP (Pe) ELSE PRINT *,'SUPG_1D_01: Exact Peclet number not given' PRINT *,'CHECK exact_reals IN keyword CONTROL' STOP 'SUPG_1D_01: Exact Peclet number not given' END IF ! DATA PRESENT END IF ! FIRST CALL EXACT_SOL (1) = (EXP(Pe * XYZ(1)) - E_Pe)/(1.d0 - E_Pe) ! allow for Hermite solution (Least Squares) IF (N_G_DOF > 1 ) EXACT_SOL (2) = Pe*EXP(Pe * XYZ(1)) / (1.d0 - E_Pe) END SUBROUTINE SUPG_1D_01 SUBROUTINE SUPG_1D_01_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE ! Pe * u,x + u,xx = 0, u(0) = 1, u(1) = 0 REAL(DP), SAVE :: Pe, E_Pe INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) PRINT *, & 'NOTE: USED SUPG_1D_01_FLUX' ! Recover global Peclet number from properties IF ( EXACT_FL > 0 ) THEN Pe = GET_REAL_EXACT (1) ! Recover saved Pe E_Pe = EXP (Pe) ELSE IF ( MISC_FL > 0 ) THEN Pe = GET_REAL_MISC (1) ! Recover saved Pe E_Pe = EXP (Pe) ELSE PRINT *,'SUPG_1D_01_FLUX: Exact Peclet number not given' PRINT *,'CHECK exact_reals IN keyword CONTROL' STOP 'SUPG_1D_01_FLUX: Exact Peclet number not given' END IF ! DATA PRESENT END IF ! FIRST CALL FLUX (1) = Pe * EXP(Pe * XYZ(1)) / (1.d0 - E_Pe) END SUBROUTINE SUPG_1D_01_FLUX SUBROUTINE SUPG_1D_01_FLUX_GRAD (XYZ, FLUX_GRAD) ! 2nd deriv ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET EXACT FLUX GRADIENTS (2ND DERIVATIVES) AT A NODE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for MAX_NP, N_R_B, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT(OUT) :: FLUX_GRAD (N_R_B * N_SPACE) ! FLUX_GRAD = EXACT FLUX GRAD (2ND DERIV) COMPONENTS AT A POINT ! HIGHER = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS ! Pe * u,x + u,xx = 0, u(0) = 1, u(1) = 0 REAL(DP), SAVE :: Pe, E_Pe INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) PRINT *, & 'NOTE: USED SUPG_1D_01_FLUX_GRAD' ! Recover global Peclet number from properties IF ( EXACT_FL > 0 ) THEN Pe = GET_REAL_EXACT (1) ! Recover saved Pe E_Pe = EXP (Pe) ELSE IF ( MISC_FL > 0 ) THEN Pe = GET_REAL_MISC (1) ! Recover saved Pe E_Pe = EXP (Pe) ELSE PRINT *,'SUPG_1D_01_FLUX_GRAD: Exact Peclet number not given' PRINT *,'CHECK exact_reals IN keyword CONTROL' STOP 'SUPG_1D_01_FLUX_GRAD: Exact Peclet number not given' END IF ! DATA PRESENT END IF ! FIRST CALL SETUP FLUX_GRAD (1) = Pe * Pe * EXP(Pe * XYZ(1)) / (1.d0 - E_Pe) END SUBROUTINE SUPG_1D_01_FLUX_GRAD ! --------------------- EXACT_CASE = 19 ------------------------ SUBROUTINE CARSLAW_ORTHOTROPIC (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 19 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! u,xx + u,yy / k_sq = - Q / K_xx, k_sq = K_xx / K_yy ! rectangle 2A by 2B in x- y-directions, axis at center INTEGER, PARAMETER :: N_SUM = 50 REAL(DP), SAVE :: K_XX, K_YY, K_XY, SOURCE, k_sq, k REAL(DP), SAVE :: A, A_SQ, B, C_1, C_2, C_3, C_4, C_5 ! sizes INTEGER, SAVE :: NOTE = 0 REAL(DP) :: RATIO, TERMS, COS_TERM INTEGER :: N, M IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY PRINT *, 'NOTE: USED CARSLAW_ORTHOTROPIC' ! Recover data (see example 202) IF ( EXACT_REALS > 5 ) THEN K_XX = GET_REAL_EXACT (1) ; K_YY = GET_REAL_EXACT (2) K_XY = GET_REAL_EXACT (3) ; SOURCE = GET_REAL_EXACT (4) A = GET_REAL_EXACT (5) ; B = GET_REAL_EXACT (6) ELSE IF ( N_LP_FLO > 3 ) THEN K_XX = GET_REAL_LP (1) ; K_YY = GET_REAL_LP (2) K_XY = GET_REAL_LP (3) ; SOURCE = GET_REAL_LP (4) IF ( MISC_FL < 2 ) STOP 'CARSLAW_ORTHOTROPIC: NO SIZES' A = GET_REAL_MISC (1) ; B = GET_REAL_MISC (2) ELSE STOP 'CARSLAW_ORTHOTROPIC: MISSING DATA' END IF ! DATA PRESENT IF (K_XY /= 0.d0 ) THEN PRINT *,'WARNING: CARSLAW_ORTHOTROPIC, K_XY NOT ZERO' N_WARN = N_WARN + 1 ; END IF ! orthotropic IF ( K_YY > 0.d0 ) THEN k_sq = K_xx / K_yy ; k = sqrt ( k_sq ) ELSE ; STOP 'CARSLAW_ORTHOTROPIC: IMPOSSIBLE K_YY' END IF ! valid data A_SQ = A*A C_1 = -SOURCE / ( 2 * K_xx) C_2 = C_1 * 32 * A * A / PI **3 C_3 = PI / (2 * A) ! PI / (2 * A) C_4 = C_3 * k ! k * PI / (2 * A) C_5 = C_4 * B ! B * k * PI / (2 * A) END IF ! FIRST CALL TERMS = 0.d0 DO N = 0, N_SUM M = 2 * N + 1 RATIO = ( COSH (M * C_4 * XYZ (2) ) / COSH (M * C_4 * B) ) COS_TERM = COS ( M * C_3 * XYZ (1)) IF ( N == (2 * ( N / 2 )) ) THEN ! EVEN TERMS = TERMS + COS_TERM * RATIO / M **3 ELSE ! ODD TERMS = TERMS - COS_TERM * RATIO / M **3 END IF END DO ! SUM EXACT_SOL (1) = C_1 * (XYZ (1) **2 - A_SQ) + C_2 * TERMS END SUBROUTINE CARSLAW_ORTHOTROPIC SUBROUTINE CARSLAW_ORTHOTROPIC_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 19 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE ! u,xx + u,yy / k_sq = - Q / K_xx, k_sq = K_xx / K_yy INTEGER, PARAMETER :: N_SUM = 50 REAL(DP), SAVE :: K_XX, K_YY, K_XY, SOURCE, k_sq, k REAL(DP), SAVE :: A, A_SQ, B, C_1, C_2, C_3, C_4, C_5 ! sizes REAL(DP), SAVE :: E (2, 2) ! E (N_R_B, N_R_B) ! constitutive INTEGER, SAVE :: NOTE = 0 REAL(DP) :: X_RATIO, Y_RATIO, X_TERM, Y_TERM, COS_TERM, SIN_TERM INTEGER :: N, M IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY PRINT *, 'NOTE: USED CARSLAW_ORTHOTROPIC FLUX' ! Recover data (see example 202) IF ( EXACT_REALS > 5 ) THEN K_XX = GET_REAL_EXACT (1) ; K_YY = GET_REAL_EXACT (2) K_XY = GET_REAL_EXACT (3) ; SOURCE = GET_REAL_EXACT (4) A = GET_REAL_EXACT (5) ; B = GET_REAL_EXACT (6) ELSE IF ( N_LP_FLO > 3 ) THEN K_XX = GET_REAL_LP (1) ; K_YY = GET_REAL_LP (2) K_XY = GET_REAL_LP (3) ; SOURCE = GET_REAL_LP (4) IF ( MISC_FL < 2 ) STOP 'CARSLAW_ORTHOTROPIC: NO SIZES' A = GET_REAL_MISC (1) ; B = GET_REAL_MISC (2) ELSE STOP 'CARSLAW_ORTHOTROPIC: MISSING DATA' END IF ! DATA PRESENT IF (K_XY /= 0.d0 ) THEN PRINT *,'WARNING: CARSLAW_ORTHOTROPIC, K_XY NOT ZERO' N_WARN = N_WARN + 1 ; END IF ! orthotropic IF ( K_YY > 0.d0 ) THEN k_sq = K_xx / K_yy ; k = sqrt ( k_sq ) ELSE ; STOP 'CARSLAW_ORTHOTROPIC: IMPOSSIBLE K_YY' END IF ! valid data A_SQ = A*A C_1 = -SOURCE / ( 2 * K_xx) C_2 = C_1 * 32 * A * A / PI **3 C_3 = PI / (2 * A) ! PI / (2 * A) C_4 = C_3 * k ! k * PI / (2 * A) C_5 = C_4 * B ! B * k * PI / (2 * A) CALL POISSON_ANISOTROPIC_2D_E_MATRIX (K_XX, K_YY, K_XY, E) END IF ! FIRST CALL X_TERM = 0.d0 ; Y_TERM = 0.d0 DO N = 0, N_SUM M = 2 * N + 1 X_RATIO = ( COSH (M * C_4 * XYZ (2) ) / COSH (M * C_4 * B) ) Y_RATIO = ( SINH (M * C_4 * XYZ (2) ) / COSH (M * C_4 * B) ) COS_TERM = COS ( M * C_3 * XYZ (1)) SIN_TERM = -SIN ( M * C_3 * XYZ (1)) IF ( N == (2 * ( N / 2 )) ) THEN ! EVEN X_TERM = X_TERM + C_3 * SIN_TERM * X_RATIO / M **2 Y_TERM = Y_TERM + C_4 * COS_TERM * Y_RATIO / M **2 ELSE ! ODD X_TERM = X_TERM - C_3 * SIN_TERM * X_RATIO / M **2 Y_TERM = Y_TERM - C_4 * COS_TERM * Y_RATIO / M **2 END IF END DO ! SUM FLUX = 0.d0 ! Actually gradient FLUX (1) = C_1 * XYZ (1) * 2 + C_2 * X_TERM FLUX (2) = C_2 * Y_TERM FLUX = MATMUL (E, FLUX) ! now flux END SUBROUTINE CARSLAW_ORTHOTROPIC_FLUX SUBROUTINE CARSLAW_ORTHOTROPIC_FLUX_GRAD (XYZ, FLUX_GRAD) ! 2nd deriv ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET EXACT FLUX GRADIENTS (2ND DERIVATIVES) AT A NODE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for MAX_NP, N_R_B, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT(OUT) :: FLUX_GRAD (N_R_B * N_SPACE) ! FLUX_GRAD = EXACT FLUX GRAD (2ND DERIV) COMPONENTS AT A POINT ! HIGHER = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS FLUX_GRAD = 0.d0 IF ( TOUCH ) print *, XYZ END SUBROUTINE CARSLAW_ORTHOTROPIC_FLUX_GRAD ! --------------------- EXACT_CASE = 20 ------------------------ SUBROUTINE CUBIC_2D_LAPLACE (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 20 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! u,xx + u,yy = 0, u=-x^3 -y^3 + 3x^2y + 3xy^2 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY PRINT *, 'NOTE: USED CUBIC_2D_LAPLACE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /= 2 ) & STOP 'CUBIC_2D_LAPLACE: INVALID DATA FOR THIS SOLUTION' END IF ! FIRST CALL EXACT_SOL (1) = - XYZ(1) **3 - XYZ(2) **3 & + 3 * XYZ (2) * XYZ (1) **2 & + 3 * XYZ (1) * XYZ (2) **2 END SUBROUTINE CUBIC_2D_LAPLACE SUBROUTINE CUBIC_2D_LAPLACE_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 19 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE ! u,xx + u,yy = 0, u=-x^3 -y^3 + 3x^2y + 3xy^2 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY PRINT *, 'NOTE: USED CUBIC_2D_LAPLACE_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /= 2 ) & STOP 'CUBIC_2D_LAPLACE_FLUX: INVALID DATA' END IF ! FIRST CALL FLUX (1) = -3 * XYZ(1) **2 + 6 * XYZ(1) * XYZ(2) + 3 * XYZ (2) **2 FLUX (2) = -3 * XYZ(2) **2 + 6 * XYZ(1) * XYZ(2) + 3 * XYZ (1) **2 END SUBROUTINE CUBIC_2D_LAPLACE_FLUX SUBROUTINE CUBIC_2D_LAPLACE_FLUX_GRAD (XYZ, FLUX_GRAD) ! 2nd deriv ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET EXACT FLUX GRADIENTS (2ND DERIVATIVES) AT A NODE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for MAX_NP, N_R_B, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT(OUT) :: FLUX_GRAD (N_R_B * N_SPACE) ! FLUX_GRAD = EXACT FLUX GRAD (2ND DERIV) COMPONENTS AT A POINT ! HIGHER = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS ! u,xx + u,yy = 0, u=-x^3 -y^3 + 3x^2y + 3xy^2 ! 1:u,xx 2:u,xy = u,yx 3:u,yy INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY PRINT *, 'NOTE: USED CUBIC_2D_LAPLACE_FLUX_GRAD' ; END IF FLUX_GRAD (1) = 6 * (-XYZ(1) + XYZ(2) ) ! u,xx FLUX_GRAD (2) = 6 * ( XYZ(1) + XYZ(2) ) ! u,xy = u,yx FLUX_GRAD (3) = 6 * ( XYZ(1) - XYZ(2) ) ! u,yy END SUBROUTINE CUBIC_2D_LAPLACE_FLUX_GRAD ! --------------------- EXACT_CASE = 21 ------------------------ SUBROUTINE CYLINDER_HEAT_EBC (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 21 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), PARAMETER :: R_in = 1.d0, R_out = 2.d0, T_in = 100.d0, & T_out = 10.d0, ratio = 0.6931471d0 ! Myers: Analytic Meth. in Heat Transfer, page 16 ! T = [T_in*ln(R_out/R) + T_out*ln(R/R_in)]/ln(R_out/R_in) ! T = T_in - (T_in - T_out)*ln(R/R_in)]/ln(R_out/R_in) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY PRINT *, 'NOTE: USED CYLINDER_HEAT_EBC PARAMETERS' !bIF ( N_SPACE /= 1 .OR. N_G_DOF /= 1 .OR. N_R_B /= 1 ) & ! allow general axisymmetric IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 .OR. N_R_B > 2 ) & STOP 'CYLINDER_HEAT_EBC: INVALID DATA FOR THIS SOLUTION' END IF ! FIRST CALL EXACT_SOL (1) = (T_in * LOG (R_out/xyz(1)) & + T_out * LOG (xyz(1)/R_in)) / ratio END SUBROUTINE CYLINDER_HEAT_EBC SUBROUTINE CYLINDER_HEAT_EBC_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE ! EXACT_CASE = 21 REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), PARAMETER :: R_in = 1.d0, R_out = 2.d0, T_in = 100.d0, & T_out = 10.d0, ratio = 0.6931471d0 ! T = T_in - (T_in - T_out)*ln(R/R_in)]/ln(R_out/R_in) ! F = TWO_PI*K*(T_in - T_out)/ln(R_out/R_in), total area ! FLUX = K*(T_in - T_out)/ln(R_out/R_in)/R, per unit area INTEGER, SAVE :: NOTE = 0, IE IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY !b IF ( N_SPACE /= 1 .OR. N_G_DOF /= 1 .OR. N_R_B /= 1 ) & IF ( N_SPACE > 2 .OR. N_G_DOF /= 1 .OR. N_R_B > 2 ) & STOP 'CYLINDER_HEAT_EBC_FLUX: INVALID DATA FOR THIS SOLUTION' END IF ! FIRST CALL IE = GET_THIS_ELEMENT_NUMBER () FLUX (1) = GET_REAL_LP (1)*(T_in - T_out) / ratio / XYZ(1) !b flux = -flux END SUBROUTINE CYLINDER_HEAT_EBC_FLUX ! --------------------- EXACT_CASE = 22 ------------------------ SUBROUTINE KREYSZIG_SPHERE_5 (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! Kreyszig 9.12.5 Unit sphere with EBC T(1, RADIANS) = cos (RADIANS)^2 REAL(DP), SAVE :: RADIANS, THIRD = 1.d0/3.d0 RADIANS = ATAN2 (XYZ(1), XYZ(2)) ! angle from Z EXACT_SOL = THIRD + (XYZ(1) **2 + XYZ(2) **2) & * ( COS(RADIANS) **2 - THIRD ) IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE KREYSZIG_SPHERE_5 SUBROUTINE KREYSZIG_SPHERE_5_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE ! include 'my_exact_flux_inc' ! Kreyszig 9.12.5 Unit sphere with EBC T(1, RADIANS) = cos (RADIANS)^2 REAL(DP), SAVE :: RADIANS, THIRD = 1.d0/3.d0, U_r, U_p, R RADIANS = ATAN2 (XYZ(1), XYZ(2)) ! angle from Z R = SQRT (XYZ(1) **2 + XYZ(2) **2) U_r = 2 * R * ( Cos(RADIANS) **2 - THIRD ) U_p = -2 * R * R * Cos(RADIANS) * Sin(RADIANS) FLUX (1) = U_r * Sin(RADIANS) + U_p * Cos(RADIANS) FLUX (2) = U_r * Cos(RADIANS) - U_p * Sin(RADIANS) IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE KREYSZIG_SPHERE_5_FLUX ! --------------------- EXACT_CASE = 23 ------------------------ SUBROUTINE HOLE_IN_INFINITE_PLATE (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP) :: ANG_X, R, Z, THICK REAL(DP) :: C1, C2, C3, C4, S1, S2, S3, S4 REAL(DP), SAVE :: T, Y_M, P_R, K, S_M, Cu, Kp, Km, K3 REAL(DP), SAVE :: A, HEIGHT = 1.d0 LOGICAL, SAVE :: FIRST_CALL = .TRUE. IF ( FIRST_CALL ) THEN ! one time calculations FIRST_CALL = .FALSE. IF ( N_SPACE /= 2 .OR. N_G_DOF /= 2 ) & STOP 'HOLE_IN_INFINITE_PLATE: INVALID CONSTANTS' IF ( EXACT_REALS > 3 ) THEN ! can use_exact_bc ! Exact properties: 1-radius of hole, A; 2-tension, T; ! 3-elastic modulus, Y_M; 4-Poisson's ratio, P_R, ! 5-height, else use 'el_reals' and 'reals' in keywords A = GET_REAL_EXACT (1) ; T = GET_REAL_EXACT (2) Y_M = GET_REAL_EXACT (3) ; P_R = GET_REAL_EXACT (4) IF ( EXACT_REALS > 4 ) HEIGHT = GET_REAL_EXACT (5) ELSE IF ( REALS > 1 ) THEN ! Misc properties: 1-radius of hole, A; 2-tension, T; ! 3-height A = GET_REAL_MISC (1) ; T = GET_REAL_MISC (2) IF ( REALS > 2 ) HEIGHT = GET_REAL_MISC (3) ! Element real properties: 1-elastic modulus, Y_M, ! 2-Poisson's ratio, P_R Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2) ELSE ! default geometry and data A = 1.d0; T=100.d0; Y_M = 1.d0; P_R=0.25d0 PRINT *, 'WARNING: EXACT HOLE_IN_INFINITE_PLATE DEFAULTED DATA TO' PRINT *, 'A = 1.d0; T=100.d0; Y_M = 1.d0; P_R=0.25d0' PRINT *, 'SET KEYWORD exact_reals 4 AND SUPPLY DATA AFTER REMARKS' N_WARN = N_WARN + 1 END IF ! property source S_M = 0.5d0 * Y_M / (1.d0 + P_R) ; THICK = 1.d0 IF ( PLANE_STRAIN ) THEN K = (3.d0 - P_R) /(1.d0 + P_R) ! plane strain ELSE K = 3.d0 - 4.d0 * P_R !plane Stress END IF Kp = 1.d0+K ; Km = 1.d0-K ; k3 = K - 3.d0 Cu = T * A / (8.d0 * S_M) END IF ! first call IF ( AREA_THICK /= 1.d0 ) THICK = AREA_THICK ! misc terms Z = SQRT ( XYZ (1) **2 + XYZ (2) **2 ) ; R = Z / A ! > 1 ANG_X = ATAN2 (XYZ(2), XYZ(1)) C1=Cos(ANG_X) ; C2=Cos(ANG_X*2); C3=Cos(ANG_X*3); C4=Cos(ANG_X*4) S1=Sin(ANG_X) ; S2=Sin(ANG_X*2); S3=Sin(ANG_X*3); S4=Sin(ANG_X*4) EXACT_SOL (1) = Cu*(R*Kp*C1 + 2.d0*(Kp*C1+C3)/r - 2.d0*C3/r/r/r) EXACT_SOL (2) = Cu*(R*K3*S1 + 2.d0*(Km*S1+S3)/r - 2.d0*S3/r/r/r) IF ( XYZ (1) == 0.d0 ) EXACT_SOL (1) = 0.d0 IF ( XYZ (2) == 0.d0 ) EXACT_SOL (2) = 0.d0 ! Scale to insure uniform traction over height EXACT_SOL = EXACT_SOL * 1.137217d0* HEIGHT / 3.d0 END SUBROUTINE HOLE_IN_INFINITE_PLATE SUBROUTINE HOLE_IN_INFINITE_PLATE_FLUX (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT FLUX FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP) :: ANG_T, ANG_X, R, R2, Z, THICK REAL(DP) :: C1, C2, C3, C4, S1, S2, S3, S4 REAL(DP) :: MOHR_2D_STRESS (3,3) REAL(DP), SAVE :: A, T !b , Y_M, P_R, K, S_M, Cu, Kp, Km, K3 LOGICAL, SAVE :: FIRST_CALL = .TRUE. IF ( FIRST_CALL ) THEN ! one time calculations FIRST_CALL = .FALSE. ! 3-elastic modulus, Y_M; 4-Poisson's ratio, P_R, ! 5-height, else use 'el_reals' and 'reals' in keywords IF ( N_SPACE /= 2 .OR. N_G_DOF /= 2 ) & STOP 'HOLE_IN_INFINITE_PLATE: INVALID CONSTANTS' IF ( EXACT_REALS > 3 ) THEN ! can use_exact_bc ! Exact properties: 1-radius of hole, A; 2-tension, T; ! 3-elastic modulus, S_M; 4-Poisson's ratio, P_R, ! 5-height, else use 'el_reals' and 'reals' in keywords A = GET_REAL_EXACT (1) ; T = GET_REAL_EXACT (2) ELSE IF ( REALS > 1 ) THEN ! Misc properties: 1-radius of hole, A; 2-tension, T; ! 3-height A = GET_REAL_MISC (1) ; T = GET_REAL_MISC (2) ELSE A = 1.d0; T=100.d0 END IF ! property source END IF ! first call THICK = 1.d0 IF ( AREA_THICK /= 1.d0 ) THICK = AREA_THICK ! misc terms Z = SQRT ( XYZ (1) **2 + XYZ (2) **2 ) R = Z / A ; R2 = R*R ! > 1 ANG_X = ATAN2 (XYZ(2), XYZ(1)) C1=Cos(ANG_X) ; C2=Cos(ANG_X*2); C3=Cos(ANG_X*3); C4=Cos(ANG_X*4) S1=Sin(ANG_X) ; S2=Sin(ANG_X*2); S3=Sin(ANG_X*3); S4=Sin(ANG_X*4) FLUX (1) = T*(1 - (1.5d0*C2+C4)/r/r + 1.5d0*C4/r/r/r/r) FLUX (2) = -T*((0.5d0*C2-C4)/r/r + 1.5d0*C4/r/r/r/r) FLUX (3) = -T*((0.5d0*S2+S4)/r/r - 1.5d0*S4/r/r/r/r) FLUX = FLUX / THICK END SUBROUTINE HOLE_IN_INFINITE_PLATE_FLUX ! ====================== ! This function provides the solution of the 2-D Steady state heat ! transfer, Babuska problem (ref: FEAD , Vol. 3, pp. 341-54, 1987) ! d = sqrt(r*R + z*Z) ! theta = atan2(Z,R) ! fofxy = 0.0700754*sqrt(d)* sin(theta/2.0) ! fofxy = 3.*R - 4.*Z**2 ! fofxy = 3.*R - 4.*Z ! strong diagonal ! p = (R+Z)/sqrt(2.0) ! fofxy=R*(1.-R)*Z*(1.-Z)*ATAN(20.*(p-0.8)) ! temp1 = 4.0 * R + 4.0 * Z ! temp2 = (EXP (temp1) - EXP(-temp1)) * 0.5d0 ! fofxy = 1.0 + R*Z + temp2 ! ====================== ! Here is the program to compute exact solution of eq.(5.1)-(5.3) in ! Codina, but the equation(5.2) is wrong. It should be exp(u/k) ! instead of exp(-u/k). ! PROGRAM MAIN !! Compute the exact solution of equ(5.1)-(5.3) in codina paper !b ! DOUBLE PRECISION, PARAMETER :: PI= 3.141592653589793d0 ! DOUBLE PRECISION, ALLOCATABLE :: X (:),F(:) ! DOUBLE PRECISION :: U, K, C1, C2 ! ! OPEN (11, FILE="input.dat") ! OPEN (12, FILE="output.dat") ! READ (11,*) N, U, K ! ALLOCATE (X(N), F(N)) ! ! X= (/((I-1.)/(N-1.), I=1,N)/) ! C2= 2*U/((U**2*PI+K**2*PI**3)*(1-EXP(U/K))) ! C1= -(1+EXP(U/K))*C2/2 ! F=C1+C2*EXP(U/K*X)+ K/(U**2+(K*PI)**2)*(SIN(PI*X)-U/(K*PI)*COS(PI*X)) ! WRITE (12,900) ! 900 FORMAT (6X,'X', 6X,'PHI') ! WRITE (12, 1000) (X(I), F(I), I=1,N) ! 1000 FORMAT (2(F10.5)) ! CLOSE (11) ! CLOSE (12) ! END PROGRAM MAIN ! !! ============input.dat !! 21 1 0.01 !! N U K ! copyright 1999, J. E. Akin, all rights reserved. ! --------------------- EXACT_CASE = 24 ------------------------ SUBROUTINE PATCH_TEST_2D_PE (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 24 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_2D_PE' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 ) & STOP 'PATCH_TEST_2D_PE: INVALID CONSTANTS' END IF ! EXACT_SOL U = 1 + 3 * X - 4 * Y ! displacements EXACT_SOL (1) = 1.d0 + 3.d0 * XYZ(1) - 4.d0 * XYZ(2) END SUBROUTINE PATCH_TEST_2D_PE SUBROUTINE PATCH_TEST_2D_PE_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 24 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP) :: E (N_R_B, N_R_B) INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED PATCH_TEST_2D_PE_FLUX' IF ( N_SPACE /= 2 .OR. N_G_DOF /= 1 .OR. N_R_B /=2 ) & STOP 'PATCH_TEST_2D_PE_FLUX: INVALID CONSTANTS' END IF ! EXACT_SOL U = V = 1 + 3 * X - 4 * Y ! STRAINS : 3, -4, -1 FLUX (1) = 3.d0 ! dU/dX FLUX (2) = -4.d0 ! dV/dY IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE PATCH_TEST_2D_PE_FLUX ! --------------------- EXACT_CASE = 25 ------------------------ SUBROUTINE COSINE_HILL_BC_ONLY (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! Define any new array or variable types, then give statements REAL(DP) :: R INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) PRINT *, & 'NOTE: USED COSINE_HILL_BC_ONLY' IF ( .NOT. USE_EXACT_BC ) THEN PRINT *, 'WARNING, COSINE_HILL_BC_ONLY VALID for BC ONLY,' PRINT *, 'PUT KEYWORD use_exact_bc IN CONTROL AND' PRINT *, 'REMOVE KEYWORDS list_exact AND list_exact_flux' N_WARN = N_WARN + 1 END IF END IF ! Codina CMAME 94,239-262, 1992, Ex 5, Cosine Hill interior BC R = SQRT (XYZ(1) **2 + XYZ(2) **2) IF ( R < 0.5d0 ) THEN EXACT_SOL (1) = SIN ( 0.5d0*TWO_PI*(1-2*R) ) ELSE EXACT_SOL (1) = 0.d0 END IF ! position END SUBROUTINE COSINE_HILL_BC_ONLY ! --------------------- EXACT_CASE = 26 ------------------------ SUBROUTINE HEINRICH_EXAMPLE_5_9 (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 26 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), PARAMETER :: DENOM = 6.3890560988d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED HEINRICH_EXAMPLE_5_9' IF ( N_R_B /= 1 .OR. N_G_DOF /= 1 ) & STOP 'HEINRICH_EXAMPLE_5_9: INVALID CONSTANTS' END IF ! - d(k * dp/dx)/dx + r*p = Q, k=1, r=1, Q=1, u(0) = 0, u(1) = 1 ! EXACT_SOL U = 1 + (EXP(X) - EXP (2-X))/(EXP (2) - 1)) EXACT_SOL (1) = 1 + (EXP(XYZ(1)) - EXP (2.d0-XYZ(1)))/DENOM END SUBROUTINE HEINRICH_EXAMPLE_5_9 SUBROUTINE HEINRICH_EXAMPLE_5_9_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 26 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE !b REAL(DP) :: E (N_R_B, N_R_B) REAL(DP), PARAMETER :: DENOM = 6.3890560988d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED HEINRICH_EXAMPLE_5_9_FLUX' IF ( N_G_DOF /= 1 .OR. N_R_B /= 1 ) & STOP 'HEINRICH_EXAMPLE_5_9_FLUX: INVALID CONSTANTS' END IF ! EXACT FLUX dU/dX = (EXP(X) + EXP (2-X))/(EXP (2) - 1) FLUX (1) = (EXP(XYZ(1)) + EXP (2.d0-XYZ(1)))/DENOM IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE HEINRICH_EXAMPLE_5_9_FLUX ! --------------------- EXACT_CASE = 27 ------------------------ SUBROUTINE ADE_1_D_TEST (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 27 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), PARAMETER :: EME3 = -17.367255095d0, DENOM = 3*EME3 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED ADE_1_D_TEST' IF ( N_R_B /= 1 .OR. N_G_DOF /= 1 ) & STOP 'ADE_1_D_TEST: INVALID CONSTANTS' END IF ! u * dp/dx - d(k * dp/dx)/dx + r*p = Q, assume u, k, r, Q constant ! 4 dp/dx - d(1 * dp/dx)/dx - 3 p = -1, so u=4, k=1, r=-3, Q=-1 ! for p(0)=1/3, p(1)=1 the exact solution is ! p(x) = [2(EXP(x) -EXP(3X) + (EXP - EXP(3))] / 3 * (EXP - EXP(3) EXACT_SOL (1) = (2*(EXP(XYZ(1)) - EXP(3*XYZ(1)))+ EME3)/DENOM END SUBROUTINE ADE_1_D_TEST SUBROUTINE ADE_1_D_TEST_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 27 IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE !b REAL(DP) :: E (N_R_B, N_R_B) REAL(DP), PARAMETER :: EME3 = -17.367255095d0, DENOM = 3*EME3 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED ADE_1_D_TEST_FLUX' IF ( N_G_DOF /= 1 .OR. N_R_B /= 1 ) & STOP 'ADE_1_D_TEST_FLUX: INVALID CONSTANTS' END IF ! EXACT FLUX dU/dX = (EXP(X) + EXP (2-X))/(EXP (2) - 1) FLUX (1) = 2*(EXP(XYZ(1)) - 3*EXP(3*XYZ(1))) / DENOM IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE ADE_1_D_TEST_FLUX ! --------------------- EXACT_CASE = 28 ------------------------ SUBROUTINE GRIFFITHS_TEST_1 (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 28 Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! u * dp/dx - d(k * dp/dx)/dx + r*p = Q, assume u, k, r, Q known ! u dp/dx - d(1 * dp/dx)/dx = 3u*x^2 so k=1, r=0, Q=3u*x^2 ! for p(0)=0, p(1)=0 the exact solution is p(x) = x^3 + 3x/u + ! + 6x/u/u -(1+3/u+6/u/u)*(EXP(ux)-1)/(EXP(u)-1) ! Here u=60 so scale = 60 and flow_case = 0 REAL(DP), SAVE :: X, U = 60.d0, EU_1 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED GRIFFITHS_TEST_1' ! IF ( N_R_B /= 1 .OR. N_G_DOF /= 1 ) & ! STOP 'GRIFFITHS_TEST_1: INVALID CONSTANTS' IF ( INTEGERS > 1 ) U = GET_INTEGER_MISC (2) EU_1 = EXP (U) - 1.d0 END IF X = XYZ (1) ! Griffiths Lorenz Example 11.1 EXACT_SOL (1) = X **3 + 3 * X **2 / U + 6 * X / U / U & - (1 + 3/U + 6/U/U)*(EXP(U*X) - 1) / EU_1 END SUBROUTINE GRIFFITHS_TEST_1 SUBROUTINE GRIFFITHS_TEST_1_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 28 Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), SAVE :: X, U = 60.d0, EU_1 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED GRIFFITHS_TEST_1_FLUX' ! IF ( N_G_DOF /= 1 .OR. N_R_B /= 1 ) & ! STOP 'GRIFFITHS_TEST_1_FLUX: INVALID CONSTANTS' IF ( INTEGERS > 1 ) U = GET_INTEGER_MISC (2) EU_1 = EXP (U) - 1.d0 END IF ! EXACT FLUX X = XYZ (1) ! Griffiths Lorenz Example 11.1 FLUX (1) = 3 * X **2 + 6 * X / U + 6 / U / U & - U * (1 + 3/U + 6/U/U)*(EXP(U*X) - 1) / EU_1 END SUBROUTINE GRIFFITHS_TEST_1_FLUX SUBROUTINE GRIFFITHS_TEST_1_SOURCE (XYZ, SOURCE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOURCE FOR ELEMENT MATRICES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE REAL(DP), SAVE :: X, U = 60.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED GRIFFITHS_TEST_1_SOURCE' ! IF ( N_G_DOF /= 1 .OR. N_R_B /= 1 ) & ! STOP 'GRIFFITHS_TEST_1_SOURCE: INVALID CONSTANTS' IF ( INTEGERS > 1 ) U = GET_INTEGER_MISC (2) END IF ! EXACT SOURCE X = XYZ (1) IF ( INTEGERS > 1 ) THEN ! Griffiths Lorenz Example 11.1 SOURCE = 3 * U * X **2 ELSE STOP 'Griffiths Lorenz Example 11.1 needs integers 2 for u=60' END IF END SUBROUTINE GRIFFITHS_TEST_1_SOURCE ! --------------------- EXACT_CASE = 29 ------------------------ SUBROUTINE GRIFFITHS_TEST_2 (XYZ, EXACT_SOL) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 29 Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! u * dp/dx - d(k * dp/dx)/dx + r*p = Q, assume u, k, r, Q known ! u dp/dx - d(1 * dp/dx)/dx = Q so k=1, r=0, with ! Q = U*(3*X **2 + TWO_PI*COS(TWO_PI*X)) + TWO_PI*TWO_PI*SIN(TWO_PI*X) ! for p(0)=0, p(1)=0 the exact solution is p(x) = x^3 + 3x/u + ! + 6x/u/u -(1+3/u+6/u/u)*(EXP(ux)-1)/(EXP(u)-1) + SIN (TWO_PI * X) ! Here u=60 so scale = 60 and flow_case = 0 REAL(DP), SAVE :: X, U = 60.d0, EU_1 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED GRIFFITHS_TEST_2' ! IF ( N_R_B /= 1 .OR. N_G_DOF /= 1 ) & ! STOP 'GRIFFITHS_TEST_2: INVALID CONSTANTS' IF ( INTEGERS > 1 ) U = GET_INTEGER_MISC (2) EU_1 = EXP (U) - 1.d0 END IF X = XYZ (1) ! Griffiths Lorenz Example 11.2 EXACT_SOL (1) = X **3 + 3 * X **2 / U + 6 * X / U / U & - (1 + 3/U + 6/U/U)*(EXP(U*X) - 1) / EU_1 & + SIN (TWO_PI * X) END SUBROUTINE GRIFFITHS_TEST_2 SUBROUTINE GRIFFITHS_TEST_2_FLUX (XYZ, FLUX) Use System_Constants ! for N_G_DOF, N_SPACE, EXACT_CASE = 29 Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT VALUE REAL(DP), SAVE :: X, U = 60.d0, EU_1 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED GRIFFITHS_TEST_2_FLUX' ! IF ( N_G_DOF /= 1 .OR. N_R_B /= 1 ) & ! STOP 'GRIFFITHS_TEST_2_FLUX: INVALID CONSTANTS' IF ( INTEGERS > 1 ) U = GET_INTEGER_MISC (2) EU_1 = EXP (U) - 1.d0 END IF ! EXACT FLUX X = XYZ (1) ! Griffiths Lorenz Example 11.2 FLUX (1) = 3 * X **2 + 6 * X / U + 6 / U / U & - U * (1 + 3/U + 6/U/U)*(EXP(U*X) - 1) / EU_1 & + TWO_PI * COS (TWO_PI * X) END SUBROUTINE GRIFFITHS_TEST_2_FLUX SUBROUTINE GRIFFITHS_TEST_2_SOURCE (XYZ, SOURCE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOURCE FOR ELEMENT MATRICES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, N_R_B Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: SOURCE ! PT EXACT VALUE REAL(DP), SAVE :: X, U = 60.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN NOTE = 1; IF ( DEBUG_EXACT ) PRINT *, 'NOTE: USED GRIFFITHS_TEST_2_SOURCE' ! IF ( N_G_DOF /= 1 .OR. N_R_B /= 1 ) & ! STOP 'GRIFFITHS_TEST_2_SOURCE: INVALID CONSTANTS' IF ( INTEGERS > 1 ) U = GET_INTEGER_MISC (2) END IF ! EXACT SOURCE X = XYZ (1) IF ( INTEGERS > 1 ) THEN ! Griffiths Lorenz Example 11.2 SOURCE = U * (3 * X **2 & + TWO_PI * COS (TWO_PI*X) ) & + TWO_PI * TWO_PI * SIN (TWO_PI*X) ELSE STOP 'Griffiths Lorenz Example 11.2 needs integers 2 for u=60' END IF END SUBROUTINE GRIFFITHS_TEST_2_SOURCE ! --------------------- EXACT_CASE = 30 B ------------------------ SUBROUTINE STOCASTIC_VOLATILITY_BC (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR BOUNDARY CONDITIONS ONLY, EXACT_CASE = 30 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! Donald Williams, for boundary conditions only REAL(DP), SAVE :: Max_S = 80.d0, Max_V = 80.d0 REAL(DP), SAVE :: Strike = 40.d0 INTEGER, SAVE :: NOTE = 0 IF ( NOTE == 0 ) THEN ; NOTE = 1; ! FIRST CALL ONLY IF ( EXACT_REALS >= 1 ) Strike = GET_REAL_EXACT (1) IF ( EXACT_REALS >= 3 ) THEN Max_S = GET_REAL_EXACT (2) Max_V = GET_REAL_EXACT (3) ELSE PRINT *,'DEFAULTED EXACT REAL PROPERTIES TO:' PRINT *,'Max_S, Max_V, Strike ',Max_S, Max_V, Strike END IF IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) PRINT *, & 'NOTE: USED 213.my_exact_inc', Strike, Max_S, Max_V END IF ! This version is correct for BC only ! IF ( XYZ (1) == 0.d0 ) EXACT_SOL (1) = Strike ! IF ( XYZ (1) > 0.d0 ) THEN ! EXACT_SOL (1) = MAX ((Strike - XYZ(1)), 0.d0) ! ELSE ! EXACT_SOL (1) = 0.d0 ! END IF ! position EXACT_SOL (1) = MAX ((Strike - XYZ(1)), 0.d0) IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) WRITE (N_BUG, *) & '213.my_exact_inc ', XYZ, EXACT_SOL (1) ! end file my_exact_inc END SUBROUTINE STOCASTIC_VOLATILITY_BC ! --------------------- EXACT_CASE = 31 ------------------------ SUBROUTINE EXACT_TRANSIENT_TEST_1 (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! SOLUTION OF Y' + A * Y = F ! y(x) = (1 - e (-Ax)) * F / A, with y(0) = 0 REAL(DP), SAVE :: A = 2.d0, F = 10.d0 LOGICAL, SAVE :: FIRST = .TRUE. IF ( FIRST ) THEN ! Get coefficients FIRST = .FALSE. IF ( REALS > 1 ) THEN A = GET_REAL_MISC (1) ; F = GET_REAL_MISC (2) ELSE ! Default to book example PRINT *,'NOTE: Default to book example A = 2, F = 10' END IF ! data given END IF ! First call EXACT_SOL (1) = (1 - EXP (-A * XYZ(1)) ) * F / A END SUBROUTINE EXACT_TRANSIENT_TEST_1 ! --------------------- EXACT_CASE = 32 ------------------------ SUBROUTINE CONVECT_DIFFUS_ON_SQ (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! Convection Diffusion of Passive Scalar on a Square ! Ref: Kobayashi, M. H., et al, J. Comp. Physics, 150, 40-75, 1999 ! Advection-Diffusion Equations: 1-D, 2-D, 3-D, Axisymmetric ! u * Del P - Del ( E Del P) + r P - Q = 0 ! u_1 = Sin (PI_BY_2L * XYZ (1)) * Cos (PI_BY_2L * XYZ (2)) ! u_2 = -Cos (PI_BY_2L * XYZ (1)) * Sin (PI_BY_2L * XYZ (2)) REAL(DP), PARAMETER :: LEN = 1.d0, PI_BY_2L = TWO_PI / (4 * LEN) EXACT_SOL (1) = SIN (PI_BY_2L * XYZ (1)) * SIN (PI_BY_2L * XYZ (2)) END SUBROUTINE CONVECT_DIFFUS_ON_SQ ! --------------------- EXACT_CASE = 33 ------------------------ SUBROUTINE RECT_BAR_BEND_DEFLECTIONS (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions Use Geometric_Properties ! for MAX_ MIN_XYZ (N_SPACE) IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE REAL(DP), SAVE :: A, B, L, I_Y, R ! GEOMETRY REAL(DP), SAVE :: M = 1.d0, Y_M = 1.d0, P_R = 0.25d0 ! DEFAULTS LOGICAL, SAVE :: FIRST = .TRUE. ! Pure Bending of a Rectangular Prismatic Solid Bar ! Ref: Timoshenko & Goodier "Theory of Elasticity" 2nd Ed, Sect 88 ! 3-D Elasticity. X vertically down, Y horizontal netural axis ! Z- horizontal center length line. Width 2B, Height 2A, Length L ! Sigma_Z = Ex/R all other stresses zero IF ( FIRST ) THEN ! Get exact constant values FIRST = .FALSE. A = ( MAX_XYZ (1) - MIN_XYZ (1) )/2 !b 0.5 B = ( MAX_XYZ (2) - MIN_XYZ (2) )/2 !b 2.0 L = MAX_XYZ (3) - MIN_XYZ (3) !b 8.0 I_Y = B * A **3 * 16.d0 / 12 ! Get the applied moment about y IF ( EXACT_REALS > 2 ) THEN Y_M = GET_REAL_EXACT (1) ; P_R = GET_REAL_EXACT (2) M = GET_REAL_EXACT (3) ELSE PRINT *, 'WARNING: DEFAULTED EXACT PROPERTIES, USE' PRINT *, 'exact_reals 3 FOR E, Nu, Moment' N_WARN = N_WARN + 1 END IF ! exact data given R = Y_M * I_Y / M END IF ! first call EXACT_SOL (1) = -0.5d0*( XYZ(3) **2 & + P_R*(XYZ(1) **2 - XYZ(2) **2) ) / R EXACT_SOL (2) = P_R * XYZ (1) * XYZ (2) / R EXACT_SOL (3) = P_R * XYZ (1) * XYZ (3) / R END SUBROUTINE RECT_BAR_BEND_DEFLECTIONS SUBROUTINE RECT_BAR_BEND_STRESSES (XYZ, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions Use Geometric_Properties ! for MAX_ MIN_XYZ (N_SPACE) IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: FLUX (N_R_B) ! PT EXACT STRESS REAL(DP), SAVE :: A, B, L, I_Y, R ! GEOMETRY REAL(DP), SAVE :: M = 1.d0, Y_M = 1.d0, P_R = 0.25d0 ! DEFAULTS LOGICAL, SAVE :: FIRST = .TRUE. ! Pure Bending of a Rectangular Prismatic Solid Bar ! Ref: Timoshenko & Goodier "Theory of Elasticity" 2nd Ed, Sect 88 ! 3-D Elasticity. X vertically down, Y horizontal netural axis ! Z- horizontal center length line. Width 2B, Height 2A, Length L ! Sigma_Z = Ex/R all other stresses zero IF ( FIRST ) THEN ! Get exact constant values FIRST = .FALSE. A = ( MAX_XYZ (1) - MIN_XYZ (1) )/2 !b 0.5 B = ( MAX_XYZ (2) - MIN_XYZ (2) )/2 !b 2.0 L = MAX_XYZ (3) - MIN_XYZ (3) !b 8.0 I_Y = B * A **3 * 16.d0 / 12 ! Get the applied moment about y IF ( EXACT_REALS > 2 ) THEN Y_M = GET_REAL_EXACT (1) ; P_R = GET_REAL_EXACT (2) M = GET_REAL_EXACT (3) END IF ! exact data given R = Y_M * I_Y / M END IF ! first call ! STRESS COMPONENT ORDER: XX, YY, XY, ZZ, XZ, YZ FLUX (1:3) = 0.d0 ; FLUX (4) = Y_M * XYZ (1) / R FLUX (5:6) = 0.d0 ! Exact strains are: -P_R*X/R, -P_R*X/R, 0, X/R, 0, 0 END SUBROUTINE RECT_BAR_BEND_STRESSES ! --------------------- EXACT_CASE = 34 ------------------------ SUBROUTINE TRANSIENT_HALF_WALL (XYZ, EXACT_SOL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT SOLUTION FOR OUTPUT OR ERROR ESTIMATES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_G_DOF, N_SPACE, MAX_XYZ Use Sys_Properties_Data ! for GET_INTEGER_* or GET_REAL_* functions IMPLICIT NONE REAL(DP), INTENT (IN) :: XYZ (N_SPACE) ! POINT IN SPACE REAL(DP), INTENT (OUT) :: EXACT_SOL (N_G_DOF) ! PT EXACT VALUE ! begin 122_01.my_exact_inc ! Myer's half plane wall transient Eq (3.1.7) page 84 INTEGER, SAVE :: N_SUM = 100 LOGICAL, SAVE :: FIRST = .TRUE. REAL(DP), SAVE :: L, A REAL(DP), SAVE :: K_e = 1.d0, Rho_e = 1.d0 ! default INTEGER :: M, N REAL(DP) :: INC, TERM_X, TERM_T, TOTAL REAL(DP) :: PI_SQ_T, PI_X_BY_L REAL(DP), PARAMETER :: TOLERANCE = 1.d-7 IF ( FIRST ) THEN ; FIRST = .FALSE. IF ( DEBUG_EXACT .OR. DEBUG_INCLUDE ) PRINT *, & 'NOTE: USED TRANSIENT_HALF_WALL in element ', THIS_EL IF ( N_LP_FLO > 1 ) THEN ! use non-default properties K_e = GET_REAL_LP (1) ! thermal conductivity Rho_e = GET_REAL_LP (2) ! rho, c_p ELSE IF ( EXACT_REALS > 1 ) THEN K_e = GET_REAL_EXACT (1) ! thermal conductivity Rho_e = GET_REAL_EXACT (2) ! rho, c_p END IF PRINT *,'Size MAX_XYZ ', SIZE (MAX_XYZ), SIZE (MIN_XYZ) A = K_e / Rho_e L = 2 * MAX_XYZ (1) !b 2 * MAX_XYZ (1) ! if NOT half-symmetry model END IF PI_X_BY_L = PI * XYZ (1) / L PI_SQ_T = A * PI_SQ * TIME / L **2 TOTAL = 0.d0 DO N = 0, N_SUM M = N + N + 1 TERM_X = SIN ( M * PI_X_BY_L ) / M TERM_T = EXP ( -M * M * PI_SQ_T ) INC = TERM_X * TERM_T TOTAL = TOTAL + INC IF ( ABS( INC ) <= TOLERANCE ) EXIT ! this loop END DO IF ( DEBUG_EXACT ) PRINT *,'TRANSIENT_HALF_WALL ', M, INC IF ( INITIAL_VALUE /= 0.d0 ) TOTAL = TOTAL * INITIAL_VALUE EXACT_SOL (1) = TOTAL * 4 / PI ! end 122_01.my_exact_inc IF ( TOUCH ) PRINT *, XYZ (1) ! TOUCH is never true END SUBROUTINE TRANSIENT_HALF_WALL ! --------------------- EXACT_CASE = 35 ------------------------