FUNCTION VOID THETA_meth(t0, t1, toldold, Taoldold, Taold, Ta, Tboldold, Tbold, Tb, n, T, TP, Qa, Qb, A, r, c, d_1, d0, d1, TOL, THETA, dtstart, Jac_approx) /* Uses a Fortran subroutine for time integration (Theta-method, theta = 0 => Backward Euler (default), theta = 0.5 => Midpoint method) with predictor-corrector procedure and stepsize control. [Qa[t1], Qb[t1]] = F(Ta[t0:t1], Tb[t0:t1]) where Ta and Tb are supplied at beginning and end of timestep (t0,t1) and also at beginning of old timestep and interpolated linearly or quadraticly in between depending on current order k of IDA method. Revisions: 020305 LE Remove previous change, change in int_pol to avoid division by zero and min dt set to 1.0d-10. 020301 LE Skip computations when t1 - t0 <= 1.0d-20. */ LANGUAGE F77 INPUT INT n; FLOAT t0, t1, toldold, Taoldold, Taold, Ta, Tboldold, Tbold, Tb, T, TP, A, r, c, d_1, d0, d1, TOL, THETA, dtstart, Jac_approx; OUTPUT FLOAT T, TP, Qa, Qb, d_1, d0, d1, dtstart; CODE * EXTERNAL THETA_meth(t0, t1, toldold, * Taoldold, Taold, Ta, Tboldold, Tbold, Tb, * nCells, T, work, Qanew, Qbnew, * A, R, C, * work(nCells+1),work(2*nCells+1),work(3*nCells+1), * TOL, THETA, dtstart, Jac_approx) SUBROUTINE THETA_meth(t0, t1, tm1, & Tam1, Ta0, Ta1, Tbm1, Tb0, Tb1, & n, T, TP, Qa, Qb, & A, r, c, & d_1, d0, d1, & TOL, THETA, dtstart, Jac_approx) INTEGER n REAL(KIND=8) t0, t1, tm1, & Tam1, Ta0, Ta1, Tbm1, Tb0, Tb1, & T(n), TP(n), Qa, Qb, & A, r(n+1), c(n), & d_1(n),d0(n),d1(n), & TOL, THETA, dtstart, Jac_approx cle> LOGICAL factor, thru LOGICAL factor, thru, dt_change INTEGER i, di, istep, & k_max, k_integ, step_type, residual REAL(KIND=8) resT(n),dTemp(n), & dTmax, dt, h, k0, k1, time, dtnew, & Ta, Tb, & Te_1, Te0, Te1, dt_1_th * Get information about step type and integration order call step_info(k_max, k_integ, step_type, residual) * If residual = 0 and Jac_approx =/ 1 return if (residual .eq. 0 .and. nint(Jac_approx) .ne. 1) then return end if k_integ = min(k_integ, 2) c write(*,*) 'k_max, k_integ, step_type, residual' c write(*,*) k_max, k_integ, step_type, residual * Initial value computation if (step_type .le. 0) then c dtstart = dtnew dtstart = -1.0 * Heatfluxes Qa = A * (Ta1 - T(1))/ r(1) Qb = A * (T(n) - Tb1)/ r(n+1) RETURN end if * Get time coordinates etc h = t1 - t0 time = t0 thru = .FALSE. di = int (n/10) * If small time step assume no change c if (h .le. 1.0d-2) then c Qa = A * (Ta1 - T(1))/ r(1) c Qb = A * (T(n) - Tb1)/ r(n+1) c RETURN c end if c write (*,*) 'tm1,t0,t1',tm1,t0,t1,h * Initialize c if (step_type .eq. 1) then c dt = h/2 c else dt = h c end if factor = .true. istep = 0 dtnew = dtstart c WRITE(*,98) 'THETA_meth tm1,t0,t1,dtnew,h'tm1,,t0 ,t1,dtnew,h c WRITE(*,*) 'THETA_meth tm1,t0,t1,dtnew,h'tm1,,t0 ,t1,dtnew,h c WRITE(*,*) ' TP(1),Ta0,Tb0',TP(1),Ta0,Tb0 * Repeat from her for each timestep cle> dt_change = .true. 50 CONTINUE istep = istep + 1 * Select timestep cle> if (abs(dt-dtnew) .lt. 0.001*max(dt,dtnew)) then if (.not. dt_change) then factor = .false. c write(*,*) 'factor = ', factor else factor = .true. c write(*,*) 'factor = ', factor end if IF (time + dt .GT. t1) THEN dt = t1 - time thru = .TRUE. factor = .true. ELSEIF (time + dt .EQ. t1) THEN thru = .TRUE. ENDIF cle> c dt = max(1.0d-20,dt) dt = max(1.0d-10,dt) dt_1_th = (1-THETA)*dt * Interpolate boundary temperatures call int_pol(tm1, t0, t1, Tam1, Ta0, Ta1, & k_integ, time+dt_1_th, Ta) call int_pol(tm1, t0, t1, Tbm1, Tb0, Tb1, & k_integ, time+dt_1_th, Tb) * Matrix if (factor) then do i = 1, n d_1(i) = - (1-THETA) /r(i) d0(i) = c(i)/dt + (1-THETA)*(1/r(i)+1/r(i+1)) d1(i) = - (1-THETA) /r(i+1) end do end if * Residual IF (n .EQ. 1) THEN Te_1= Ta Te0 = T(1) + dt_1_th*TP(1) Te1 = Tb resT(1) = c(1)*TP(1) - & (Te_1/r(1) - Te0*(1/r(1)+1/r(2)) + Te1/r(2)) ELSE do i = 1, n IF (i .EQ. 1) THEN Te_1= Ta Te0 = T(1) + dt_1_th*TP(1) Te1 = T(2) + dt_1_th*TP(2) ELSE IF (i .EQ. n) THEN Te_1= T(n-1) + dt_1_th*TP(n-1) Te0 = T(n) + dt_1_th*TP(n) Te1 = Tb ELSE Te_1= T(i-1) + dt_1_th*TP(i-1) Te0 = T(i) + dt_1_th*TP(i) Te1 = T(i+1) + dt_1_th*TP(i+1) ENDIF resT(i) = c(i)*TP(i) - & (Te_1/r(i) - Te0*(1/r(i)+1/r(i+1)) + Te1/r(i+1)) end do ENDIF * Prediction do i = 1,n T(i) = T(i) + dt * TP(i) end do * Solution of tridiagonal system c write(*,*) '0 D0(1):',d0(1),c(i),dt call trisol(d_1,d0,d1, resT, dTemp, n, factor) c write(*,*) '1 D0(1):',d0(1) * Update variables and their time derivatives dTmax = 1e-9 do i = 1,n T(i) = T(i) - dTemp(i) IF (dt .GT. 0) THEN TP(i) = TP(i) - dTemp(i) / dt END IF dTmax = MAX(dTmax,ABS(dTemp(i))) end do * Estimate new timestep dtnew = SQRT(TOL / dTmax) * dt c WRITE(*,*) 'BDF1wall t, dt, dtnew ',t, dt, dtnew c WRITE(7,*) '2 t, dt, dtnew, dTmax ',t, dt, dtnew, dTmax c write (7,99) (dT(i),i=1,n,di) c 99 format (1x,10f6.2) * Update current time time = time + dt c WRITE(*,98) ' t, T(n), dt, dtnew ',time, T(n), dt,dtnew c 98 format (a,10g9.2) c if (istep .eq. 1) then c write(*,*) 'Factor ' , factor, Ta, Tb c write(*,*) 'ResT ' c WRITE(*,990) (resT(i), i = 1,n) c write(*,*) 'dTemp ' c WRITE(*,990) (dTemp(i), i = 1,n) c write(*,*) 'T ' c WRITE(*,990) (T(i), i = 1,n) c end if c WRITE(*,990) time, dt, T(1), T(10), T(20) c 990 format (5f13.5) cle> IF (.NOT. thru) GOTO 50 IF (.NOT. thru) then if (abs(dt-dtnew) .ge. 0.001*max(dt,dtnew)) then dt_change = .true. end if dt = dtnew GOTO 50 end if * Last used time step dtstart = dt * Heatfluxes Qa = A * (Ta1 - T(1))/ r(1) Qb = A * (T(n) - Tb1)/ r(n+1) c WRITE(*,*) ' THETA_meth klar Qa, Qb ', Qa, Qb c write(*,*) ' istep, dt ', istep, dt RETURN END *---------------------------------------------------------- subroutine trisol(a,b,c, h, x, n, factor) c Numerical recipes page 40 INTEGER n,i REAL(KIND=8) a(n),b(n),c(n),h(n),x(n) LOGICAL factor if (n .eq. 1) then x(1) = h(1) / b(1) else if (factor) then c(1) = c(1) / b(1) do i = 2,n-1 b(i) = b(i) - a(i) * c(i-1) c(i) = c(i) / b(i) end do b(n) = b(n) - a(n) * c(n-1) endif x(1) = h(1) / b(1) do i = 2,n x(i) = (h(i) - a(i) * x(i-1)) / b(i) end do do i = n-1,1,-1 x(i) = x(i) - c(i) * x(i+1) end do end if return end subroutine *---------------------------------------------------------- subroutine int_pol( t0, t1, t2, v0, v1, v2, & k_pol, time, value) * Polynomial interpolation of order k_pol (= 1, 2) INTEGER k_pol REAL(KIND=8) t0, t1, t2, v0, v1, v2, & time, value, diff1_21, diff1_32, diff2_31 if ( t2 - t1 .lt. 1.0d-10) then if (abs(time - t1) .lt. abs(time - t2)) then value = v1 else value = v2 end if else if (k_pol .eq. 1) then diff1_21 = (v2 - v1) / (t2 - t1) value = v1 + (time - t1) * diff1_21 else if (k_pol .eq. 2) then if (t1 - t0 .lt. 1.0d-10) then diff1_21 = (v2 - v1) / (t2 - t1) value = v1 + (time - t1) * diff1_21 else diff1_21 = (v1 - v0) / (t1 - t0) diff1_32 = (v2 - v1) / (t2 - t1) diff2_31 = (diff1_32 - diff1_21) / (t2 - t0) value = v1 + (time - t1) * (diff1_21 + & (time - t0) * diff2_31) end if end if c write(*,*) 'k_pol ', k_pol return end subroutine *---------------------------------------------------------- END_CODE EXTENSIONS JACOBIAN THETA_methJ; JACOBIAN_OVER Qa, Qb; JACOBIAN_UNDER Ta, Tb; END_EXTENSIONS END_FUNCTION FUNCTION VOID THETA_methJ(t0, t1, toldold, Taoldold, Taold, Ta, Tboldold, Tbold, Tb, n, T, TP, A, r, c, d_1, d0, d1, TOL, THETA, dtstart, Jac_approx, J) LANGUAGE F77 INPUT INT n; FLOAT t0, t1, toldold, Taoldold, Taold, Ta, Tboldold, Tbold, Tb, T, TP, A, r, c, d_1, d0, d1, TOL, THETA, dtstart, Jac_approx; OUTPUT FLOAT J; CODE SUBROUTINE THETA_methJ(t0, t1, tm1, & Tam1, Ta0, Ta1, Tbm1, Tb0, Tb1, & n, T, TP, & A, r, c, & d_1, d0, d1, & TOL, THETA, dtstart, Jac_approx, & J) * Compute Jacobian: * * Jac_approx = 3 : Zero Jacobian * Jac_approx = 2 : Approximate diagonal Jacobian * Jac_approx = 1 : Jacobian evaluated by integration using * THETA routine * Jac_approx = 0 : Jacobian evaluated with one time step (default) INTEGER n REAL(KIND=8) t0, t1, tm1, & Tam1, Ta0, Ta1, Tbm1, Tb0, Tb1, & T(n), TP(n), & A, r(n+1), c(n), & d_1(n),d0(n),d1(n), & TOL, dtstart, THETA, Jac_approx, & J(2,2) LOGICAL factor INTEGER k_max, k_integ, step_type, residual, & i REAL(KIND=8) resT(n),dTemp(n), & Qa, Qb, & Tam1work, Ta0work, Ta1work, & Tbm1work, Tb0work, Tb1work, & dt, onemth * Zero Jacobian if (nint(Jac_approx) .eq. 3) then J(1,1) = 0 J(2,1) = 0 J(1,2) = 0 J(2,2) = 0 * Approximate diagonal Jacobian elseif (nint(Jac_approx) .eq. 2) then J(1,1) = A/r(1) J(2,1) = 0 J(1,2) = 0 J(2,2) =-A/r(n+1) * Jacobian evaluated by integration using THETA routine elseif (nint(Jac_approx) .eq. 1) then do i = 1, n resT(i) = 0 dTemp(i) = 0 c T(i) = 0 c TP(i) = 0 end do * Perturb Ta and integrate Tam1work= 0 Tbm1work= 0 Ta0work = 0 Tb0work = 0 Ta1work = 1 Tb1work = 0 call THETA_meth(t0, t1, tm1, & Tam1work, Ta0work, Ta1work, & Tbm1work, Tb0work, Tb1work, & n, resT, dTemp, Qa, Qb, c & n, T, TP, Qa, Qb, & A, r, c, & d_1, d0, d1, & TOL, THETA, dtstart, Jac_approx) J(1,1) = Qa J(2,1) = Qb * Perturb Tb and integrate do i = 1, n resT(i) = 0 dTemp(i) = 0 c T(i) = 0 c TP(i) = 0 end do Ta1work = 0 Tb1work = 1 call THETA_meth(t0, t1, tm1, & Tam1work, Ta0work, Ta1work, & Tbm1work, Tb0work, Tb1work, & n, resT, dTemp, Qa, Qb, c & n, T, TP, Qa, Qb, & A, r, c, & d_1, d0, d1, & TOL, THETA, dtstart, Jac_approx) J(1,2) = Qa J(2,2) = Qb * Jacobian evaluated with one time step else * Get information about step type and integration order call step_info(k_max, k_integ, step_type, residual) * Initial value computation if (step_type .le. 0) then J(1,1) = A/r(1) J(2,1) = 0 J(1,2) = 0 J(2,2) =-A/r(n+1) RETURN end if onemth = 1-THETA * Matrix dt = t1 - t0 do i = 1, n d_1(i) = - onemth /r(i) d0(i) = c(i)/dt + onemth*(1/r(i)+1/r(i+1)) d1(i) = - onemth /r(i+1) resT(i) = 0 end do * Perturb Ta and solve resT(1) = -onemth/r(1) factor = .true. call trisol(d_1,d0,d1, resT, dTemp, n, factor) J(1,1) = A/r(1)*(1+dTemp(1)) J(2,1) = -A/r(n+1)*dTemp(n) * Perturb Tb and solve resT(1) = 0 resT(n) = -onemth/r(n+1) factor = .false. call trisol(d_1,d0,d1, resT, dTemp, n, factor) J(1,2) = A/r(1)*dTemp(1) J(2,2) = -A/r(n+1)*(1+dTemp(n)) end if RETURN END END_CODE END_FUNCTION