|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 C 00002 C 00003 C ********************************************************************** 00004 C * * 00005 C * ADAMS-BASHFORTH-MOULTON (abmlib.f) * 00006 C * * 00007 C * Library module containing an Adams-Bashforth-Moulton predictor- * 00008 C * corrector integrator. * 00009 C * * 00010 C * Contains: * 00011 C * ABM: The ABM integration routine. * 00012 C * AbsABMError: Computes the absolute error of the current ABM * 00013 C * integration step. * 00014 C * RelABMError: Computes the relative error of the current ABM * 00015 C * integration step. * 00016 C * WriteABMStep: In the current form WriteABMStep is just a dummy * 00017 C * routine doing absolutely nothing; it is included * 00018 C * here for formal reasons, but can (if desired) be * 00019 C * extended easily such that it writes the size and * 00020 C * error of the current integration step to a file. * 00021 C * ABMErrorMsg: Returns for a given error number a corresponding * 00022 C * error message. * 00023 C * * 00024 C ********************************************************************** 00025 00026 00027 C ********************************************************************** 00028 C * * 00029 C * SUBROUTINE ABM * 00030 C * * 00031 C * Integrates a system of complex first order differential equations * 00032 C * employing the Adams-Bashforth-Moulton predictor-corrector method. * 00033 C * The ABM routine runs with variable step sizes and (except for the * 00034 C * beginning) a fixed integration order. The ODE is of the form * 00035 C * dPsi/dt = Func(AbsTime,Psi) =: DtPsi. All computations are * 00036 C * performed with double precision. The routine allows a previous * 00037 C * integration to be resumed (e. g. after the integration had been * 00038 C * stopped to write Psi to a file), as long as the AuxPsi array * 00039 C * hasn't been overwritten and the integration order and error * 00040 C * tolerance hasn't been changed. To do so simply set the restart * 00041 C * flag "RestartABM" to ".false.". * 00042 C * * 00043 C * Input parameters: * 00044 C * Psi: The (complex) initial-value vector. * 00045 C * DtPsi: Time derivative of the initial-value vector (must * 00046 C * be passed to the integrator regardless of the * 00047 C * flag "RestartABM"). * 00048 C * PsiDim Length of Psi and DtPsi vectors. * 00049 C * IntPeriod: Lenght of time interval to be integrated. * 00050 C * AbsTime: Absolute time, i. e. Psi_initial=Psi(AbsTime). * 00051 C * IntOrder: Desired integration order. * 00052 C * InitStep: Size of first integration step. * 00053 C * TolError: Maximum error that is tolerated. * 00054 C * RestartABM: Restart flag; if false, a previous integration is * 00055 C * continued, otherwise a new integration is started * 00056 C * (see comment above for details). * 00057 C * * 00058 C * Output parameters: * 00059 C * Psi: Solution of the ODE at time AbsTime+IntPeriod. * 00060 C * DtPsi: DtPsi is undetermined. * 00061 C * Steps: Counter for the overall number of integration * 00062 C * steps. (Note that Steps isn't set to zero at the * 00063 C * beginning.) * 00064 C * RepeatedSteps: Counter for the number of failed and thus * 00065 C * repeated integration steps. (Note that * 00066 C * RepeatedSteps isn't set to zero at the * 00067 C * beginning.) * 00068 C * ErrorCode: Error code having the following meaning: * 00069 C * 0: everything was o. k., * 00070 C * 1: illegal integration order, * 00071 C * 2: stepsize underflow. * 00072 C * * 00073 C * Other parameters: * 00074 C * AuxPsi: Auxiliary array of (minimum) size * 00075 C * PsiDim*(IntOrder+1). * 00076 C * CData: Complex*16 data needed in the Func (and possibly * 00077 C * CalcError) routine. * 00078 C * RData: Real*8 data needed in the Func (and possibly * 00079 C * CalcError) routine. * 00080 C * IData: Integer data needed in the Func (and possibly * 00081 C * CalcError) routine. * 00082 C * LData: Logical data needed in the Func (and possibly * 00083 C * CalcError) routine. * 00084 C * * 00085 C * External routines: * 00086 C * Func: Computes the time derivative DtPsi of Psi at time * 00087 C * AbsTime. Called as Func(AbsTime,Psi,DtPsi,CData, * 00088 C * RData,IData,LData). * 00089 C * CalcError: Determines the error of the current integration * 00090 C * step. Called as CalcError(Predicted_Psi, * 00091 C * Corrected_Psi,PsiDim,Error,CData,RData,IData, * 00092 C * LData). * 00093 C * WriteStep: Writes the stepsize and error to a file. * 00094 C * Called as WriteStep(Step,Order,Stepsize,Error). * 00095 C * * 00096 C * V6.0 MB * 00097 C * V7.0 GW addition of CData,RData,IData,LData arrays * 00098 C * * 00099 C ********************************************************************** 00100 00103 Subroutine ABM (Psi,DtPsi,PsiDim,IntPeriod,AbsTime,IntOrder, 00104 + InitStep,TolError,RestartABM,Steps,RepeatedSteps, 00105 + ErrorCode,AuxPsi,Func,CalcError,CData, 00106 + RData,IData,LData) 00107 00108 Implicit None 00109 00110 Real*8 One6th,One30th,One210th,RelativeMinStep 00111 Integer MaxOrder 00112 Parameter (One6th = 1.0D0/6.0D0,One30th = 1.0D0/30.0D0, 00113 + One210th = 1.0D0/210.0D0,RelativeMinStep = 1.0D-12, 00114 + MaxOrder = 8) 00115 00116 Logical RestartABM,LData(*) 00117 Integer PsiDim,IntOrder,Steps,RepeatedSteps,ErrorCode,IData(*) 00118 Real*8 IntPeriod,AbsTime,InitStep,TolError,RData(*) 00119 Complex*16 Psi(PsiDim),DtPsi(PsiDim),AuxPsi(PsiDim,IntOrder+1), 00120 + CData(*) 00121 External Func,CalcError 00122 00123 Logical StepIsRepeated,CurOrdEqIntOrd,UseOldH1 00124 Integer CurrentOrder,D,P,I 00125 Real*8 Time,NextTime,IntError,MinError,Error,D2,D3,D4,D5,D6, 00126 + D7,H1,H1Sqr,MinusH1Cube,MinStepSize,Distance(MaxOrder), 00127 + H(MaxOrder-1),PredCoef(MaxOrder),CorrCoef(MaxOrder), 00128 + SumOfD(2:MaxOrder-2,2:MaxOrder-1), 00129 + ProdOfD(2:MaxOrder-2,2:MaxOrder-1) 00130 Complex*16 InterimAuxPsi(MaxOrder) 00131 00132 Save IntError,Error,MinError,H,Distance,CurrentOrder, 00133 + CurOrdEqIntOrd,StepIsRepeated,UseOldH1 00134 00135 00136 C --- CHECK INTEGRATION ORDER --- 00137 00138 If (IntOrder .GT. MaxOrder) Then 00139 ErrorCode = 1 00140 Return 00141 EndIf 00142 00143 C --- INITIALISE VARIABLES --- 00144 00145 ErrorCode = 0 00146 Time = AbsTime 00147 MinStepSize = RelativeMinStep*IntPeriod 00148 00149 C --- CONTINUE INTEGRATION --- 00150 00151 If (.Not. RestartABM) Then 00152 Goto 200 00153 EndIf 00154 00155 C --- INITIALISE VARIABLES --- 00156 00157 C "UseOldH1" is a flag set to true if the step size is artificially 00158 C shortened in order to fit it to the remaining integration period 00159 C before leaving the routine. This allows the next stepsize to be 00160 C corrected back to the previous value. 00161 00162 IntError = 0.40D0*TolError 00163 MinError = 0.01D0*TolError 00164 Do P = 1,IntOrder-1 00165 H(P) = InitStep 00166 EndDo 00167 CurrentOrder = 2 00168 CurOrdEqIntOrd = .False. 00169 StepIsRepeated = .False. 00170 UseOldH1 = .False. 00171 00172 C --- INITIALISE AUXPSI --- 00173 00174 Do D = 1,PsiDim 00175 AuxPsi(D,1) = Psi(D) 00176 AuxPsi(D,2) = DtPsi(D) 00177 EndDo 00178 00179 C --- INTEGRATION LOOP --- 00180 00181 100 Continue 00182 00183 C --- CHECK WETHER STEPSIZE IS TOO SMALL --- 00184 00185 If ((H(1) .LT. MinStepSize) .And. 00186 + (H(1) .LT. AbsTime+IntPeriod-Time)) Then 00187 ErrorCode = 2 00188 Return 00189 EndIf 00190 00191 C --- CALCULATE PREDICTOR AND CORRECTOR COEFFICIENTS --- 00192 00193 Do P = 1,CurrentOrder 00194 If (P .Eq. 1) Then 00195 H1 = H(1) 00196 Distance(1) = H1 00197 PredCoef(1) = H1 00198 CorrCoef(1) = H1 00199 ElseIf (P .Eq. 2) Then 00200 H1Sqr = H1*H1 00201 Distance(2) = Distance(1)+H(2) 00202 PredCoef(2) = 0.5*H1Sqr 00203 CorrCoef(2) = -PredCoef(2) 00204 ElseIf (P .Eq. 3) Then 00205 MinusH1Cube = -H1Sqr*H1 00206 D2 = H(2) 00207 Distance(3) = Distance(2)+H(3) 00208 PredCoef(3) = One6th*H1Sqr*(3.0*D2+2.0*H1) 00209 CorrCoef(3) = One6th*MinusH1Cube 00210 ElseIf (P .Eq. 4) Then 00211 D3 = D2+H(3) 00212 Distance(4) = Distance(3)+H(4) 00213 SumOfD(2,3) = D2+D3 00214 ProdOfD(2,3) = D2*D3 00215 PredCoef(4) = One6th*H1Sqr*(3.0*ProdOfD(2,3) 00216 + +(2.0*SumOfD(2,3)+1.5*H1)*H1) 00217 CorrCoef(4) = One6th*MinusH1Cube*(D2+0.5*H1) 00218 ElseIf (P .Eq. 5) Then 00219 D4 = D3+H(4) 00220 Distance(5) = Distance(4)+H(5) 00221 SumOfD(3,3) = D3 00222 ProdOfD(3,3) = D3 00223 Do I = 2,3 00224 SumOfD(I,4) = SumOfD(I,3)+D4 00225 ProdOfD(I,4) = ProdOfD(I,3)*D4 00226 EndDo 00227 PredCoef(5) = One30th*H1Sqr*(15.0*ProdOfD(2,4) 00228 + +(10.0*(D2*SumOfD(3,4)+ProdOfD(3,4)) 00229 + +(7.5*SumOfD(2,4)+6.0*H1)*H1)*H1) 00230 CorrCoef(5) = One30th*MinusH1Cube*(5.0*ProdOfD(2,3) 00231 + +(2.5*SumOfD(2,3)+1.5*H1)*H1) 00232 ElseIf (P .Eq. 6) Then 00233 D5 = D4+H(5) 00234 Distance(6) = Distance(5)+H(6) 00235 SumOfD(4,4) = D4 00236 ProdOfD(4,4) = D4 00237 Do I = 2,4 00238 SumOfD(I,5) = SumOfD(I,4)+D5 00239 ProdOfD(I,5) = ProdOfD(I,4)*D5 00240 EndDo 00241 PredCoef(6) = One30th*H1Sqr*(15.0*ProdOfD(2,5) 00242 + +(10.0*(D2*(D3*SumOfD(4,5)+ProdOfD(4,5)) 00243 + +ProdOfD(3,5)) 00244 + +(7.5*(D2*SumOfD(3,5)+D3*SumOfD(4,5) 00245 + +ProdOfD(4,5)) 00246 + +(6.0*SumOfD(2,5)+5.0*H1)*H1)*H1)*H1) 00247 CorrCoef(6) = One30th*MinusH1Cube*(5.0*ProdOfD(2,4) 00248 + +(2.5*(D2*SumOfD(3,4)+ProdOfD(3,4)) 00249 + +(1.5*SumOfD(2,4)+H1)*H1)*H1) 00250 ElseIf (P .Eq. 7) Then 00251 D6 = D5+H(6) 00252 Distance(7) = Distance(6)+H(7) 00253 SumOfD(5,5) = D5 00254 ProdOfD(5,5) = D5 00255 Do I = 2,5 00256 SumOfD(I,6) = SumOfD(I,5)+D6 00257 ProdOfD(I,6) = ProdOfD(I,5)*D6 00258 EndDo 00259 PredCoef(7) = One210th*H1Sqr*(105.0*ProdOfD(2,6) 00260 + +(70.0*(D2*(D3*(D4*SumOfD(5,6)+ProdOfD(5,6)) 00261 + +ProdOfD(4,6))+ProdOfD(3,6)) 00262 + +(52.5*(D2*(D3*SumOfD(4,6)+D4*SumOfD(5,6) 00263 + +ProdOfD(5,6))+D3*(D4*SumOfD(5,6) 00264 + +ProdOfD(5,6))+ProdOfD(4,6)) 00265 + +(42.0*(D2*SumOfD(3,6)+D3*SumOfD(4,6) 00266 + +D4*SumOfD(5,6)+ProdOfD(5,6)) 00267 + +(35.0*SumOfD(2,6)+30.0*H1)*H1)*H1)*H1)*H1) 00268 CorrCoef(7) = One210th*MinusH1Cube*(35.0*ProdOfD(2,5) 00269 + +(17.5*(D2*(D3*SumOfD(4,5)+ProdOfD(4,5)) 00270 + +ProdOfD(3,5)) 00271 + +(10.5*(D2*SumOfD(3,5)+D3*SumOfD(4,5) 00272 + +ProdOfD(4,5))+(7.0*SumOfD(2,5) 00273 + +5.0*H1)*H1)*H1)*H1) 00274 ElseIf (P .Eq. 8) Then 00275 D7 = D6+H(7) 00276 SumOfD(6,6) = D6 00277 ProdOfD(6,6) = D6 00278 Do I = 2,6 00279 SumOfD(I,7) = SumOfD(I,6)+D7 00280 ProdOfD(I,7) = ProdOfD(I,6)*D7 00281 EndDo 00282 PredCoef(8) = One210th*H1Sqr*(105.0*ProdOfD(2,7) 00283 + +(60.0*(D2*(D3*(D4*(D5*SumOfD(6,7) 00284 + +ProdOfD(6,7))+ProdOfD(5,7))+ProdOfD(4,7)) 00285 + +ProdOfD(3,7)) 00286 + +(52.5*(D2*(D3*(D4*SumOfD(5,7)+D5*SumOfD(6,7)) 00287 + +ProdOfD(6,7))+D4*(D5*SumOfD(6,7) 00288 + +ProdOfD(6,7))+ProdOfD(5,7)) 00289 + +D3*(D4*(D5*SumOfD(6,7)+ProdOfD(6,7) 00290 + +ProdOfD(5,7))+ProdOfD(4,7)) 00291 + +(42.0*(D2*(D3*SumOfD(4,7)+D4*SumOfD(5,7) 00292 + +D5*SumOfD(6,7)+ProdOfD(6,7)) 00293 + +D3*(D4*SumOfD(5,7)+D5*SumOfD(6,7) 00294 + +ProdOfD(6,7))+D4*(D5*SumOfD(6,7) 00295 + +ProdOfD(6,7))+ProdOfD(5,7)) 00296 + +(35.0*(D2*SumOfD(3,7)+D3*SumOfD(4,7) 00297 + +D4*SumOfD(5,7)+D5*SumOfD(6,7)+ProdOfD(6,7)) 00298 + +(30.0*SumOfD(2,7) 00299 + +26.25*H1)*H1)*H1)*H1)*H1)*H1) 00300 CorrCoef(8) = One210th*MinusH1Cube*(35.0*ProdOfD(2,6) 00301 + +(17.5*(D2*(D3*(D4*SumOfD(5,6)+ProdOfD(5,6)) 00302 + +ProdOfD(4,6))+ProdOfD(3,6)) 00303 + +(10.5*(D2*(D3*SumOfD(4,6)+D4*SumOfD(5,6) 00304 + +ProdOfD(5,6))+D3*(D4*SumOfD(5,6) 00305 + +ProdOfD(5,6))+ProdOfD(4,6)) 00306 + +(7.0*(D2*SumOfD(3,6)+D3*SumOfD(4,6) 00307 + +D4*SumOfD(5,6)+ProdOfD(5,6))+(5.0*SumOfD(2,6) 00308 + +3.75*H1)*H1)*H1)*H1)*H1) 00309 Else 00310 ErrorCode = 1 00311 Return 00312 EndIf 00313 EndDo 00314 00315 C --- RESTORE ORIGINAL PSI WHEN STEP IS REPEATED --- 00316 00317 If (StepIsRepeated) Then 00318 Do D = 1,PsiDim 00319 Psi(D) = AuxPsi(D,1) 00320 EndDo 00321 StepIsRepeated = .False. 00322 EndIf 00323 00324 C --- PREDICT PSI --- 00325 00326 C As long as the order p is increased, the predictor must be used as a 00327 C p-step method. Later on, the predictor is run as a (p+1)-step method. 00328 00329 Do D = 1,PsiDim 00330 If (CurOrdEqIntOrd) Then 00331 Do P = 1,CurrentOrder 00332 Psi(D) = Psi(D)+PredCoef(P)*AuxPsi(D,P+1) 00333 EndDo 00334 Else 00335 Do P = 1,CurrentOrder-1 00336 Psi(D) = Psi(D)+PredCoef(P)*AuxPsi(D,P+1) 00337 EndDo 00338 EndIf 00339 EndDo 00340 00341 C --- EVALUATE FUNCTION WITH PREDICTED PSI --- 00342 00343 NextTime = Time+H(1) 00344 Call Func(NextTime,Psi,DtPsi,CData,RData,IData,LData) 00345 00346 C --- SAVE PREDICTED PSI IN DTPSI AND CORRECT PSI --- 00347 00348 Do D = 1,PsiDim 00349 InterimAuxPsi(1) = DtPsi(D) 00350 Do P = 2,CurrentOrder 00351 InterimAuxPsi(P) = (InterimAuxPsi(P-1)-AuxPsi(D,P)) 00352 + /Distance(P-1) 00353 EndDo 00354 DtPsi(D) = Psi(D) 00355 Psi(D) = AuxPsi(D,1) 00356 Do P = 1,CurrentOrder 00357 Psi(D) = Psi(D)+CorrCoef(P)*InterimAuxPsi(P) 00358 EndDo 00359 EndDo 00360 00361 C --- CALCULATE ERROR AND (IF REQUIRED) NEW STEP SIZE --- 00362 00363 Call CalcError(DtPsi,Psi,PsiDim,Error,CData,RData,IData,LData) 00364 If (CurrentOrder .LT. IntOrder) Then 00365 Error = Error*(2*IntOrder+1.0D0)/(2*CurrentOrder+1.0D0) 00366 EndIf 00367 Steps = Steps+1 00368 If (Error .GT. TolError) Then 00369 H(1) = 0.8D0*H(1)*(IntError/Error)**(1.0D0/(CurrentOrder+1)) 00370 RepeatedSteps = RepeatedSteps+1 00371 StepIsRepeated = .True. 00372 Goto 100 00373 EndIf 00374 00375 C --- WRITE STEPSIZE AND ERROR --- 00376 00377 ! Call WriteStep(Steps,CurrentOrder,H(1),Error,Time) 00378 00379 C --- RETURN WHEN FINISHED --- 00380 00381 Time = Time+H(1) 00382 If (Time .GE. AbsTime+IntPeriod-MinStepSize) Then 00383 Return 00384 EndIf 00385 00386 C --- EVALUATE FUNCTION WITH CORRECTED PSI --- 00387 00388 Call Func(NextTime,Psi,DtPsi,CData,RData,IData,LData) 00389 00390 200 Continue 00391 00392 C --- INCREASE INTEGRATION ORDER --- 00393 00394 If (CurrentOrder .LT. IntOrder) Then 00395 CurrentOrder = CurrentOrder+1 00396 Else 00397 CurOrdEqIntOrd = .True. 00398 EndIf 00399 00400 C --- COMPUTE NEW DEVIDED DIFFERENCES --- 00401 00402 C As long as the order p is increased, the predictor must be used as a 00403 C p-step method. Later on, the predictor is run as a (p+1)-step method. 00404 00405 If (CurOrdEqIntOrd) Then 00406 Do D = 1,PsiDim 00407 AuxPsi(D,1) = Psi(D) 00408 InterimAuxPsi(1) = DtPsi(D) 00409 Do P = 2,CurrentOrder 00410 InterimAuxPsi(P) = (InterimAuxPsi(P-1)-AuxPsi(D,P)) 00411 + /Distance(P-1) 00412 EndDo 00413 Do P = 2,CurrentOrder+1 00414 AuxPsi(D,P) = InterimAuxPsi(P-1) 00415 EndDo 00416 EndDo 00417 Else 00418 Do D = 1,PsiDim 00419 AuxPsi(D,1) = Psi(D) 00420 InterimAuxPsi(1) = DtPsi(D) 00421 Do P = 2,CurrentOrder-1 00422 InterimAuxPsi(P) = (InterimAuxPsi(P-1)-AuxPsi(D,P)) 00423 + /Distance(P-1) 00424 EndDo 00425 Do P = 2,CurrentOrder 00426 AuxPsi(D,P) = InterimAuxPsi(P-1) 00427 EndDo 00428 EndDo 00429 EndIf 00430 00431 C --- CALCULATE NEW STEP SIZE --- 00432 00433 C If the previous step had been shortened (see above), the next stepsize 00434 C may be set (close) to the stepsize before that one. 00435 00436 Error = Max(Error,MinError) 00437 Do P = CurrentOrder-1,2,-1 00438 H(P) = H(P-1) 00439 EndDo 00440 H(1) = H(1)*(IntError/Error)**(1.0D0/(CurrentOrder+1)) 00441 If (UseOldH1) Then 00442 If (H(1) .LT. 0.98D0*H(3)) Then 00443 H(1) = 0.98D0*H(3) 00444 EndIf 00445 UseOldH1 = .False. 00446 EndIf 00447 If (Time+H(1) .GT. AbsTime+IntPeriod) Then 00448 H(1) = AbsTime+IntPeriod-Time 00449 UseOldH1 = .True. 00450 EndIf 00451 00452 C --- CONTINUE INTEGRATION --- 00453 00454 Goto 100 00455 00456 End Subroutine ABM 00457 00458 C ********************************************************************** 00459 C * * 00460 C * SUBROUTINE ABSABMERROR * 00461 C * * 00462 C * Computes the absolute error of the current integration step of the * 00463 C * Adams-Bashforth-Moulton predictor-corrector integrator. * 00464 C * * 00465 C * Input parameters: * 00466 C * Predicted: Predicted solution vector. * 00467 C * Corrected: Corrected soultion vector. * 00468 C * PsiDim: Length of solution vectors. * 00469 C * Cdata: Complex data array (not used in this error routine) * 00470 C * * 00471 C * Output parameters: * 00472 C * Error: Computed error. * 00473 C * * 00474 C * Other parameters: * 00475 C * CData: Complex data array (not used in this error routine). * 00476 C * RData: Real data array (not used in this error routine). * 00477 C * IData: Integer data array (not used in this error routine). * 00478 C * LData: Logical data array (not used in this error routine). * 00479 C * * 00480 C * V6.0 MB * 00481 C * V7.0 GW addition of Data arrays * 00482 C * * 00483 C ********************************************************************** 00484 00487 Subroutine AbsABMError (Predicted,Corrected,Psidim,Error,CData, 00488 + RData,IData,LData) 00489 00490 Implicit None 00491 00492 Logical LData(*) 00493 Integer PsiDim,IData(*) 00494 Real*8 Error,RData(*) 00495 Complex*16 Predicted(PsiDim),Corrected(PsiDim),CData(*) 00496 00497 Integer D 00498 Real*8 SingleError,Boundary 00499 00500 C --- INITIALISE VARIABLES --- 00501 00502 C If a component of the corrected solution is less than "Boundary" the 00503 C absolute error is taken, else the relative error is computed. 00504 00505 Boundary = 1.0D0 00506 Error = 0.0D0 00507 00508 C --- COMPUTE ERROR OF BOTH REAL AND IMAGINARY PART --- 00509 00510 Do D = 1,PsiDim 00511 SingleError = DAbs(Dble(Predicted(D))-Dble(Corrected(D))) 00512 If (DAbs(Dble(Corrected(D))) .GT. Boundary) Then 00513 SingleError = SingleError/DAbs(Dble(Corrected(D))) 00514 EndIf 00515 Error = Max(Error,SingleError) 00516 SingleError = DAbs(DImag(Predicted(D))-DImag(Corrected(D))) 00517 If (DAbs(DImag(Corrected(D))) .GT. Boundary) Then 00518 SingleError = SingleError/DAbs(DImag(Corrected(D))) 00519 EndIf 00520 Error = Max(Error,SingleError) 00521 EndDo 00522 00523 Return 00524 End Subroutine AbsABMError 00525 00526 C ********************************************************************** 00527 C * * 00528 C * SUBROUTINE RELABMERROR * 00529 C * * 00530 C * Computes the relative error of the current integration step of the * 00531 C * Adams-Bashforth-Moulton predictor-corrector integrator. * 00532 C * * 00533 C * Input parameters: * 00534 C * Predicted: Predicted solution vector. * 00535 C * Corrected: Corrected soultion vector. * 00536 C * PsiDim: Length of solution vectors. * 00537 C * * 00538 C * Output parameters: * 00539 C * Error: Computed error. * 00540 C * * 00541 C * Other parameters: * 00542 C * CData: Complex data array (not used in this error routine). * 00543 C * RData: Real data array (not used in this error routine). * 00544 C * IData: Integer data array (not used in this error routine). * 00545 C * LData: Logical data array (not used in this error routine). * 00546 C * * 00547 C * V6.0 MB * 00548 C * V7.0 MB addition of Data arrays * 00549 C * * 00550 C ********************************************************************** 00551 00554 Subroutine RelABMError (Predicted,Corrected,Psidim,Error,CData, 00555 + RData,IData,LData) 00556 00557 Implicit None 00558 00559 Real*8 Tiny 00560 Parameter (Tiny = 1.0D-30) 00561 00562 Logical LData(*) 00563 Integer PsiDim,IData(*) 00564 Real*8 Error,RData(*) 00565 Complex*16 Predicted(PsiDim),Corrected(PsiDim),CData(*) 00566 00567 Integer D 00568 Real*8 Maximum 00569 00570 C --- COMPUTE ERROR OF BOTH REAL AND IMAGINARY PART --- 00571 00572 Error = 0.0D0 00573 Maximum = Tiny 00574 Do D = 1,PsiDim 00575 Maximum = Max(Maximum,DAbs(Dble(Corrected(D))), 00576 + DAbs(DImag(Corrected(D)))) 00577 Error = Max(Error,DAbs(Dble(Predicted(D))-Dble(Corrected(D))), 00578 + DAbs(DImag(Predicted(D))-DImag(Corrected(D)))) 00579 EndDo 00580 Error = Error/Maximum 00581 00582 Return 00583 End Subroutine RelABMError 00584 00585 C ********************************************************************** 00586 C * * 00587 C * SUBROUTINE WRITEABMSTEP * 00588 C * * 00589 C * Writes for the Adams-Bashforth-Moulton predictor-corrector * 00590 C * integrator the stepsize and error of the current integration step * 00591 C * to a file. (In the current form WriteABMStep is a dummy routine * 00592 C * doing absolutely nothing; it is included here for formal reasons, * 00593 C * but can (if desired) be extended easily such that it writes the * 00594 C * size and error of the current integration step to a file. * 00595 C * * 00596 C * Input parameters: * 00597 C * Step: Counter for the number of steps. * 00598 C * Order: Current integration order. * 00599 C * Stepsize: Size of current integration step. * 00600 C * Error: Error of current integration step. * 00601 C * * 00602 C * Output parameters: * 00603 C * none * 00604 C * * 00605 C * V6.0 MB * 00606 C * * 00607 C ********************************************************************** 00608 00611 Subroutine WriteABMStep (Step,Order,Stepsize,Error) 00612 00613 Implicit None 00614 00615 Integer Step,Order 00616 Real*8 Stepsize,Error 00617 00618 C --- WRITE STEPSIZE AND ERROR --- 00619 00620 Return 00621 End Subroutine WriteABMStep 00622 00623 C ********************************************************************** 00624 C * * 00625 C * SUBROUTINE ABMERRORMSG * 00626 C * * 00627 C * Generates for a given error number returned by "ABM" a * 00628 C * corresponding error message. * 00629 C * * 00630 C * Input parameters: * 00631 C * Error: Error code returned by ABM. * 00632 C * * 00633 C * Output parameters: * 00634 C * Msg: Error message. * 00635 C * * 00636 C * V7.0 MB * 00637 C * * 00638 C ********************************************************************** 00639 00642 Subroutine ABMErrorMsg (Error,Msg) 00643 00644 Implicit None 00645 00646 Integer Error 00647 Character*(*) Msg 00648 00649 C --- GENERATE ERROR MESSAGE --- 00650 00651 If (Error .Eq. 1) Then 00652 Msg = 'Illegal integration order' 00653 ElseIf (Error .Eq. 2) Then 00654 Msg = 'Stepsize underflow' 00655 Else 00656 Msg = 'Unknown error occurred' 00657 EndIf 00658 00659 Return 00660 End Subroutine ABMErrorMsg 00661
1.7.3
