summaryrefslogtreecommitdiff
path: root/Physikunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'Physikunit.pas')
-rw-r--r--Physikunit.pas608
1 files changed, 608 insertions, 0 deletions
diff --git a/Physikunit.pas b/Physikunit.pas
new file mode 100644
index 0000000..0aee452
--- /dev/null
+++ b/Physikunit.pas
@@ -0,0 +1,608 @@
+unit Physikunit;
+
+{$mode objfpc}{$H+}
+
+{ $DEFINE Zeitschrittueberwachung}
+{ $DEFINE Dichteueberwachung}
+
+interface
+
+uses
+ Classes, SysUtils, Math;
+
+type
+ TZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier);
+ TVerteilungsfunktion = function(x: extended): extended;
+ TFeldInhalt = (
+ fiA,fiAX,fiAY,fiAZ,
+ fidAXDT,fidAYDT,fidAZDT,
+ fiN,fidPhiDX,fidPsiDX,
+ fiP,fiPX,fiPY,fiPZ,
+ fiGamma,fiiGamma);
+
+ TAusgabeDatei = record
+ Ableitung: Boolean;
+ Inhalt: TFeldInhalt;
+ Datei: File;
+ end;
+
+ TProtokolldatei = record
+ datei: textfile;
+ einrueckung: longint;
+ end;
+
+ { TProtokollant }
+
+ TProtokollant = class(TObject)
+ private
+ sDat: ^TProtokolldatei;
+ Bes: String;
+ Kinder: array of TProtokollant;
+ Elter: TProtokollant;
+ public
+ constructor create(Dateiname: string); overload;
+ constructor create(Elter_: TProtokollant; Besitzer: string); overload;
+ destructor destroy; override;
+ procedure schreibe(s: string);
+ procedure schreibe(s: string; tee: boolean);
+ procedure destroyall;
+ end;
+
+ { TFeld }
+
+ pTFeld = ^TFeld;
+ TFeld = class(TObject)
+ public
+ Groessen: array[Boolean,TFeldInhalt] of extended;
+ // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
+ dX,iDX,x,pDNMax: extended;
+ lN,rN: pTFeld;
+ constructor create(deltaX,pDNMa: extended);
+ procedure berechneGammaUndP;
+ procedure berechneNAbleitung;
+ procedure berechneAbleitungen(dT: extended);
+ end;
+
+ { TGitter }
+
+ TGitter = class(TObject)
+ private
+ groesse: longint;
+ Prot: TProtokollant;
+ procedure berechneAbleitungen(Felder: longint; dt: extended); overload;
+ procedure berechneAbleitungen(Felder: longint; dt: extended; out dTMax: extended); overload;
+ procedure setzeRaender(Felder: longint);
+ public
+ aktuelleFelder: longint;
+ Felders: array of array of TFeld;
+ dTMaximum,dX,iDX,xl,xr,t,pDNMax: extended;
+ pDNMaxDynamisch: boolean;
+ Zeitverfahren: TZeitverfahren;
+ Al,Ar: TVerteilungsfunktion;
+ constructor create(size: longint; deltaT,deltaX,pDNMa: extended; n0: TVerteilungsfunktion; ZV: TZeitverfahren; Protokollant: TProtokollant);
+ destructor destroy; override;
+ procedure iteriereSchritt(var dT: extended);
+ procedure macheAusgabe(AD: TAusgabedatei; sDX: extended);
+ function gibErhaltungsgroessen: string;
+ end;
+
+function Nullfunktion(x: extended): extended;
+function timetostr(t: extended): string;
+
+const
+ OptionNamen: array[TFeldInhalt] of string = (
+ 'A','AX','AY','AZ','dAXDT','dAYDT','dAZDT','N','dPhiDX','dPsiDX','P','PX','PY','PZ','Gamma','iGamma');
+ AusgabeNamen: array[TFeldInhalt] of string = (
+ 'fiA','fiAX','fiAY','fiAZ','fidAXDT','fidAYDT','fidAZDT','fiN','fidPhiDX','fidPsiDX','fiP','fiPX','fiPY','fiPZ','fiGamma','fiiGamma');
+ AusgabeHilfe = '--d?A[XYZ]? --d?dA[XYZ]dt --d?N --d?dPhiDX --d?dPsiDX --P[XYZ]? --Gamma --iGamma';
+
+implementation
+
+{ TProtokollant }
+
+constructor TProtokollant.create(Dateiname: string);
+begin
+ inherited create;
+ getmem(sDat,sizeof(TProtokolldatei));
+ sDat^.einrueckung:=10;
+ assignfile(sDat^.datei,Dateiname);
+ rewrite(sDat^.datei);
+ Bes:='';
+ setlength(Kinder,0);
+ Elter:=nil;
+end;
+
+constructor TProtokollant.create(Elter_: TProtokollant; Besitzer: string);
+begin
+ inherited create;
+ sDat:=Elter_.sDat;
+ setlength(Elter_.Kinder,length(Elter_.Kinder)+1);
+ Elter_.Kinder[length(Elter_.Kinder)-1]:=self;
+ Elter:=Elter_;
+ setlength(Kinder,0);
+ Bes:=Elter.Bes+'.'+Besitzer;
+ if pos('.',Bes)=1 then delete(Bes,1,1);
+ if length(Bes)+4 > sDat^.einrueckung then
+ sDat^.einrueckung:=length(Bes)+4;
+end;
+
+destructor TProtokollant.destroy;
+var i: longint;
+begin
+ while length(Kinder)>0 do
+ Kinder[length(Kinder)-1].free;
+ if assigned(Elter) then begin
+ flush(sDat^.datei);
+ i:=0;
+ while (i<length(Elter.Kinder)) and (Elter.Kinder[i]<>self) do
+ inc(i);
+ if i>=length(Elter.Kinder) then
+ schreibe('destroy fehlgeschlagen, ich kann mich nicht sehen.');
+ inc(i);
+ while i<length(Elter.Kinder) do begin
+ Elter.Kinder[i-1]:=Elter.Kinder[i];
+ inc(i);
+ end;
+ setlength(Elter.Kinder,length(Elter.Kinder)-1);
+ end
+ else begin
+ closefile(sDat^.datei);
+ freemem(sDat,sizeof(TProtokolldatei));
+ end;
+ inherited destroy;
+end;
+
+procedure TProtokollant.destroyall;
+begin
+ if Assigned(Elter) then Elter.destroyall
+ else destroy;
+end;
+
+procedure TProtokollant.schreibe(s: string);
+begin
+ schreibe(s,false);
+end;
+
+procedure TProtokollant.schreibe(s: string; tee: boolean);
+var i: longint;
+begin
+ for i:=length(Bes)+1 to sDat^.einrueckung do
+ s:=' '+s;
+ writeln(sDat^.datei,Bes+s);
+ if tee then
+ writeln(Bes+s);
+end;
+
+{ TFeld }
+
+constructor TFeld.create(deltaX,pDNMa: extended);
+var FI: TFeldinhalt;
+ Abl: Boolean;
+begin
+ inherited create;
+ lN:=nil;
+ rN:=nil;
+ for Abl:=False to True do
+ for FI:=low(TFeldinhalt) to high(TFeldinhalt) do
+ Groessen[Abl,FI]:=0;
+ dX:=deltaX;
+ iDX:=1/dX;
+ x:=0;
+ pDNMax:=pDNMa;
+ Groessen[false,fiN]:=1;
+ Groessen[false,fiGamma]:=1;
+ Groessen[false,fiiGamma]:=1;
+end;
+
+procedure TFeld.berechneGammaUndP; // aus aktuellem dPsi/dX und A
+var b: byte;
+ tmp: extended;
+begin
+ tmp:=0;
+ for b:=0 to 2 do begin
+ Groessen[false,TFeldinhalt(Byte(fiPX)+b)]:=
+ Groessen[false,TFeldinhalt(Byte(fiAX)+b)] + Groessen[false,fidPsiDX]*Byte(b=0);
+ tmp:=tmp+sqr(Groessen[false,TFeldinhalt(Byte(fiPX)+b)]);
+ end;
+ Groessen[false,fiGamma]:=sqrt(1+tmp);
+ Groessen[false,fiiGamma]:=1/Groessen[false,fiGamma];
+end;
+
+procedure TFeld.berechneNAbleitung; // Zeitableitung von n berechnen
+begin
+(* // dn/dt = -d(n/Gamma p_x)/dx
+ Groessen[true,fiN]:=
+ -( rN^.rN^.Groessen[false,fiN]*rN^.rN^.Groessen[false,fiiGamma]*rN^.rN^.Groessen[false,fiPX]
+ + rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*rN^.Groessen[false,fiPX]
+ - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*lN^.Groessen[false,fiPX]
+ - lN^.lN^.Groessen[false,fiN]*lN^.lN^.Groessen[false,fiiGamma]*lN^.lN^.Groessen[false,fiPX]
+ )*iDX/6;
+ // es wird über die beiden nächsten linken bzw. rechten Nachbarn gemittelt *)
+
+ // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
+ Groessen[true,fiN]:=
+ ( rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*(pDNMax - rN^.Groessen[false,fiPX])
+ - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*(-pDNMax - lN^.Groessen[false,fiPX])
+ - Groessen[false,fiN]*Groessen[false,fiiGamma]*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 TFeld.berechneAbleitungen(dT: extended); // Zeitableitungen berechnen
+var b: byte;
+begin
+ // d2Psi/dxdt = dPhi/dx - dGamma/dx
+ Groessen[true,fidPsiDX]:=
+ Groessen[false,fidPhiDX]
+ - (rN^.Groessen[false,fiGamma] - lN^.Groessen[false,fiGamma])*iDX/2;
+ // d3Phi/dx2dt = dn/dt ... hier muss integriert werden:
+ // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + <dn/dt> * deltax
+ Groessen[true,fidPhiDX]:=
+ lN^.Groessen[true,fidPhiDX]
+ + (lN^.Groessen[true,fiN] + Groessen[true,fiN])*dX/2;
+ // d2A/dt2 = - n/gamma p + Laplace(A) ...
+ for b:=0 to 2 do
+ Groessen[true,TFeldinhalt(Byte(fidAXDT)+b)]:=
+ - (Groessen[false,fiN]*Groessen[false,fiiGamma]*Groessen[false,TFeldinhalt(Byte(fiPX)+b)])
+ + (rN^.Groessen[false,TFeldinhalt(Byte(fiAX)+b)]
+ - 2*Groessen[false,TFeldinhalt(Byte(fiAX)+b)]
+ + lN^.Groessen[false,TFeldinhalt(Byte(fiAX)+b)])*sqr(iDX);
+ Groessen[true,fidAXDT]:=Groessen[true,fidAXDT] - Groessen[false,fidPhiDX]; // ... - Nabla(phi)
+ // dA/dt = dA/dt
+ for b:=0 to 2 do
+ Groessen[true,TFeldInhalt(Byte(fiAX)+b)]:=
+ Groessen[false,TFeldInhalt(Byte(fidAXDT)+b)] +
+ Groessen[true,TFeldInhalt(Byte(fidAXDT)+b)]*dT;
+end;
+
+{ TGitter }
+
+procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended; out dTMax: extended);
+var i: longint;
+begin
+ berechneAbleitungen(Felder,dT);
+ dTMax:=dTMaximum;
+ for i:=0 to groesse+3 do
+ if Felders[Felder,i].Groessen[true,fiN]<0 then begin
+ dTMax:=min(dTMax,Abs(Felders[Felder,i].Groessen[false,fiN]/Felders[Felder,i].Groessen[true,fiN]));
+ end;
+end;
+
+procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended);
+var i: longint;
+begin
+ for i:=0 to groesse+3 do begin
+(* Felders[Felder,i].Groessen[false,fiN]:=
+ max(Felders[Felder,i].Groessen[false,fiN],0); // n >= 0 *)
+ Felders[Felder,i].berechneGammaUndP;
+ end;
+ Felders[Felder,1].Groessen[false,fiPX]:= // Teilchen werden reflektiert
+ Abs(Felders[Felder,1].Groessen[false,fiPX]);
+ Felders[Felder,groesse+2].Groessen[false,fiPX]:=
+ -Abs(Felders[Felder,groesse+2].Groessen[false,fiPX]);
+ setzeRaender(Felder);
+ for i:=1 to groesse+2 do
+ Felders[Felder,i].berechneNAbleitung;
+ for i:=1 to groesse+2 do
+ Felders[Felder,i].berechneAbleitungen(dT);
+end;
+
+procedure TGitter.setzeRaender(Felder: longint);
+var FI: TFeldInhalt;
+begin
+ for FI:=fiAX to fiAZ do begin // Vakuumrandbedingungen für das A-Feld
+ Felders[Felder,0].Groessen[true,FI]:=
+ (Felders[Felder,1].Groessen[false,FI] -
+ Felders[Felder,0].Groessen[false,FI])*iDX;
+ Felders[Felder,groesse+3].Groessen[true,FI]:=
+ (Felders[Felder,groesse+2].Groessen[false,FI] -
+ Felders[Felder,groesse+3].Groessen[false,FI])*iDX;
+ end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler)
+ Felders[Felder,0].Groessen[true,fiAy]:=
+ Felders[Felder,0].Groessen[true,fiAy]
+ + (Al(t+dTMaximum)-Al(t))/dTMaximum;
+ Felders[Felder,groesse+3].Groessen[true,fiAy]:=
+ Felders[Felder,groesse+3].Groessen[true,fiAy]
+ + (Ar(t+dTMaximum)-Ar(t))/dTMaximum;
+end;
+
+constructor TGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; n0: TVerteilungsfunktion; ZV: TZeitverfahren; Protokollant: TProtokollant);
+var i,j: longint;
+begin
+ inherited create;
+ Ar:=@Nullfunktion;
+ Al:=@Nullfunktion;
+ Zeitverfahren:=ZV;
+ groesse:=size;
+ Prot:=TProtokollant.create(Protokollant,'TGitter');
+ dTMaximum:=deltaT;
+ dX:=deltaX;
+ iDX:=1/dX;
+ if pDNMa < 0 then begin
+ pDNMax:=0;
+ pDNMaxDynamisch:=true;
+ end
+ else pDNMax:=pDNMa;
+ case Zeitverfahren of
+ zfEulerVorwaerts: Setlength(Felders,2);
+ zfRungeKuttaVier: Setlength(Felders,5);
+ end{of Case};
+ xl:=dX/2;
+ for i:=0 to length(Felders)-1 do begin
+ setlength(Felders[i],size+4); // auf jeder Seite je zwei Felder Rand extra
+ for j:=0 to length(felders[i])-1 do begin
+ felders[i,j]:=TFeld.create(dX,pDNMax);
+ felders[i,j].x:=xl+j*dX;
+ felders[i,j].Groessen[false,fiN]:=n0(felders[i,j].x);
+ felders[i,j].Groessen[false,fidPsiDX]:=0.0001*Cos(j/(length(felders[i])-1)*pi*15);
+ end;
+ felders[i,0].lN:=@felders[i,0];
+ for j:=1 to length(felders[i])-1 do begin
+ felders[i,j].lN:=@felders[i,j-1];
+ felders[i,j-1].rN:=@felders[i,j];
+ end;
+ felders[i,length(felders[i])-1].rN:=@felders[i,length(felders[i])-1];
+ end;
+ xr:=xl+dX*(length(Felders[0])-1);
+ aktuelleFelder:=0;
+ t:=0;
+end;
+
+destructor TGitter.destroy;
+var i,j: longint;
+begin
+ for i:=0 to 1 do
+ for j:=0 to length(Felders[i])-1 do
+ Felders[i,j].free;
+ inherited destroy;
+end;
+
+procedure TGitter.iteriereSchritt(var dT: extended);
+var i: longint;
+ FI: TFeldInhalt;
+ dTMax: extended;
+ {$IFDEF Zeitschrittueberwachung}
+ Pro: TProtokollant;
+ {$ENDIF}
+const zeitObergrenze = 1/4;
+ zeitSollwert = 1/8;
+ zeitUntergrenze = 1/16;
+begin
+ berechneAbleitungen(aktuelleFelder,dT,dTMax); // y' = y'(t,y(t))
+ {$IFDEF Zeitschrittueberwachung}
+ Pro:=TProtokollant.create(Prot,'iteriereSchritt');
+ {$ENDIF}
+ while dT>=dTMax*zeitObergrenze do begin
+ dT:=dTMax*zeitSollwert;
+ berechneAbleitungen(aktuelleFelder,dT,dTMax); // y' = y'(t,y(t))
+ {$IFDEF Zeitschrittueberwachung}
+ Pro.schreibe('+ neues dT: '+floattostr(dT)+' (t='+floattostr(t)+')');
+ {$ENDIF}
+ end;
+ {$IFDEF Zeitschrittueberwachung}
+ if dT<1e-30 then begin
+ for i:=1 to groesse+2 do
+ if (Felders[aktuelleFelder,i].Groessen[true,fiN]<0) and
+ (Abs(Felders[aktuelleFelder,i].Groessen[false,fiN]/
+ Felders[aktuelleFelder,i].Groessen[true,fiN])<1e-15) then
+ Prot.schreibe(
+ floattostr(t)+' '+
+ inttostr(i)+' '+
+ floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fidPsiDX]*
+ Felders[aktuelleFelder,i-1].Groessen[false,fiiGamma])+' '+
+ floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i-1].Groessen[true,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i].Groessen[false,fidPsiDX]*
+ Felders[aktuelleFelder,i].Groessen[false,fiiGamma])+' '+
+ floattostr(Felders[aktuelleFelder,i].Groessen[false,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i].Groessen[true,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fidPsiDX]*
+ Felders[aktuelleFelder,i+1].Groessen[false,fiiGamma])+' '+
+ floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i+1].Groessen[true,fiN]),
+ true);
+ Pro.destroyall;
+ halt(1);
+ end;
+ {$ENDIF}
+
+ case Zeitverfahren of
+ zfEulerVorwaerts:
+ begin // y(t+dt) = y(t) + y' dt
+ for i:=0 to groesse+3 do
+ for FI:=low(TFeldinhalt) to fidPsiDX do
+ Felders[1-aktuelleFelder,i].Groessen[false,FI]:=
+ Felders[aktuelleFelder,i].Groessen[false,FI]
+ + Felders[aktuelleFelder,i].Groessen[true,FI]*dT;
+ end;
+ zfRungeKuttaVier:
+ begin
+ for i:=0 to groesse+3 do // ya = y(t) + y' dt/2
+ for FI:=low(TFeldinhalt) to fidPsiDX do
+ Felders[2,i].Groessen[false,FI]:=
+ Felders[aktuelleFelder,i].Groessen[false,FI]
+ + Felders[aktuelleFelder,i].Groessen[true,FI]*dT/2;
+ berechneAbleitungen(2,dT/2,dTMax); // ya' = y'(t+dt/2,ya)
+ if dT/2>dTMax*zeitObergrenze then begin
+ {$IFDEF Zeitschrittueberwachung}
+ Pro.schreibe('1. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true);
+ Pro.free;
+ {$ENDIF}
+ dT:=dTMax*zeitSollwert;
+ exit;
+ end;
+ for i:=0 to groesse+3 do // yb = y(t) + ya' dt/2
+ for FI:=low(TFeldinhalt) to fidPsiDX do
+ Felders[3,i].Groessen[false,FI]:=
+ Felders[aktuelleFelder,i].Groessen[false,FI]
+ + Felders[2,i].Groessen[true,FI]*dT/2;
+ berechneAbleitungen(3,dT/2,dTMax); // yb' = y'(t+dt/2,yb)
+ if dT/2>dTMax*zeitObergrenze then begin
+ {$IFDEF Zeitschrittueberwachung}
+ Pro.schreibe('2. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true);
+ Pro.free;
+ {$ENDIF}
+ dT:=dTMax*zeitSollwert;
+ exit;
+ end;
+ for i:=0 to groesse+3 do // yc = y(t) + yb' dt
+ for FI:=low(TFeldinhalt) to fidPsiDX do
+ Felders[4,i].Groessen[false,FI]:=
+ Felders[aktuelleFelder,i].Groessen[false,FI]
+ + Felders[3,i].Groessen[true,FI]*dT;
+ berechneAbleitungen(4,dT/2,dTMax); // yc' = y'(t+dt,yc)
+ if dT/2>dTMax*zeitObergrenze then begin
+ {$IFDEF Zeitschrittueberwachung}
+ Pro.schreibe('3. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true);
+ Pro.free;
+ {$ENDIF}
+ dT:=dTMax*zeitSollwert;
+ exit;
+ end;
+ for i:=0 to groesse+3 do // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
+ for FI:=low(TFeldinhalt) to fidPsiDX do
+ Felders[1-aktuelleFelder,i].Groessen[false,FI]:=
+ Felders[aktuelleFelder,i].Groessen[false,FI]
+ + (Felders[aktuelleFelder,i].Groessen[true,FI] +
+ 2*(Felders[2,i].Groessen[true,FI] +
+ Felders[3,i].Groessen[true,FI]) +
+ Felders[4,i].Groessen[true,FI])*dT/6;
+ end;
+ end{of case};
+
+ t:=t+dT;
+ aktuelleFelder:=1-aktuelleFelder;
+ if dT<dTMax*zeitUntergrenze then begin
+ dT:=dTMax*zeitSollwert;
+ {$IFDEF Zeitschrittueberwachung}
+ Pro.schreibe('- neues dT: '+floattostr(dT)+' (t='+floattostr(t)+')');
+ {$ENDIF}
+ end;
+ {$IFDEF Zeitschrittueberwachung}
+ for i:=0 to groesse+3 do
+ if Felders[aktuelleFelder,i].Groessen[false,fiN]<0 then begin
+ Pro.schreibe(
+ 'n<=0 bei '+
+ floattostr(t)+' '+
+ inttostr(i)+' '+
+ floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fidPsiDX]*
+ Felders[aktuelleFelder,i-1].Groessen[false,fiiGamma])+' '+
+ floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i-1].Groessen[true,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i].Groessen[false,fidPsiDX]*
+ Felders[aktuelleFelder,i].Groessen[false,fiiGamma])+' '+
+ floattostr(Felders[aktuelleFelder,i].Groessen[false,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i].Groessen[true,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fidPsiDX]*
+ Felders[aktuelleFelder,i+1].Groessen[false,fiiGamma])+' '+
+ floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fiN])+' '+
+ floattostr(Felders[aktuelleFelder,i+1].Groessen[true,fiN]),true);
+ Pro.destroyall;
+ halt(1);
+ end;
+ {$ENDIF}
+end;
+
+procedure TGitter.macheAusgabe(AD: TAusgabedatei; sDX: extended);
+var i: longint;
+ tmp,sX: extended;
+ b: byte;
+begin
+ if AD.Inhalt in [fiA,fiP] then
+ for i:=0 to length(Felders[aktuelleFelder])-1 do begin
+ tmp:=0;
+ for b:=1 to 3 do
+ tmp:=tmp + sqr(Felders[aktuelleFelder,i].Groessen[AD.Ableitung,TFeldinhalt(Byte(AD.Inhalt)+b)]);
+ Felders[aktuelleFelder,i].Groessen[AD.Ableitung,AD.Inhalt]:=sqrt(tmp);
+ end;
+ i:=floor((xr-xl)/sDX+1);
+ tmp:=xl+(i-1)*sDX;
+ BlockWrite(AD.Datei,xl,sizeof(extended));
+ BlockWrite(AD.Datei,tmp,sizeof(extended));
+ BlockWrite(AD.Datei,i,sizeof(longint));
+ sX:=xl-Min(dX,sDX)/2;
+ for i:=0 to length(Felders[aktuelleFelder])-1 do
+ while Felders[aktuelleFelder,i].x>=sX do begin
+ BlockWrite(AD.Datei,Felders[aktuelleFelder,i].Groessen[AD.Ableitung,AD.Inhalt],sizeof(extended));
+ sX:=sX+sDX;
+ end;
+end;
+
+function TGitter.gibErhaltungsgroessen: string;
+var i,j,k: integer;
+ n: extended;
+ Pro: TProtokollant;
+begin
+ Pro:=TProtokollant.create(Prot,'gibErhaltungsgroessen');
+ n:=0;
+ for i:=0 to groesse+3 do
+ n:=n+Felders[aktuelleFelder,i].Groessen[false,fiN];
+ for i:=1 to groesse+2 do
+ if abs(Felders[aktuelleFelder,i].Groessen[false,fiPx]*
+ (Felders[aktuelleFelder,i+1].Groessen[false,fiN]-Felders[aktuelleFelder,i-1].Groessen[false,fiN])/
+ Max(Felders[aktuelleFelder,i+1].Groessen[false,fiN]+Felders[aktuelleFelder,i-1].Groessen[false,fiN],1e-10))
+ > pDNMax then begin
+ if pDNMaxDynamisch then begin
+ pDNMax:=
+ 2*
+ abs(Felders[aktuelleFelder,i].Groessen[false,fiPx]*
+ (Felders[aktuelleFelder,i+1].Groessen[false,fiN]-Felders[aktuelleFelder,i-1].Groessen[false,fiN])/
+ Max(Felders[aktuelleFelder,i+1].Groessen[false,fiN]+Felders[aktuelleFelder,i-1].Groessen[false,fiN],1e-10));
+ for j:=0 to length(felders)-1 do
+ for k:=0 to length(felders[j])-1 do
+ felders[j,k].pDNMax:=pDNMax;
+ 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;
+ result:='n='+floattostr(n);
+ {$IFDEF Dichteueberwachung}
+ if n>1000 then begin
+ errorCode:=2;
+ Pro.schreibe(' n > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')');
+ end;
+ {$ENDIF}
+ Pro.free;
+end;
+
+function Nullfunktion(x: extended): extended;
+begin
+ result:=0*x;
+end;
+
+function timetostr(t: extended): string;
+var b: boolean;
+begin
+ result:='';
+ b:=t>=1;
+ if b then begin
+ result:=result+inttostr(floor(t))+' Tage ';
+ t:=t-floor(t);
+ end;
+ t:=t*24;
+ b:=b or (t>=1);
+ if b then begin
+ result:=result+inttostr(floor(t))+'h ';
+ t:=t-floor(t);
+ end;
+ t:=t*60;
+ b:=b or (t>=1);
+ if b then begin
+ result:=result+inttostr(floor(t))+'min ';
+ t:=t-floor(t);
+ end;
+ t:=t*60;
+ result:=result+inttostr(round(t))+'s';
+end;
+
+end.
+
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-11-24 16:20:37 +0000