summaryrefslogtreecommitdiff
path: root/Physikunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'Physikunit.pas')
-rw-r--r--Physikunit.pas202
1 files changed, 78 insertions, 124 deletions
diff --git a/Physikunit.pas b/Physikunit.pas
index 0a2ca5e..142ee4b 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -52,7 +52,7 @@ type
destructor destroy; override;
function dump: string;
procedure schreibeKopf;
- procedure speichereWerte(gitter: tGitter; sDX: extended);
+ procedure speichereWerte(gitter: tGitter; sT, sDX: extended);
end;
{ tTeilchenSpezies }
@@ -121,14 +121,14 @@ type
matAnz: longint;
spezLadungen: tExtendedArray;
lichters: array[boolean] of tMyStringlist;
- procedure setzeRaender(dT,iDX,iDT: extended);
+ procedure setzeRaender(iDX: extended);
public
inhalt: array of tWertePunkt;
par: tGitter;
constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
destructor destroy; override;
- procedure berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
+ procedure berechneAbleitungen(dX,iDX,pDNMax: extended);
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;
@@ -163,7 +163,8 @@ type
dX,iDX,pDNMax,dTMaximum,xl,t: extended;
kvs: tKnownValues;
procedure diffusionsTermAnpassen(pro: tProtokollant);
- function pruefeMaxDT(aF: longint; var mDT,dT: extended; dTMin: extended): boolean; // wurde verkleinert?
+ function pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert?
+ procedure abbrechen;
public
aktuelleFelder: longint;
@@ -334,7 +335,7 @@ begin
closefile(datei);
end;
-procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sDX: extended);
+procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: extended);
var i,cnt: longint;
sX,cX: extended;
begin
@@ -348,7 +349,7 @@ begin
if (teilchen>=0) and (matInhalt=mfP) then
gitter.berechne('P',teilchen);
- if gitter.t>=nNum+sDT/2 then
+ if sT>=nNum+sDT/2 then
schreibeKopf;
cnt:=floor((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)/sDX*gitter.dX+1);
@@ -749,12 +750,11 @@ begin
inherited destroy;
end;
-procedure tFelder.setzeRaender(dT,iDX,iDT: extended);
+procedure tFelder.setzeRaender(iDX: extended);
var
emF: tEMFeldInhalt;
rechts: boolean;
i: longint;
- nVal: extended;
begin
for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld
inhalt[0].emWerte[emF,true]:=
@@ -766,17 +766,13 @@ begin
end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler bzw. nicht beachtete numerische Dispersion)
for rechts:=false to true do
- for i:=0 to lichters[rechts].count-1 do begin
- par.kvs.add('t',par.t+dT);
- nVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil);
- par.kvs.add('t',par.t);
+ for i:=0 to lichters[rechts].count-1 do
inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true]:=
inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true] +
- (nVal - exprToFloat(false,lichters[rechts][i],par.kvs,nil))*iDT;
- end;
+ exprToFloat(false,lichters[rechts][i],par.kvs,nil);
end;
-procedure tFelder.berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
+procedure tFelder.berechneAbleitungen(dX,iDX,pDNMax: extended);
var
i: longint;
begin
@@ -795,7 +791,7 @@ begin
inhalt[1].berechneNAbleitung(iDX,pDNMax,false);
inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true);
- setzeRaender(dT,iDX,iDT);
+ setzeRaender(iDX);
for i:=1 to length(inhalt)-2 do
inhalt[i].berechneAbleitungen(dX,iDX);
@@ -863,49 +859,47 @@ end;
procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
var
- i,j: longint;
+ i,j: longint;
+ curMax: extended;
begin
+ curMax:=0;
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]*
+ curMax:=
+ max(curMax,
+ 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;
+ max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-100)));
+ if curMax > pDNMax then begin
+ if pDNMaxDynamisch then begin
+ pDNMax:=2*curMax;
+ 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',true);
+ pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern',true);
+ pDNMax:=0; // Term komplett abschalten
+ end;
+ end;
end;
-function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended; dTMin: extended): boolean; // wurde verkleinert?
+function tGitter.pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert?
begin
- mDT:=min(mDT,felders[aF].maxDT);
- // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
+ dTMax:=min(dTMax,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;
+ result:=dTMax<dT*2;
if result then begin // beabsichtigter Zeitschritt ist zu groß,
- dT:=mDT/4; // dann machen wir ihn doch kleiner!
+ dT:=dTMax/4; // dann machen wir ihn doch kleiner!
{$IFDEF Zeitschrittueberwachung}
- prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT));
+ prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')');
if dT<dTMin then begin
prot.schreibe('pruefeMaxDT: Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < '+floattostr(dTMin)+') - irgendwas scheint grundsätzlich kaputt zu sein!',true);
- abbruch:=true;
-// prot.destroyall;
-// halt(1);
+ abbrechen;
end;
prot.spuelen;
{$ENDIF}
@@ -913,6 +907,26 @@ begin
end;
end;
+procedure tGitter.abbrechen;
+var
+ i,j,k: longint;
+ maF: tMaterieFeldInhalt;
+ emF: tEMFeldInhalt;
+ abl: boolean;
+begin
+ abbruch:=true;
+ for i:=0 to length(felders)-1 do
+ for j:=0 to length(felders[i].inhalt)-1 do begin
+ for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do
+ for abl:=false to true do
+ felders[i].inhalt[j].emWerte[emF,false]:=0;
+ for k:=0 to felders[i].matAnz-1 do
+ for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do
+ for abl:=false to true do
+ felders[i].inhalt[j].matWerte[k,maF,false]:=0;
+ end;
+end;
+
procedure tGitter.iteriereSchritt(var dT: extended; dTMin: extended);
var i: longint;
{$IFDEF Zeitschrittueberwachung}
@@ -921,7 +935,7 @@ var i: longint;
{$IFDEF Dichteueberwachung}
pro: tProtokollant;
{$ENDIF}
- iDT,mDT: extended;
+ dTMax: extended;
begin
if abbruch then begin
t:=t+dT;
@@ -929,97 +943,37 @@ begin
end;
kvs.add('t',t);
diffusionsTermAnpassen(prot);
-
- mDT:=infinity;
+ dTMax:=dX*100;
repeat
- iDT:=1/dT;
- felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t))
+ felders[aktuelleFelder].berechneAbleitungen(dX,iDX,pDNMax); // y' = y'(t,y(t))
- if pruefeMaxDT(aktuelleFelder,mDT,dT,dTMin) then
+ if pruefeMaxDT(aktuelleFelder,dTMax,dT,dTMin) then
continue;
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,dTMin) 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,dTMin) 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,dTMin) 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,dTMin) 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,dTMin) 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,dTMin) 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;
- zfRungeKuttaZehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk108.txt
+ zfRungeKuttaDreiAchtel: begin
+ {$INCLUDE rk3_8.inc}
+ end;
+ zfRungeKuttaVier: begin
+ {$INCLUDE rk4.inc}
+ end;
+ zfRungeKuttaZehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk108.txt
{$INCLUDE rk108.inc}
end;
- zfRungeKuttaZwoelf: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1210.txt
+ zfRungeKuttaZwoelf: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1210.txt
{$INCLUDE rk1210.inc}
end;
- zfRungeKuttaVierzehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1412.txt
+ zfRungeKuttaVierzehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1412.txt
{$INCLUDE rk1412.inc}
end;
end{of case};
break;
until abbruch;
-
+ prot.schreibe('dT = '+floattostr(dT));
aktuelleFelder:=1-aktuelleFelder;
{$IFDEF Dichteueberwachung}
for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
@@ -1047,7 +1001,7 @@ begin
end;
{$ENDIF}
t:=t+dT;
- dT:=max(dT,min(1,mDT)/100);
+ dT:=max(dT,dTMax/100);
end;
function tGitter.dumpErhaltungsgroessen: boolean;
@@ -1462,7 +1416,7 @@ var
begin
result:=false;
- if gitter.abbruch then
+ if gitter.abbruch or (dT>sDT) then
dT:=sDT;
zeitPhysik:=zeitPhysik-now;
@@ -1471,9 +1425,9 @@ begin
zeitPhysik:=zeitPhysik+now;
zeitDatei:=zeitDatei-now;
while (gitter.t>=sT) and (sT<tEnde) do begin
- sT:=sT+sDT;
for i:=0 to length(ausgabeDateien)-1 do
- ausgabeDateien[i].speichereWerte(gitter,sDX);
+ ausgabeDateien[i].speichereWerte(gitter,sT,sDX);
+ sT:=sT+sDT;
end;
zeitDatei:=zeitDatei+now;
generated by cgit v1.2.3-54-g00ecf (git 2.44.0) at 2024-04-29 14:30:57 +0000