C*********************************************************************** C Module: profil.f C C Copyright (C) 2000 Mark Drela C C This program is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the License, or C (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public License C along with this program; if not, write to the Free Software C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. C*********************************************************************** SUBROUTINE PRWALL(DSTAR,THETA,UO,RT,MS,CT, BB, & DO, DO_DS, DO_TH, DO_UO, DO_RT, DO_MS, & UI, UI_DS, UI_TH, UI_UO, UI_RT, UI_MS, & HS, HS_DS, HS_TH, HS_UO, HS_RT, HS_MS, & CF, CF_DS, CF_TH, CF_UO, CF_RT, CF_MS, & CD, CD_DS, CD_TH, CD_UO, CD_RT, CD_MS, CD_CT ) IMPLICIT REAL (A-H,M,O-Z) C================================================================ C Returns wall slip velocity and thickness of wall BL profile C C Input: C DSTAR kinematic displacement thickness C THETA kinematic momentum thickness C RT momentum thickness based on ue and THETA C MS Mach^2 based on ue C C UO uo/ue outer velocity; assumed = 1 in this version C C Output: C BB outer profile exponent C DO thickness of profile deck C UI inner "slip" velocity C CF wall skin friction C================================================================ C PARAMETER (N=65) DIMENSION ETA(N), UIP(N), UIP_DP(N), G(N), G_BB(N) C C---- pi/2 , 2/pi DATA HPI, TOPI / 1.570796327 , 0.6366197723 / C DATA T, SQT / 0.28 , 0.5291502622 / C C---- TCON = ( atan(T^1/2) / T^1/2 - 1/(1+T)) / (2T) + 0.5/(1.0 + T) C - atan(1/T^1/2) / 2T^1/2 DATA TCON / -0.3864027035 / C C---- slip velocity coefficient DATA AK / 0.09 / C C---- log-law constants DATA VKAP, VB / 0.40 , 5.0 / C HK = DSTAR/THETA C UO = 1.0 BB = 1.0 C C---- initialize variables CALL CFT(HK,RT,MS,CF,CF_HK,CF_RT,CF_MS) SGN = SIGN( 1.0 , CF ) UT = SGN * SQRT(0.5*ABS(CF)) C UI = MIN( UT/AK * HPI , 0.90 ) DO = HK*THETA / (1.0 - 0.5*(UO+UI)) C EBK = EXP(-VB*VKAP) C DO 1000 ITER=1, 12 C SGN = SIGN( 1.0 , UT ) C C------ set d+ = DP(UT DO ; RT TH) DP = SGN * UT*RT*DO/THETA DP_DO = SGN * UT*RT /THETA DP_UT = SGN * RT*DO/THETA C DP_TH = -DP/THETA DP_RT = SGN * UT *DO/THETA DP_MS = 0.0 C C------ estimate inner profile edge velocity Ui+ using log-law UPE = LOG(DP)/VKAP + VB C C------ converge exact Ui+ using Spalding formula DO 10 ITUP=1, 5 UK = UPE*VKAP ARG = UK - VB*VKAP EXU = EXP(ARG) REZ = UPE + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) & - DP DP_U = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP C IF(ABS(REZ/DP) .LT. 1.0E-5) GO TO 11 C DUPE = -REZ/DP_U UPE = UPE + DUPE 10 CONTINUE WRITE(*,*) 'PRWALL: Ue+ convergence failed, Res =', REZ/DP 11 CONTINUE C UPE_DP = 1.0/DP_U C C 2 2 3 3 C------ set d y+/du+ and d y+/du+ at BL edge DP_UU = (EXU - EBK*(1.0 + UK))*VKAP**2 DP_UUU = (EXU - EBK )*VKAP**3 C C------ set du+/dy+ at BL edge UPD = 1.0/DP_U UPD_DP = (-1.0/DP_U**3) * DP_UU C C 2 2 C------ set d u+/dy+ at BL edge CCC UPD_DP = (-1.0/DP_U**3) * DP_UU UPDD = UPD_DP UPDD_DP = (-1.0/DP_U**4) * DP_UUU & + ( 3.0/DP_U**5) * DP_UU**2 C C------ set coefficients for Spalding profile correction polynomial DC2 = 0.5*DP*DP*UPDD - DP*UPD DC2_DP = DP*UPDD - UPD & + 0.5*DP*DP*UPDD_DP - DP*UPD_DP C DC3 = -( DP*DP*UPDD - DP*UPD ) / 3.0 DC3_DP = -(2.0 *DP*UPDD - UPD ) / 3.0 & -( DP*DP*UPDD_DP - DP*UPD_DP) / 3.0 C C------ set outer profile amplitude DUO DUO = UO - UT*(UPE + DC2 + DC3 ) DUO_DP = - UT*(UPE_DP + DC2_DP + DC3_DP) C DUO_UT = - (UPE + DC2 + DC3 ) & + DUO_DP*DP_UT DUO_DO = DUO_DP*DP_DO C DUO_TH = DUO_DP*DP_TH DUO_RT = DUO_DP*DP_RT DUO_MS = DUO_DP*DP_MS c c write(*,*) 'dUo', duo, duo_dp, duo_ut, duo_do c read(*,*) ddo, dut c if(ddo.ne.0.0 .or. dut.ne.0.0) then c do = do+ddo c ut = ut+dut c write(*,*) 'new', duo + duo_do*ddo + duo_ut*dut c go to 666 c endif C C------ set wake profile coefficients BB1 = 3.0*(BB +2.0)*(BB+3.0)/(BB+7.0) BB1_BB = 3.0*(BB*2.0+5.0 )/(BB+7.0) - BB1/(BB+7.0) BB2 = -5.0*(BB +1.0)*(BB+3.0)/(BB+7.0) BB2_BB = -5.0*(BB*2.0+4.0 )/(BB+7.0) - BB2/(BB+7.0) BB3 = 2.0*(BB +1.0)*(BB+2.0)/(BB+7.0) BB3_BB = 2.0*(BB*2.0+3.0 )/(BB+7.0) - BB3/(BB+7.0) C C------ fill eta coordinate and inner profile arrays CCC EXUPE = EXP(UPE*VKAP - VB*VKAP) EXUPE = EXU C DEXU = (EXUPE - EBK)/FLOAT(N-1) C I = 1 UIP(I) = 0.0 UIP_DP(I) = 0.0 G(I) = 0.0 G_BB(I) = 0.0 C DO 20 I=2, N ccc EXU = EBK + DEXU*FLOAT(I-1) EXU = EBK + (DEXU - 0.75*DEXU*FLOAT(N-I)/FLOAT(N-1)) & *FLOAT(I-1) C CCC UK = UP*VKAP UK = LOG(EXU) + VB*VKAP C UP = UK/VKAP C C-------- set "inverse" Spalding profile y+(u+) and derivatives YP = UP + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) YP_U = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP YP_UU = (EXU - EBK*(1.0 + UK ))*VKAP**2 C ET = YP/DP C C-------- set final inner profile (fudged Spalding) UIP(I) = UP + DC2 *ET**2 + DC3 *ET**3 UIP_DP(I) = DC2_DP*ET**2 + DC3_DP*ET**3 C ccc UIPD(I) = 1.0/YP_U + 2.0*DC2 *ET + 3.0*DC3 *ET**2 ccc UIPD_DP(I) = (-1.0/YP_U**3)*YPUU ccc & + 2.0*DC2_DP*ET + 3.0*DC3_DP*ET**2 C C-------- set outer profile ETB = ET**BB ALE = LOG(ET) C ccc G(I) = 2.0*ETB - ETB**2 ccc G_BB(I) = (2.0*ETB - 2.0*ETB**2)*ALE C G(I) = (BB1 *ET + BB2 *ET**2 + BB3 *ET**3)*ETB G_BB(I) = (BB1_BB*ET + BB2_BB*ET**2 + BB3_BB*ET**3)*ETB & + G(I)*ALE C ETA(I) = ET C 20 CONTINUE C C DSN = 0.0 DSN_DO = 0.0 DSN_UT = 0.0 DSN_BB = 0.0 C DSN_TH = 0.0 DSN_RT = 0.0 DSN_MS = 0.0 C C THN = 0.0 THN_DO = 0.0 THN_UT = 0.0 THN_BB = 0.0 C THN_TH = 0.0 THN_RT = 0.0 THN_MS = 0.0 C c TSN = 0.0 c TSN_DO = 0.0 c TSN_UT = 0.0 c TSN_BB = 0.0 c c TSN_TH = 0.0 c TSN_RT = 0.0 c TSN_MS = 0.0 C C------ perform integration DO 100 I=1, N-1 DETA = ETA(I+1) - ETA(I) GA = 0.5*(G(I+1) + G(I) ) GA_BB = 0.5*(G_BB(I+1) + G_BB(I)) C UIPA = 0.5*(UIP(I+1) + UIP(I) ) UIPA_DP = 0.5*(UIP_DP(I+1) + UIP_DP(I)) C U = UT*UIPA + DUO *GA U_DP = UT*UIPA_DP C U_DO = DUO_DO*GA + U_DP*DP_DO U_UT = UIPA + DUO_UT*GA + U_DP*DP_UT U_BB = DUO *GA_BB C U_TH = DUO_TH*GA + U_DP*DP_TH U_RT = DUO_RT*GA + U_DP*DP_RT U_MS = DUO_MS*GA + U_DP*DP_MS C C DSN = DSN + (1.0 - U )*DETA DSN_DO = DSN_DO - U_DO *DETA DSN_UT = DSN_UT - U_UT *DETA DSN_BB = DSN_BB - U_BB *DETA C DSN_TH = DSN_TH - U_TH *DETA DSN_RT = DSN_RT - U_RT *DETA DSN_MS = DSN_MS - U_MS *DETA C C THN = THN + (U - U*U) *DETA THN_DO = THN_DO + (1.0 - 2.0*U)*U_DO*DETA THN_UT = THN_UT + (1.0 - 2.0*U)*U_UT*DETA THN_BB = THN_BB + (1.0 - 2.0*U)*U_BB*DETA C THN_TH = THN_TH + (1.0 - 2.0*U)*U_TH*DETA THN_RT = THN_RT + (1.0 - 2.0*U)*U_RT*DETA THN_MS = THN_MS + (1.0 - 2.0*U)*U_MS*DETA C c TSN = TSN + (U - U*U*U) *DETA c TSN_DO = TSN_DO + (1.0 - 3.0*U*U)*U_DO*DETA c TSN_UT = TSN_UT + (1.0 - 3.0*U*U)*U_UT*DETA c TSN_BB = TSN_BB + (1.0 - 3.0*U*U)*U_BB*DETA C c TSN_TH = TSN_TH + (1.0 - 3.0*U*U)*U_TH*DETA c TSN_RT = TSN_RT + (1.0 - 3.0*U*U)*U_RT*DETA c TSN_MS = TSN_MS + (1.0 - 3.0*U*U)*U_MS*DETA C 100 CONTINUE C C------ set up 2x2 system for DO UT REZ1 = DO*DSN - THETA*HK A11 = DO*DSN_DO + DSN A12 = DO*DSN_UT cc A12 = DO*DSN_BB C REZ2 = DO*THN - THETA A21 = DO*THN_DO + THN A22 = DO*THN_UT cc A22 = DO*THN_BB C cc IF(ABS(REZ1/THETA) .LT. 2.0E-5 .AND. cc & ABS(REZ2/THETA) .LT. 2.0E-5 ) GO TO 1010 IF(ABS(REZ1/THETA) .LT. 1.0E-3 .AND. & ABS(REZ2/THETA) .LT. 1.0E-3 ) GO TO 1010 C DET = A11*A22 - A12*A21 B11 = A22/DET B12 = -A12/DET B21 = -A21/DET B22 = A11/DET C DDO = -(B11*REZ1 + B12*REZ2) DUT = -(B21*REZ1 + B22*REZ2) cc DBB = -(B21*REZ1 + B22*REZ2) C DMAX = MAX( ABS(DDO/DO) , ABS(DUT/0.05) ) cc DMAX = MAX( ABS(DDO/DO) , ABS(DBB/BB ) ) RLX = 1.0 IF(DMAX.GT.0.5) RLX = 0.5/DMAX C DO = DO + RLX*DDO UT = UT + RLX*DUT cc BB = BB + RLX*DBB c cc write(*,*) iter, do, ut, rez1, rez2 cc write(*,*) iter, do, bb, rez1, rez2 C 1000 CONTINUE C WRITE(*,*) 'PRWALL: Convergence failed. Res =', REZ1, REZ2 C 1010 CONTINUE C C C CCC REZ1 = DO*DSN - THETA*HK Z1_HK = - THETA Z1_TH = DO*DSN_TH - HK Z1_RT = DO*DSN_RT Z1_MS = DO*DSN_MS C CCC REZ2 = DO*THN - THETA Z2_HK = 0.0 Z2_TH = DO*THN_TH - 1.0 Z2_RT = DO*THN_RT Z2_MS = DO*THN_MS C DO_HK = -(B11*Z1_HK + B12*Z2_HK) DO_TH = -(B11*Z1_TH + B12*Z2_TH) DO_RT = -(B11*Z1_RT + B12*Z2_RT) DO_MS = -(B11*Z1_MS + B12*Z2_MS) C UT_HK = 0.0 UT_TH = 0.0 UT_RT = 0.0 UT_MS = 0.0 C BB_HK = 0.0 BB_TH = 0.0 BB_RT = 0.0 BB_MS = 0.0 C UT_HK = -(B21*Z1_HK + B22*Z2_HK) UT_TH = -(B21*Z1_TH + B22*Z2_TH) UT_RT = -(B21*Z1_RT + B22*Z2_RT) UT_MS = -(B21*Z1_MS + B22*Z2_MS) C cc BB_HK = -(B21*Z1_HK + B22*Z2_HK) cc BB_TH = -(B21*Z1_TH + B22*Z2_TH) cc BB_RT = -(B21*Z1_RT + B22*Z2_RT) cc BB_MS = -(B21*Z1_MS + B22*Z2_MS) C C C---- set and linearize Cf CF = SGN*2.0*UT**2 CF_UT = SGN*4.0*UT CF_DO = 0.0 C CF_HK = CF_UT*UT_HK + CF_DO*DO_HK CF_TH = CF_UT*UT_TH + CF_DO*DO_TH CF_RT = CF_UT*UT_RT + CF_DO*DO_RT CF_MS = CF_UT*UT_MS + CF_DO*DO_MS C C C---- set and linearize "slip" velocity UI = UI( DUO(DO UT TH RT MS) ) UI = UO - DUO UI_UT = - DUO_UT UI_DO = - DUO_DO C UI_HK = UI_UT*UT_HK + UI_DO*DO_HK UI_TH = UI_UT*UT_TH + UI_DO*DO_TH - DUO_TH UI_RT = UI_UT*UT_RT + UI_DO*DO_RT - DUO_RT UI_MS = UI_UT*UT_MS + UI_DO*DO_MS - DUO_MS C RETURN END ! PRWALL SUBROUTINE UWALL(TH,UO,DO,UI,RT,CF,BB, Y,U,N) C------------------------------------------ C Returns wall BL profile U(Y). C C Input: C TH kinematic momentum thickness C UO uo/ue outer velocity (= 1 for normal BL) C DO BL thickness C UI inner "slip" velocity C RT momentum thickness based on ue and THETA C CF wall skin friction C BB outer profile exponent C N number of profile array points C C Output: C Y(i) normal coordinate array C U(i) u/ue velocity profile array C------------------------------------------- C IMPLICIT REAL (A-H,M,O-Z) DIMENSION Y(N), U(N) DATA HPI / 1.570796327 / DATA AK / 0.09 / C C---- log-law constants DATA VKAP, VB / 0.40 , 5.0 / C EBK = EXP(-VB*VKAP) C SGN = SIGN( 1.0 , CF ) UT = SGN * SQRT(0.5*ABS(CF)) C C C---- set d+ = DP(UT DO ; RT TH) DP = SGN * UT*RT*DO/TH C C---- estimate inner profile edge velocity Ui+ using log-law UPE = LOG(DP)/VKAP + VB C C---- converge exact Ui+ using Spalding formula DO 10 ITUP=1, 5 UK = UPE*VKAP ARG = UK - VB*VKAP EXU = EXP(ARG) REZ = UPE + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) & - DP DP_U = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP C IF(ABS(REZ/DP) .LT. 1.0E-5) GO TO 11 C DUPE = -REZ/DP_U UPE = UPE + DUPE 10 CONTINUE WRITE(*,*) 'UWALL: Ue+ convergence failed, Res =', REZ/DP 11 CONTINUE C C 2 2 3 3 C---- set d y+/du+ and d y+/du+ at BL edge DP_UU = (EXU - EBK*(1.0 + UK))*VKAP**2 DP_UUU = (EXU - EBK )*VKAP**3 C C 2 2 C---- set du+/dy+ and d u+/dy+ at BL edge UPD = 1.0/DP_U UPDD = (-1.0/DP_U**3) * DP_UU C C---- set coefficients for Spalding profile correction polynomial DC2 = 0.5*DP*DP*UPDD - DP*UPD DC3 = -( DP*DP*UPDD - DP*UPD ) / 3.0 C C---- set outer profile amplitude DUO DUO = UO - UT*(UPE + DC2 + DC3) C C BB1 = 3.0*(BB +2.0)*(BB+3.0)/(BB+7.0) BB2 = -5.0*(BB +1.0)*(BB+3.0)/(BB+7.0) BB3 = 2.0*(BB +1.0)*(BB+2.0)/(BB+7.0) C c NE = (N*9)/10 NE = N C C---- fill Y coordinate and U profile arrays CCC EXUPE = EXP(UPE*VKAP - VB*VKAP) EXUPE = EXU C DEXU = (EXUPE - EBK)/FLOAT(NE-1) C I = 1 Y(I) = 0.0 U(I) = 0.0 DO 20 I=2, NE ccc EXU = EBK + DEXU*FLOAT(I-1) EXU = EBK + (DEXU - 0.75*DEXU*FLOAT(NE-I)/FLOAT(NE-1)) & *FLOAT(I-1) C CCC UK = UP*VKAP UK = LOG(EXU) + VB*VKAP C UP = UK/VKAP C C------ set "inverse" Spalding profile y+(u+) YP = UP + EXU - EBK*(1.0 + UK + UK*UK/2.0 + UK*UK*UK/6.0) YP_UP = 1.0 + (EXU - EBK*(1.0 + UK + UK*UK/2.0))*VKAP C ET = YP/DP C C------ set final inner profile (fudged Spalding) UIP = UP + DC2*ET**2 + DC3*ET**3 C C------ set outer profile SQE = SQRT(ET) ETB = ET**BB C ccc G = 2.0*ETB - ETB**2 C G = (BB1 *ET + BB2 *ET**2 + BB3 *ET**3)*ETB C Y(I) = ET*DO U(I) = UT*UIP + DUO*G c 20 CONTINUE C c DETA = 0.1 / FLOAT(N - NE - 1) c DO 300 I=NE+1, N c ETA = 1.0 + DETA*FLOAT(I-NE) c Y(I) = DO*ETA c U(I) = 1.0 c 300 CONTINUE C RETURN END ! UWALL SUBROUTINE FS(INORM,ISPEC,BSPEC,HSPEC,N,ETAE,GEO,ETA,F,U,S,DELTA) IMPLICIT REAL (A-H,M,O-Z) DIMENSION ETA(N), F(N), U(N), S(N) C----------------------------------------------------- C Routine for solving the Falkner-Skan equation. C C Input: C ------ C INORM 1: eta = y / sqrt(vx/Ue) "standard" Falkner-Skan coordinate C 2: eta = y / sqrt(2vx/(m+1)Ue) Hartree's coordinate C 3: eta = y / Theta momentum thickness normalized coordinate C ISPEC 1: BU = x/Ue dUe/dx ( = "m") specified C 2: H12 = Dstar/Theta specified C BSPEC specified pressure gradient parameter (if ISPEC = 1) C HSPEC specified shape parameter of U profile (if ISPEC = 2) C N total number of points in profiles C ETAE edge value of normal coordinate C GEO exponential stretching factor for ETA: C C Output: C ------- C BSPEC calculated pressure gradient parameter (if ISPEC = 2) C HSPEC calculated shape parameter of U profile (if ISPEC = 1) C ETA normal BL coordinate C F,U,S Falkner Skan profiles C DELTA normal coordinate scale y = eta * Delta C----------------------------------------------------- C PARAMETER (NMAX=257,NRMAX=3) DIMENSION A(3,3,NMAX),B(3,3,NMAX),C(3,3,NMAX), & R(3,NRMAX,NMAX) C C---- set number of righthand sides. DATA NRHS / 3 / C ITMAX = 20 C IF(N.GT.NMAX) STOP 'FS: Array overflow.' C PI = 4.0*ATAN(1.0) C CCC if(u(n) .ne. 0.0) go to 9991 c C---- initialize H or BetaU with empirical curve fits IF(ISPEC.EQ.1) THEN H = 2.6 BU = BSPEC ELSE H = HSPEC IF(H .LE. 2.17) THEN WRITE(*,*) 'FS: Specified H too low' H = 2.17 ENDIF BU = (0.058*(H-4.0)**2/(H-1.0) - 0.068) / (6.54*H - 14.07) * H**2 IF(H .GT. 4.0) BU = MIN( BU , 0.0 ) ENDIF C C---- initialize TN = Delta^2 Ue / vx IF(INORM.EQ.3) THEN TN = (6.54*H - 14.07) / H**2 ELSE TN = 1.0 ENDIF C C---- set eta array DETA = 1.0 ETA(1) = 0.0 DO 5 I=2, N ETA(I) = ETA(I-1) + DETA DETA = GEO*DETA 5 CONTINUE C DO 6 I=1, N ETA(I) = ETA(I) * ETAE/ETA(N) 6 CONTINUE C C C---- initial guess for profiles using a sine loop for U for half near wall IF(H .LE. 3.0) THEN C IF(INORM.EQ.3) THEN ETJOIN = 7.3 ELSE ETJOIN = 5.0 ENDIF C EFAC = 0.5*PI/ETJOIN DO 10 I=1, N U(I) = SIN(EFAC*ETA(I)) F(I) = 1.0/EFAC * (1.0 - COS(EFAC*ETA(I))) S(I) = EFAC*COS(EFAC*ETA(I)) IF(ETA(I) .GT. ETJOIN) GO TO 11 10 CONTINUE 11 CONTINUE IJOIN = I C C----- constant U for outer half DO 12 I=IJOIN+1, N U(I) = 1.0 F(I) = F(IJOIN) + ETA(I) - ETA(IJOIN) S(I) = 0. 12 CONTINUE C ELSE C IF(INORM.EQ.3) THEN ETJOIN1 = 0.0 ETJOIN2 = 8.0 IF(H .GT. 4.0) THEN ETJOIN1 = H - 4.0 ETJOIN2 = ETJOIN1 + 8.0 ENDIF ELSE ETJOIN1 = 0.0 ETJOIN2 = 8.0 ENDIF C DO 13 I=1, N U(I) = 0.0 S(I) = 0.0 F(I) = 0.0 IF(ETA(I) .GE. ETJOIN1) GO TO 14 13 CONTINUE 14 CONTINUE IJOIN = I C EFAC = 0.5*PI/(ETJOIN2-ETJOIN1) DO 15 I=IJOIN+1, N EBAR = ETA(I) - ETJOIN1 U(I) = 0.5 - 0.5*COS(2.0*EFAC*EBAR) F(I) = 0.5*EBAR - 0.25/EFAC * SIN(2.0*EFAC*EBAR) S(I) = EFAC*SIN(2.0*EFAC*EBAR) IF(ETA(I) .GE. ETJOIN2) GO TO 16 15 CONTINUE 16 CONTINUE IJOIN = I C C----- constant U for outer half DO 17 I=IJOIN+1, N U(I) = 1.0 F(I) = F(IJOIN) + ETA(I) - ETA(IJOIN) S(I) = 0. 17 CONTINUE C ENDIF c 9991 continue C C C---- Newton iteration loop DO 100 ITER=1, ITMAX C C------ zero out A,B,C blocks and righthand sides R DO 20 I=1, N DO 201 II=1,3 DO 2001 III=1,3 A(II,III,I) = 0. B(II,III,I) = 0. C(II,III,I) = 0. 2001 CONTINUE R(II,1,I) = 0. R(II,2,I) = 0. R(II,3,I) = 0. 201 CONTINUE 20 CONTINUE C C................................................... C A(1,1,1) = 1.0 A(2,2,1) = 1.0 A(3,2,N) = 1.0 R(1,1,1) = F(1) R(2,1,1) = U(1) R(3,1,N) = U(N) - 1.0 C IF(INORM.EQ.2) THEN BETU = 2.0*BU/(BU+1.0) BETU_BU = (2.0 - BETU/(BU+1.0))/(BU+1.0) BETN = 1.0 BETN_BU = 0.0 ELSE BETU = BU BETU_BU = 1.0 BETN = 0.5*(1.0 + BU) BETN_BU = 0.5 ENDIF C DO 30 I=1,N-1 C DETA = ETA(I+1) - ETA(I) R(1,1,I+1) = F(I+1) - F(I) - 0.5*DETA*(U(I+1)+U(I)) R(2,1,I+1) = U(I+1) - U(I) - 0.5*DETA*(S(I+1)+S(I)) R(3,1,I) = S(I+1) - S(I) & + TN * ( BETN*DETA*0.5*(F(I+1)*S(I+1) + F(I)*S(I)) & + BETU*DETA*(1.0 - 0.5*(U(I+1)**2 + U(I)**2)) ) C A(3,1,I) = TN * BETN*0.5*DETA*S(I) C(3,1,I) = TN * BETN*0.5*DETA*S(I+1) A(3,2,I) = -TN * BETU *DETA*U(I) C(3,2,I) = -TN * BETU *DETA*U(I+1) A(3,3,I) = TN * BETN*0.5*DETA*F(I) - 1.0 C(3,3,I) = TN * BETN*0.5*DETA*F(I+1) + 1.0 C B(1,1,I+1) = -1.0 A(1,1,I+1) = 1.0 B(1,2,I+1) = -0.5*DETA A(1,2,I+1) = -0.5*DETA C B(2,2,I+1) = -1.0 A(2,2,I+1) = 1.0 B(2,3,I+1) = -0.5*DETA A(2,3,I+1) = -0.5*DETA C R(3,2,I) = TN & * ( BETN_BU*DETA*0.5*(F(I+1)*S(I+1) + F(I)*S(I)) & + BETU_BU*DETA*(1.0 - 0.5*(U(I+1)**2 + U(I)**2))) R(3,3,I) = ( BETN*DETA*0.5*(F(I+1)*S(I+1) + F(I)*S(I)) & + BETU*DETA*(1.0 - 0.5*(U(I+1)**2 + U(I)**2)) ) C 30 CONTINUE C........................................................... C C C---- solve Newton system for the three solution vectors CALL B3SOLV(A,B,C,R,N,NRHS,NRMAX) C C C---- calculate and linearize Dstar, Theta, in computational space DSI = 0. DSI1 = 0. DSI2 = 0. DSI3 = 0. C THI = 0. THI1 = 0. THI2 = 0. THI3 = 0. C DO 40 I=1,N-1 US = U(I) + U(I+1) DETA = ETA(I+1) - ETA(I) C DSI = DSI + (1.0 - 0.5*US)*DETA DSI_US = -0.5*DETA C THI = THI + (1.0 - 0.5*US)*0.5*US*DETA THI_US = (0.5 - 0.5*US)*DETA C DSI1 = DSI1 + DSI_US*(R(2,1,I) + R(2,1,I+1)) DSI2 = DSI2 + DSI_US*(R(2,2,I) + R(2,2,I+1)) DSI3 = DSI3 + DSI_US*(R(2,3,I) + R(2,3,I+1)) C THI1 = THI1 + THI_US*(R(2,1,I) + R(2,1,I+1)) THI2 = THI2 + THI_US*(R(2,2,I) + R(2,2,I+1)) THI3 = THI3 + THI_US*(R(2,3,I) + R(2,3,I+1)) 40 CONTINUE C C IF(ISPEC.EQ.1) THEN C C----- set and linearize Bu = Bspec residual R1 = BSPEC - BU Q11 = 1.0 Q12 = 0.0 C ELSE C C----- set and linearize H = Hspec residual R1 = DSI - HSPEC*THI & -DSI1 + HSPEC*THI1 Q11 = -DSI2 + HSPEC*THI2 Q12 = -DSI3 + HSPEC*THI3 C ENDIF C C IF(INORM.EQ.3) THEN C C----- set and linearize normalized Theta = 1 residual R2 = THI - 1.0 & -THI1 Q21 = -THI2 Q22 = -THI3 C ELSE C C----- set eta scaling coefficient to unity R2 = 1.0 - TN Q21 = 0.0 Q22 = 1.0 C ENDIF C C DET = Q11*Q22 - Q12*Q21 DBU = -(R1 *Q22 - Q12*R2 ) / DET DTN = -(Q11*R2 - R1 *Q21) / DET C C C---- calculate changes in F,U,S, and the max and rms change RMAX = 0. RMS = 0. DO 50 I=1,N DF = -R(1,1,I) - DBU*R(1,2,I) - DTN*R(1,3,I) DU = -R(2,1,I) - DBU*R(2,2,I) - DTN*R(2,3,I) DS = -R(3,1,I) - DBU*R(3,2,I) - DTN*R(3,3,I) C RMAX = MAX(RMAX,ABS(DF),ABS(DU),ABS(DS)) RMS = DF**2 + DU**2 + DS**2 + RMS 50 CONTINUE RMS = SQRT(RMS/(3.0*FLOAT(N) + 3.0)) C RMAX = MAX(RMAX,ABS(DBU/0.5),ABS(DTN/TN)) C C---- set underrelaxation factor if necessary by limiting max change to 0.5 RLX = 1.0 IF(RMAX.GT.0.5) RLX = 0.5/RMAX C C---- update F,U,S DO 60 I=1,N DF = -R(1,1,I) - DBU*R(1,2,I) - DTN*R(1,3,I) DU = -R(2,1,I) - DBU*R(2,2,I) - DTN*R(2,3,I) DS = -R(3,1,I) - DBU*R(3,2,I) - DTN*R(3,3,I) C F(I) = F(I) + RLX*DF U(I) = U(I) + RLX*DU S(I) = S(I) + RLX*DS 60 CONTINUE C C---- update BetaU and Theta BU = BU + RLX*DBU TN = TN + RLX*DTN C C---- check for convergence IF(ITER.GT.3 .AND. RMS.LT.1.E-5) GO TO 105 C 100 CONTINUE WRITE(*,*) 'FS: Convergence failed' C 105 CONTINUE C HSPEC = DSI/THI BSPEC = BU C DELTA = SQRT(TN) C RETURN C C The END SUBROUTINE B3SOLV(A,B,C,R,N,NRHS,NRMAX) IMPLICIT REAL (A-H,M,O-Z) DIMENSION A(3,3,N), B(3,3,N), C(3,3,N), R(3,NRMAX,N) C ********************************************************************** C This routine solves a 3x3 block-tridiagonal system with an arbitrary C number of righthand sides by a standard block elimination scheme. C The solutions are returned in the Rj vectors. C C |A C ||d| |R..| C |B A C ||d| |R..| C | B . . ||.| = |R..| C | . . C||.| |R..| C | B A||d| |R..| C Mark Drela 10 March 86 C ********************************************************************** C CCC** Forward sweep: Elimination of lower block diagonal (B's). DO 1 I=1, N C IM = I-1 C C------ don't eliminate B1 block because it doesn't exist IF(I.EQ.1) GO TO 12 C C------ eliminate Bi block, thus modifying Ai and Ci blocks DO 11 K=1, 3 DO 111 L=1, 3 A(K,L,I) = A(K,L,I) & - (B(K,1,I)*C(1,L,IM) + B(K,2,I)*C(2,L,IM) + B(K,3,I)*C(3,L,IM)) 111 CONTINUE DO 112 L=1, NRHS R(K,L,I) = R(K,L,I) & - (B(K,1,I)*R(1,L,IM) + B(K,2,I)*R(2,L,IM) + B(K,3,I)*R(3,L,IM)) 112 CONTINUE 11 CONTINUE C C -1 CCC---- multiply Ci block and righthand side Ri vectors by (Ai) C using Gaussian elimination. C 12 DO 13 KPIV=1, 2 KP1 = KPIV+1 C C-------- find max pivot index KX KX = KPIV DO 131 K=KP1, 3 IF(ABS(A(K,KPIV,I)) .LE. ABS(A(KX,KPIV,I))) THEN GO TO 131 ELSE GO TO 1311 ENDIF 1311 KX = K 131 CONTINUE C IF(A(KX,KPIV,I).EQ.0.0) THEN WRITE(*,*) 'Singular A block, i = ',I STOP ENDIF C PIVOT = 1.0/A(KX,KPIV,I) C C-------- switch pivots A(KX,KPIV,I) = A(KPIV,KPIV,I) C C-------- switch rows & normalize pivot row DO 132 L=KP1, 3 TEMP = A(KX,L,I)*PIVOT A(KX,L,I) = A(KPIV,L,I) A(KPIV,L,I) = TEMP 132 CONTINUE C DO 133 L=1, 3 TEMP = C(KX,L,I)*PIVOT C(KX,L,I) = C(KPIV,L,I) C(KPIV,L,I) = TEMP 133 CONTINUE C DO 134 L=1, NRHS TEMP = R(KX,L,I)*PIVOT R(KX,L,I) = R(KPIV,L,I) R(KPIV,L,I) = TEMP 134 CONTINUE CB C-------- forward eliminate everything DO 135 K=KP1, 3 DO 1351 L=KP1, 3 A(K,L,I) = A(K,L,I) - A(K,KPIV,I)*A(KPIV,L,I) 1351 CONTINUE C(K,1,I) = C(K,1,I) - A(K,KPIV,I)*C(KPIV,1,I) C(K,2,I) = C(K,2,I) - A(K,KPIV,I)*C(KPIV,2,I) C(K,3,I) = C(K,3,I) - A(K,KPIV,I)*C(KPIV,3,I) DO 1352 L=1, NRHS R(K,L,I) = R(K,L,I) - A(K,KPIV,I)*R(KPIV,L,I) 1352 CONTINUE 135 CONTINUE C 13 CONTINUE C C------ solve for last row IF(A(3,3,I).EQ.0.0) THEN WRITE(*,*) 'Singular A block, i = ',I STOP ENDIF PIVOT = 1.0/A(3,3,I) C(3,1,I) = C(3,1,I)*PIVOT C(3,2,I) = C(3,2,I)*PIVOT C(3,3,I) = C(3,3,I)*PIVOT DO 14 L=1, NRHS R(3,L,I) = R(3,L,I)*PIVOT 14 CONTINUE C C------ back substitute everything DO 15 KPIV=2, 1, -1 KP1 = KPIV+1 DO 151 K=KP1, 3 C(KPIV,1,I) = C(KPIV,1,I) - A(KPIV,K,I)*C(K,1,I) C(KPIV,2,I) = C(KPIV,2,I) - A(KPIV,K,I)*C(K,2,I) C(KPIV,3,I) = C(KPIV,3,I) - A(KPIV,K,I)*C(K,3,I) DO 1511 L=1, NRHS R(KPIV,L,I) = R(KPIV,L,I) - A(KPIV,K,I)*R(K,L,I) 1511 CONTINUE 151 CONTINUE 15 CONTINUE 1 CONTINUE C CCC** Backward sweep: Back substitution using upper block diagonal (Ci's). DO 2 I=N-1, 1, -1 IP = I+1 DO 21 L=1, NRHS DO 211 K=1, 3 R(K,L,I) = R(K,L,I) & - (R(1,L,IP)*C(K,1,I) + R(2,L,IP)*C(K,2,I) + R(3,L,IP)*C(K,3,I)) 211 CONTINUE 21 CONTINUE 2 CONTINUE C RETURN END ! B3SOLV