193 lines
6.5 KiB
C++
193 lines
6.5 KiB
C++
|
|
c===============================================================================
|
|
c RUNGE_KUTTA.INC
|
|
c===============================================================================
|
|
|
|
c Dieses Includefile erledigt fuer die Subroutinen 'INTEGRATIONSSTEP_RUNGE_KUTTA'
|
|
c die Fehlerbetrachtung, das Ertasten des Uebergabebereiches zur naechsten Mappe,
|
|
c die damit verbundenen Variationen des Zeitschrittes dt sowie die letztendliche
|
|
c Festlegung des neuen Ortes, der neuen Geschwindigkeit und der neuen Zeit.
|
|
|
|
|
|
c Zaehle die Schritte:
|
|
|
|
steps = steps + 1
|
|
|
|
|
|
c Fehlerbetrachtung:
|
|
|
|
c Fehlerbetrachtung:
|
|
c der groesste (absolute bzw. relative) Fehler im Ort soll kleiner als eps_x
|
|
c sein, der groesste Fehler in der Geschwindigkeit kleiner als eps_v:
|
|
c -> Bestimme den jeweils groessten Fehler der drei Komponenten des Ortes und
|
|
c der Geschwindigkeit (dh. die groesste Differenz der Aederungen):
|
|
|
|
maxErr_x = 0.
|
|
maxErr_v = 0.
|
|
|
|
do i = 1, 3
|
|
xDifferenz(i) = Dx1(i)-Dx2(i)
|
|
vDifferenz(i) = Dv1(i)-Dv2(i)
|
|
if (log_relativ) then
|
|
if (Dx1(i).NE.0.) maxErr_x = Max(maxErr_x,Abs(xDifferenz(i)/Dx1(i)))
|
|
if (Dv1(i).NE.0.) maxErr_v = Max(maxErr_v,Abs(vDifferenz(i)/Dv1(i)))
|
|
else
|
|
maxErr_x = Max( maxErr_x, Abs( xDifferenz(i) ) )
|
|
maxErr_v = Max( maxErr_v, Abs( vDifferenz(i) ) )
|
|
endif
|
|
enddo
|
|
|
|
c - Skaliere den jeweils groessten relativen Fehler auf das jeweilige Epsilon:
|
|
|
|
maxErr_x = maxErr_x / eps_x
|
|
maxErr_v = maxErr_v / eps_v
|
|
|
|
c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
|
|
c der groessere der beiden reskalierten Fehler bestimmt, ob der Integrations-
|
|
c schritt mit kleinerem Zeitintervall wiederholt werden muss, bzw. um welchen
|
|
c Faktor das Zeitintervall fuer den naechsten Schritt vergroessert werden kann:
|
|
|
|
c Liegt der Fehler ausserhalb des Toleranzbereiches und ist dt bereits jetzt
|
|
c kleiner als dtsmall, so mache keinen neuen Versuch sondern akzeptiere als Not-
|
|
c loesung den bestehenden Naeherungswert. Setze dt in diesem Fall als Default
|
|
c fuer den kommenden Integrationsschritt auf dtsmall. Setze aber auch das flag
|
|
c 'flag_dtsmall', damit gezaehlt werden kann, wie oft dieses Prozedur fuer ein
|
|
c bestimmtes Teilchen angewendet werden muss. Ist dies zu oft der Fall, so brich
|
|
c diese Trajektorienberechnung ganz ab (-> destiny = code_dtsmall).
|
|
c (2. Teil erfolgt weiter unten)
|
|
c
|
|
c Es kam vor, dass ohne Ueberschreitung der Fehlertoleranz ein 'dtupper'
|
|
c und ein 'dtlower' gefunden wurde, beim Ertasten des Uebergabebereiches
|
|
c die Fehlergrenze bei mittleren dt-Werten dann aber ueberschritten wurde,
|
|
c wodurch dt immer wieder verkuerzt wurde, ohne dass der Uebergabebereich
|
|
c erreicht werden konnte. Letztlich bildete das ganze eine unendliche Schleife.
|
|
c Daher werden jetzt jedesmal, wenn die Fehlergrenze ueberschritten wird
|
|
c 'found_upper' und 'found_lower' resettet.
|
|
|
|
|
|
maxErr = Max(maxErr_x,maxErr_v)
|
|
|
|
if (maxErr.GT.1.) then
|
|
found_upper_upp = .false.
|
|
found_lower_upp = .false.
|
|
found_upper_low = .false.
|
|
found_lower_low = .false.
|
|
if (dt.LT.dtsmall) then ! Fehler immer noch zu gross, obwohl
|
|
flag_dtsmall = .true. ! dtsmall schon unterschritten ist
|
|
else
|
|
!c Bestimme kuerzeren Zeitschritt fuer neuen Versuch (vgl. Referenz):
|
|
dt = safety * dt * (maxErr**pShrink)
|
|
goto 10
|
|
endif
|
|
endif
|
|
|
|
|
|
x_1 = x(1) + Dx1(1) + xDifferenz(1) / 15.
|
|
|
|
|
|
c Falls x(1) (== x_1) jetzt jenseits des Mappenendes liegen sollte, behalte
|
|
c dieses Faktum im Gedaechtnis und verkuerze den aktuell verwendeten Zeitschritt
|
|
c so lange um Faktor 0.5, bis x(1) innerhalb oder vor dem Uebergabebereich liegt.
|
|
c Liegt es dann davor, suche einen mittleren Zeitschritt, bei dem es innerhalb
|
|
c liegt.
|
|
c Hat das Teilchen danach (oder nachdem es direkt in den Uebergabebereich traf)
|
|
c positives v(1), so setze das Logical 'reachedEndOfMap' fuer die Berechnung
|
|
c des Schnittpunkts der Trajektorie mit dem Mappenende.
|
|
c (v(1)<0. ist entweder moeglich falls es bereits vor dem Mappenende reflektiert
|
|
c wurde oder gerade aus Mappe mit hoeherer Nummer kam).
|
|
|
|
if (x_1.GT.xStartUeberUpp) then
|
|
if (.NOT.found_upper_upp) dt_save = dt
|
|
if (x_1.LE.xMax .AND. v(1).GT.0.) then
|
|
reachedEndOfMap = .true.
|
|
elseif (x_1.GT.xMax) then
|
|
dtupper = dt
|
|
found_upper_upp = .true.
|
|
if (.NOT.found_lower_upp) then
|
|
dt = min(0.5*dt,(xStartUeberUpp-x(1))/(x_1-x(1))*dt)
|
|
else
|
|
dt = (dtlower+dtupper)/2.
|
|
endif
|
|
goto 10 ! neue Berechnung
|
|
endif
|
|
elseif (found_upper_upp) then
|
|
found_lower_upp = .true.
|
|
dtlower = dt
|
|
dt = (dtlower+dtupper)/2.
|
|
goto 10 ! neue Berechnung
|
|
endif
|
|
|
|
|
|
c entsprechende Behandlung wie oben fuer den Fall, dass x(1) (== x_1) jetzt im
|
|
c Bereich des Mappenanfangs liegt:
|
|
|
|
if (x_1.LT.xStartUeberLow) then
|
|
if (.NOT.found_upper_low) dt_save = dt
|
|
if (x_1.GE.xMin .AND. v(1).LT.0) then
|
|
backOneMap = .true.
|
|
elseif (x_1.LT.xmin) then
|
|
found_upper_low = .true.
|
|
dtupper = dt
|
|
if (.NOT.found_lower_low) then
|
|
dt = min(0.5*dt,(xStartUeberLow-x(1))/(x_1-x(1))*dt)
|
|
else
|
|
dt = (dtlower+dtupper)/2.
|
|
endif
|
|
goto 10 ! neue Berechnung
|
|
endif
|
|
elseif (found_upper_low) then
|
|
found_lower_low = .true.
|
|
dtlower = dt
|
|
dt = (dtlower+dtupper)/2.
|
|
goto 10 ! neue Berechnung
|
|
endif
|
|
|
|
|
|
c Nimm die Ergebnisse aus dem dt-Schritt und den beiden dt/2-Schritten und
|
|
c berechne damit den neuen Ort und die neue Geschwindigkeit mit Genauigkeit
|
|
c fuenfter Ordnung in dt:
|
|
|
|
x(1) = x_1
|
|
x(2) = x(2) + Dx1(2) + xDifferenz(2) / 15.
|
|
x(3) = x(3) + Dx1(3) + xDifferenz(3) / 15.
|
|
|
|
v(1) = v(1) + Dv1(1) + vDifferenz(1) / 15.
|
|
v(2) = v(2) + Dv1(2) + vDifferenz(2) / 15.
|
|
v(3) = v(3) + Dv1(3) + vDifferenz(3) / 15.
|
|
|
|
|
|
c alten Zeitschritt addieren:
|
|
|
|
t = t + dt
|
|
|
|
|
|
c Falls ein Uebergabebereich erreicht wurde, berechne Schnittpunkt der
|
|
c Trajektorie mit x=xmax (Mappenende) bzw. mit x=xmin (Mappenanfang):
|
|
|
|
if (reachedEndOfMap) goto 7766
|
|
if (backOneMap) goto 7767
|
|
|
|
|
|
c neuen Zeitschritt so gross wie sinnvoller weise moeglich machen:
|
|
|
|
3454 if (flag_dtSmall) then
|
|
if (n_dtsmall.LT.maxBelowDtSmall) then
|
|
dt = dtSmall ! fuer naechsten RK-Schritt
|
|
n_dtsmall = n_dtsmall + 1
|
|
else
|
|
destiny = code_dtsmall ! gib Teilchen verloren
|
|
RETURN
|
|
endif
|
|
else
|
|
if (maxErr.GT.errCon) then
|
|
dt = safety * dt * (maxErr**pGrow) ! vgl. Referenz
|
|
else
|
|
dt = 4. * dt ! <- Vergroesserung des Zeitschritts max. um
|
|
endif ! Faktor 4!
|
|
! pruefen, ob Maximallaenge fuer ersten Testschritt nicht ueberschritten ist:
|
|
if (log_confine) dt = min(dt,dl_max/sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)))
|
|
endif
|
|
|
|
|