summaryrefslogtreecommitdiff
path: root/Physikunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'Physikunit.pas')
-rw-r--r--Physikunit.pas520
1 files changed, 376 insertions, 144 deletions
diff --git a/Physikunit.pas b/Physikunit.pas
index 8c216c9..da9d56e 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -2,16 +2,20 @@ unit Physikunit;
{$mode objfpc}{$H+}
+{$IFNDEF UNIX}
+ {$ERROR This program can be compiled only on/for Unix/Linux based systems.}
+{$ENDIF}
+
{$DEFINE Zeitschrittueberwachung}
{$DEFINE Dichteueberwachung}
interface
uses
- Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit;
+ Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix;
type
- tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier);
+ tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier);
tVerteilungsfunktion = function(x: extended): extended;
tEMFeldInhalt = (
efA,efAX,efAY,efAZ,
@@ -54,7 +58,7 @@ type
ortsGrenzen: array of extended;
dichteStuecke: array of string;
public
- nMax, // maximale Massendichte in nc
+ nMin,nMax, // minimale und maximale Massendichte in nc
spezifischeLadung: extended; // q/m in e/me
constructor create;
destructor destroy; override;
@@ -66,9 +70,9 @@ type
{ tWertePunkt }
tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen
- matWerte: array of array[tMaterieFeldInhalt,boolean] of double;
+ matWerte: array of array[tMaterieFeldInhalt,boolean] of extended;
// Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen
- emWerte: array[tEMFeldInhalt,boolean] of double;
+ emWerte: array[tEMFeldInhalt,boolean] of extended;
// EM-Felder und deren Zeitableitungen
// A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
lN,rN: tWertePunkt;
@@ -76,10 +80,13 @@ type
constructor create(linkerNachbar: tWertePunkt; materien: longint);
destructor destroy; override;
procedure berechneGammaUndP;
- procedure berechneNAbleitung(iDX,pDNMax: extended);
+ procedure berechneNAbleitung(iDX,pDNMax: extended); overload;
+ procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload;
procedure berechneAbleitungen(dX,iDX: extended);
- procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt
- procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); overload;
+ procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload;
+ procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload;
function maxDT: extended;
end;
@@ -87,8 +94,9 @@ type
tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
private
- matAnz: longint;
- lichters: array[boolean] of tMyStringlist;
+ matAnz: longint;
+ spezLadungen: tExtendedArray;
+ lichters: array[boolean] of tMyStringlist;
procedure setzeRaender(dT,iDX,iDT: extended);
public
inhalt: array of tWertePunkt;
@@ -97,8 +105,10 @@ type
constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
destructor destroy; override;
procedure berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
- procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt
- procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2,in3: tFelder; fak2,fak3: extended); overload;
+ procedure liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); overload;
+ procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload;
function maxDT: extended;
end;
@@ -112,6 +122,7 @@ type
dX,iDX,pDNMax,dTMaximum,xl,t: extended;
kvs: tKnownValues;
procedure diffusionsTermAnpassen(pro: tProtokollant);
+ function pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert?
public
aktuelleFelder: longint;
@@ -134,11 +145,14 @@ type
fortschrittsAnzeige: boolean;
ausgabeDateien: array of tAusgabeDatei;
public
+ gotSighup: boolean;
constructor create(inName: string; protokollant: tProtokollant; name: string);
destructor destroy; override;
function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende?
end;
+procedure SignalCapture(signal : longint); cdecl;
+
implementation
const
@@ -149,6 +163,9 @@ const
'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA'
);
+var
+ sighupSimulationen: array of tSimulation;
+
{ tAusgabeDatei }
constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
@@ -284,7 +301,7 @@ begin
for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
while cX>=sX do begin
dec(cnt);
- blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(double));
+ blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(extended));
sX:=sX+sDX;
end;
cX:=cX+gitter.dX;
@@ -294,7 +311,7 @@ begin
for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
while cX>=sX do begin
dec(cnt);
- blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(double));
+ blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(extended));
sX:=sX+sDX;
end;
cX:=cX+gitter.dX;
@@ -317,7 +334,8 @@ begin
setlength(ortsGrenzen,0);
fillchar(dichteStuecke,sizeof(dichteStuecke),#0);
setlength(dichteStuecke,0);
- nMax:=0;
+ nMin:=1E-3;
+ nMax:=nMin;
spezifischeLadung:=0;
end;
@@ -350,7 +368,7 @@ begin
while (i<length(ortsGrenzen)) and (ortsGrenzen[i]<x) do
inc(i);
kvs.add('x',x);
- result:=nMax*exprToFloat(false,dichteStuecke[i],kvs,nil);
+ result:=(nMax-nMin)*exprToFloat(false,dichteStuecke[i],kvs,nil)+nMin;
kvs.rem('x');
end;
@@ -479,15 +497,35 @@ begin
for i:=0 to length(matWerte)-1 do
// dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
matWerte[i,mfN,true]:=
- ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false])
- - ln.matWerte[i,mfN,false] * ln.matWerte[i,mfIGamma,false] * (-pDNMax - rN.matWerte[i,mfPX,false])
- - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2;
+ ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false])
+ + lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false])
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2;
// Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung
// und soll die numerische Stabilität bis zu p * dn/dx / n von 2*pDNMax gewährleisten.
// Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit
// alpha = pDNMax * dX
end;
+procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean);
+var
+ i: longint;
+begin
+ if rechts then begin
+ for i:=0 to length(matWerte)-1 do
+ // rechter Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
+ matWerte[i,mfN,true]:=
+ ( lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false])
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax - matWerte[i,mfPX,false])) * iDX/2;
+ end
+ else begin
+ for i:=0 to length(matWerte)-1 do
+ // linker Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
+ matWerte[i,mfN,true]:=
+ ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false])
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax + matWerte[i,mfPX,false])) * iDX/2;
+ end
+end;
+
procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen
var
i: longint;
@@ -558,6 +596,68 @@ begin
matWerte[i,maF,false]:= in1.matWerte[i,maF,false] + fak2 * in2.matWerte[i,maF,true];
end;
+procedure tWertePunkt.liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ emF: tEMFeldInhalt;
+ maF: tMaterieFeldInhalt;
+ i: longint;
+begin
+(* tEMFeldInhalt = (
+ efA,efAX,efAY,efAZ,
+ efDAXDT,efDAYDT,efDAZDT,
+ efDPhiDX
+ ); *)
+ for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
+ emWerte[emF,false]:=
+ in1.emWerte[emF,false]
+ + fak2 * in2.emWerte[emF,true]
+ + fak3 * in3.emWerte[emF,true];
+
+(* tMaterieFeldInhalt = (
+ mfN,mfDPsiDX,
+ mfP,mfPX,mfPY,mfPZ,
+ mfGamma,mfIGamma
+ ); *)
+ for i:=0 to length(matWerte)-1 do // siehe oben
+ for maF:=mfN to mfDPsiDX do
+ matWerte[i,maF,false]:=
+ in1.matWerte[i,maF,false]
+ + fak2 * in2.matWerte[i,maF,true]
+ + fak3 * in3.matWerte[i,maF,true];
+end;
+
+procedure tWertePunkt.liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ emF: tEMFeldInhalt;
+ maF: tMaterieFeldInhalt;
+ i: longint;
+begin
+(* tEMFeldInhalt = (
+ efA,efAX,efAY,efAZ,
+ efDAXDT,efDAYDT,efDAZDT,
+ efDPhiDX
+ ); *)
+ for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
+ emWerte[emF,false]:=
+ in1.emWerte[emF,false]
+ + fak2 * in2.emWerte[emF,true]
+ + fak3 * in3.emWerte[emF,true]
+ + fak4 * in4.emWerte[emF,true];
+
+(* tMaterieFeldInhalt = (
+ mfN,mfDPsiDX,
+ mfP,mfPX,mfPY,mfPZ,
+ mfGamma,mfIGamma
+ ); *)
+ for i:=0 to length(matWerte)-1 do // siehe oben
+ for maF:=mfN to mfDPsiDX do
+ matWerte[i,maF,false]:=
+ in1.matWerte[i,maF,false]
+ + fak2 * in2.matWerte[i,maF,true]
+ + fak3 * in3.matWerte[i,maF,true]
+ + fak4 * in4.matWerte[i,maF,true];
+end;
+
procedure tWertePunkt.liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
var
emF: tEMFeldInhalt;
@@ -596,11 +696,11 @@ function tWertePunkt.maxDT: extended;
var
i: longint;
begin
- result:=1;
+ result:=infinity;
for i:=0 to length(matWerte)-1 do
if (matWerte[i,mfN,true]<0) and
(matWerte[i,mfN,false]>0) then
- result:=min(result,matWerte[i,mfN,false]/matWerte[i,mfN,true]);
+ result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]);
end;
{ tFelder }
@@ -613,6 +713,10 @@ begin
inherited create;
par:=parent;
matAnz:=length(teilchen);
+ fillchar(spezLadungen,sizeof(spezLadungen),#0);
+ setlength(spezLadungen,matAnz);
+ for i:=0 to length(spezLadungen)-1 do
+ spezLadungen[i]:=teilchen[i].spezifischeLadung;
fillchar(inhalt,sizeof(inhalt),#0);
setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen
inhalt[0]:=tWertePunkt.create(nil,matAnz);
@@ -634,7 +738,13 @@ begin
end;
destructor tFelder.destroy;
+var
+ i: longint;
begin
+ setlength(spezLadungen,0);
+ for i:=0 to length(inhalt)-1 do
+ inhalt[i].free;
+ setlength(inhalt,0);
lichters[false].free;
lichters[true].free;
inherited destroy;
@@ -645,7 +755,7 @@ var
emF: tEMFeldInhalt;
rechts: boolean;
i: longint;
- oVal: extended;
+ nVal: extended;
begin
for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld
inhalt[0].emWerte[emF,true]:=
@@ -658,12 +768,12 @@ begin
for rechts:=false to true do
for i:=0 to lichters[rechts].count-1 do begin
- par.kvs.add('t',par.t-dT);
- oVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil);
+ par.kvs.add('t',par.t+dT);
+ nVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil);
par.kvs.add('t',par.t);
inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true]:=
inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true] +
- (exprToFloat(false,lichters[rechts][i],par.kvs,nil) - oVal)*iDT;
+ (nVal - exprToFloat(false,lichters[rechts][i],par.kvs,nil))*iDT;
end;
end;
@@ -680,11 +790,15 @@ begin
-abs(inhalt[length(inhalt)-2].matWerte[i,mfPX,false]);
end;
+ for i:=2 to length(inhalt)-3 do
+ inhalt[i].berechneNAbleitung(iDX,pDNMax);
+
+ inhalt[1].berechneNAbleitung(iDX,pDNMax,false);
+ inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true);
+
setzeRaender(dT,iDX,iDT);
for i:=1 to length(inhalt)-2 do
- inhalt[i].berechneNAbleitung(iDX,pDNMax);
- for i:=1 to length(inhalt)-2 do
inhalt[i].berechneAbleitungen(dX,iDX);
end;
@@ -696,6 +810,22 @@ begin
inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],fak2);
end;
+procedure tFelder.liKo(in1,in2,in3: tFelder; fak2,fak3: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ i: longint;
+begin
+ for i:=0 to length(inhalt)-1 do
+ inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],fak2,fak3);
+end;
+
+procedure tFelder.liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ i: longint;
+begin
+ for i:=0 to length(inhalt)-1 do
+ inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],in4.inhalt[i],fak2,fak3,fak4);
+end;
+
procedure tFelder.liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
var
i: longint;
@@ -732,8 +862,10 @@ begin
end
else pDNMax:=pDNMa;
case Zeitverfahren of
- zfEulerVorwaerts: Setlength(Felders,2);
- zfRungeKuttaVier: Setlength(Felders,5);
+ zfEulerVorwaerts:
+ Setlength(Felders,2);
+ zfRungeKuttaDreiAchtel,zfRungeKuttaVier:
+ Setlength(Felders,5);
end{of Case};
xl:=dX/2;
@@ -755,102 +887,178 @@ begin
inherited destroy;
end;
+procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
+var
+ i,j: longint;
+begin
+ for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
+ for j:=0 to felders[aktuelleFelder].matAnz-1 do
+ if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
+ (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
+ max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10))
+ > pDNMax then begin
+ if pDNMaxDynamisch then begin
+ pDNMax:=
+ 2*
+ abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]*
+ (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
+ Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10));
+ pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')');
+ end
+ else
+ if pDNMax>0 then begin
+ pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+
+ 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true);
+ pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true);
+ pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true);
+ pDNMax:=-1;
+ end;
+ end;
+end;
+
+function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert?
+begin
+ mDT:=min(mDT,felders[aF].maxDT);
+ // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
+ // unphysikalischen Effekten (z.B. negativen Dichten) führt
+
+ result:=mDT<dT*2;
+ if result then begin // beabsichtigter Zeitschritt ist zu groß,
+ dT:=mDT/4; // dann machen wir ihn doch kleiner!
+ {$IFDEF Zeitschrittueberwachung}
+ prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT));
+ if dT<1E-30 then begin
+ prot.schreibe('pruefeMaxDT: Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true);
+ prot.destroyall;
+ halt(1);
+ end;
+ {$ENDIF}
+ exit;
+ end;
+end;
+
procedure tGitter.iteriereSchritt(var dT: extended);
var i: longint;
{$IFDEF Zeitschrittueberwachung}
- pro: tProtokollant;
j: longint;
{$ENDIF}
+ {$IFDEF Dichteueberwachung}
+ pro: tProtokollant;
+ {$ENDIF}
iDT,mDT: extended;
begin
- kvs.add('t',t+dT);
+ kvs.add('t',t);
diffusionsTermAnpassen(prot);
+ mDT:=infinity;
+
repeat
iDT:=1/dT;
felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t))
- mDT:=felders[aktuelleFelder].maxDT;
- // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
- // unphysikalischen Effekten (z.B. negativen Dichten) führt
-
- if mDT<dT*2 then begin // beabsichtigter Zeitschritt ist zu groß
- dT:=mDT/4; // machen wir ihn kleiner!
- {$IFDEF Zeitschrittueberwachung}
- pro:=tProtokollant.create(prot,'iteriereSchritt');
- pro.schreibe('dT -> '+floattostr(dT));
- if dT<1E-30 then begin
- pro.schreibe('Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true);
- pro.destroyall;
- halt(1);
- end;
- pro.free;
- {$ENDIF}
+ if pruefeMaxDT(aktuelleFelder,mDT,dT) then
continue;
- end;
+
+ case zeitverfahren of
+ zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
+ felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT);
+ zfRungeKuttaDreiAchtel:
+ begin
+ felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/3); // ya = y(t) + y' dt/3
+ felders[2].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/3,ya)
+
+ if pruefeMaxDT(2,mDT,dT) then
+ continue;
+
+ felders[3].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],-dT/3,dT); // yb = y(t) - y' dt/3 + ya' dt
+ felders[3].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // yb' = y'(t+2dt/3,yb)
+
+ if pruefeMaxDT(3,mDT,dT) then
+ continue;
+
+ felders[4].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[3],dT,-dT,dT); // yc = y(t) + a' dt - ya' dt + yb' dt
+ felders[4].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc)
+
+ if pruefeMaxDT(4,mDT,dT) then
+ continue;
+
+ felders[1-aktuelleFelder].liKo(
+ felders[aktuelleFelder],
+ felders[aktuelleFelder],
+ felders[2],
+ felders[3],
+ felders[4],
+ dT/8,
+ dT*3/8,
+ dT*3/8,
+ dT/8
+ ); // y(t+dt) = y(t) + (y' + 3(ya' + yb') + yc') dt/8
+ end;
+ zfRungeKuttaVier:
+ begin
+ felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2
+ felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/2,ya)
+
+ if pruefeMaxDT(2,mDT,dT) then
+ continue;
+
+ felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2
+ felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yb' = y'(t+dt/2,yb)
+
+ if pruefeMaxDT(3,mDT,dT) then
+ continue;
+
+ felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt
+ felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc)
+
+ if pruefeMaxDT(4,mDT,dT) then
+ continue;
+
+ felders[1-aktuelleFelder].liKo(
+ felders[aktuelleFelder],
+ felders[aktuelleFelder],
+ felders[2],
+ felders[3],
+ felders[4],
+ dT/6,
+ dT/3,
+ dT/3,
+ dT/6
+ ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
+ end;
+ end{of case};
break;
until false;
- if dT<mDT/1000 then
- dT:=min(1,mDT)/1000;
-
- case zeitverfahren of
- zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
- felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT);
- zfRungeKuttaVier:
- begin
- felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2
- felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/2,ya)
-
- felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2
- felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yb' = y'(t+dt/2,yb)
-
- felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt
- felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc)
-
- felders[1-aktuelleFelder].liKo(
- felders[aktuelleFelder],
- felders[aktuelleFelder],
- felders[2],
- felders[3],
- felders[4],
- dT/6,
- dT/3,
- dT/3,
- dT/6
- ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
- end;
- end{of case};
-
aktuelleFelder:=1-aktuelleFelder;
{$IFDEF Dichteueberwachung}
- with felders[aktuelleFelder] do
- for i:=0 to length(inhalt)-1 do
- for j:=0 to matAnz-1 do
- if inhalt[i].matWerte[j,mfN,false]<0 then begin
- pro:=tProtokollant.create(prot,'iteriereSchritt');
- pro.schreibe(
- 'n<=0 bei '+
- floattostr(t)+' '+
- inttostr(i)+' '+
- floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
- inhalt[i-1].matWerte[j,mfIGamma,false])+' '+
- floattostr(inhalt[i-1].matWerte[j,mfN,false])+' '+
- floattostr(inhalt[i-1].matWerte[j,mfN,true])+' '+
- floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
- inhalt[i].matWerte[j,mfIGamma,false])+' '+
- floattostr(inhalt[i].matWerte[j,mfN,false])+' '+
- floattostr(inhalt[i].matWerte[j,mfN,true])+' '+
- floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
- inhalt[i+1].matWerte[j,mfIGamma,false])+' '+
- floattostr(inhalt[i+1].matWerte[j,mfN,false])+' '+
- floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
- pro.destroyall;
- halt(1);
- end;
+ for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
+ for j:=0 to felders[aktuelleFelder].matAnz-1 do
+ if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then
+ with felders[1-aktuelleFelder] do begin
+ pro:=tProtokollant.create(prot,'iteriereSchritt');
+ pro.schreibe('n<0 bei:',true);
+ pro.schreibe(' t = ' +floattostr(t)+',',true);
+ pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true);
+ pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
+ inhalt[i-1].matWerte[j,mfIGamma,false])+',',true);
+ pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true);
+ pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true);
+ pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
+ inhalt[i].matWerte[j,mfIGamma,false])+',',true);
+ pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true);
+ pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true);
+ pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
+ inhalt[i+1].matWerte[j,mfIGamma,false])+',',true);
+ pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true);
+ pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
+ pro.destroyall;
+ halt(1);
+ end;
{$ENDIF}
t:=t+dT;
+ dT:=max(dT,min(1,mDT)/100);
end;
function tGitter.dumpErhaltungsgroessen: boolean;
@@ -885,35 +1093,6 @@ begin
pro.free;
end;
-procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
-var
- i,j: longint;
-begin
- for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
- for j:=0 to felders[aktuelleFelder].matAnz-1 do
- if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
- (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
- max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10))
- > pDNMax then begin
- if pDNMaxDynamisch then begin
- pDNMax:=
- 2*
- abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]*
- (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
- Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10));
- pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')');
- end
- else
- if pDNMax>0 then begin
- pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+
- 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true);
- pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true);
- pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true);
- pDNMax:=-1;
- end;
- end;
-end;
-
procedure tGitter.berechne(was: char; teilchen: longint);
var
i: longint;
@@ -957,6 +1136,7 @@ var
teilchen: array of tTeilchenSpezies;
lichter: tMyStringlist;
pro: tProtokollant;
+ na: pSigActionRec;
begin
inherited create;
prot:=tProtokollant.create(protokollant,name);
@@ -986,6 +1166,7 @@ begin
breite:=10.0;
pDNMax:=-1;
fortschrittsAnzeige:=false;
+ gotSighup:=false;
// Standardeinstellungen Bereich 'ausgaben'
aPrefix:=extractfilepath(paramstr(0));
@@ -1016,6 +1197,10 @@ begin
halt(1);
end;
if s='allgemeinEnde' then break;
+ if s='runge-Kutta-3/8' then begin
+ Zeitverfahren:=zfRungeKuttaDreiAchtel;
+ continue;
+ end;
if s='runge-Kutta-4' then begin
Zeitverfahren:=zfRungeKuttaVier;
continue;
@@ -1140,6 +1325,11 @@ begin
kvs.add('nmax'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMax);
continue;
end;
+ if startetMit('minimaldichte ',s) then begin
+ teilchen[length(teilchen)-1].nMin:=exprToFloat(false,s,kvs,nil);
+ kvs.add('nmin'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMin);
+ continue;
+ end;
if startetMit('verteilung ',s) then begin
teilchen[length(teilchen)-1].liesDichteFunktion(s,ifile,kvs,teilchen,prot);
continue;
@@ -1178,6 +1368,8 @@ begin
case zeitverfahren of
zfEulerVorwaerts:
pro.schreibe('Iteration mittels Euler-Vorwärts');
+ zfRungeKuttaDreiAchtel:
+ pro.schreibe('Iteration mittels Runge-Kutta-3/8');
zfRungeKuttaVier:
pro.schreibe('Iteration mittels Runge-Kutta-4');
else
@@ -1208,6 +1400,23 @@ begin
end;
pro.free;
+ if length(sighupSimulationen)=0 then begin
+ new(na);
+ na^.sa_Handler := sigActionHandler(@signalCapture);
+ fillchar(na^.sa_Mask,sizeof(na^.sa_mask),#0);
+ na^.sa_Flags := 0;
+ {$ifdef Linux} // Linux specific
+ na^.sa_Restorer := Nil;
+ {$endif}
+ if (fPSigaction(SIGHUP, na, nil) <> 0) then begin
+ writeln(stderr, 'Error: ', fpgeterrno, '.');
+ halt(1);
+ end;
+ dispose(na);
+ end;
+ setlength(sighupSimulationen,length(sighupSimulationen)+1);
+ sighupSimulationen[length(sighupSimulationen)-1]:=self;
+
gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter');
for i:=0 to length(teilchen)-1 do
teilchen[i].free;
@@ -1216,9 +1425,17 @@ begin
end;
destructor tSimulation.destroy;
+var
+ i,j: longint;
begin
gitter.free;
prot.free;
+ for i:=length(sighupSimulationen)-1 downto 0 do
+ if sighupSimulationen[i]=self then begin
+ for j:=i+1 to length(sighupSimulationen)-1 do
+ sighupSimulationen[j-1]:=sighupSimulationen[j];
+ setlength(sighupSimulationen,length(sighupSimulationen)-1);
+ end;
inherited destroy;
end;
@@ -1240,22 +1457,37 @@ begin
end;
zeitDatei:=zeitDatei+now;
- if fortschrittsAnzeige then begin
- if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) then begin
- prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
- prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
- prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
- prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
- if errorCode<2 then
- if not gitter.dumpErhaltungsgroessen then begin
- errorCode:=2;
- exit;
- end;
- end;
+ if gotSighup or
+ (fortschrittsAnzeige and (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde))) then begin
+ gotSighup:=false;
+ prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
+ prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
+ prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
+ prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
+ if errorCode<2 then
+ if not gitter.dumpErhaltungsgroessen then begin
+ errorCode:=2;
+ exit;
+ end;
end;
result:=gitter.t<tEnde;
end;
+// allgemeine Funktionen *******************************************************
+
+procedure SignalCapture(signal : longint); cdecl;
+var
+ i: longint;
+begin
+ if signal=SIGHUP then
+ for i:=0 to length(sighupSimulationen)-1 do
+ sighupSimulationen[i].gotSighup:=true;
+end;
+
+initialization
+ setlength(sighupSimulationen,0);
+finalization
+ setlength(sighupSimulationen,0);
end.
generated by cgit v1.2.3-54-g00ecf (git 2.44.0) at 2024-04-29 02:48:52 +0000