summaryrefslogtreecommitdiff
path: root/Physikunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'Physikunit.pas')
-rw-r--r--Physikunit.pas208
1 files changed, 115 insertions, 93 deletions
diff --git a/Physikunit.pas b/Physikunit.pas
index 27cca80..ef77ee4 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -97,6 +97,7 @@ type
public
werte: array of array[boolean] of extended; // Besetzungsdichten und deren Ableitungen für verschiedene Teilchenspezies
iMGamma: array of extended; // 1/m/gamma (= v/p)
+ grad: array of array[0..1] of extended; // Ableitung der Besetzungsdichte nach x und px
constructor create(anzTeilchen: longint); overload; // NullPunkt erstellen
constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); overload;
destructor destroy; override;
@@ -107,10 +108,8 @@ type
procedure akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended); inline;
function gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; inline; // ich habe jemanden erreicht
procedure setzeNull; inline;
- procedure berechneAbleitungen(pX,iDX,iDPX: extended); inline;
- function gradient(rtg,tlc: longint; iD,fak: extended): extended; inline;
+ procedure berechneAbleitungen(pX: extended); inline;
property istNullPunkt: boolean read istNullPkt;
- function gibNullPunkt: tImpulsPunkt;
end;
{ tRaumPunkt }
@@ -134,7 +133,7 @@ type
procedure berechneEMFelder(dX,iDX: extended); overload; inline;
procedure berechneEMQuellen;
procedure berechneEMAbleitungen(iDX: extended); inline;
- procedure berechnePhasenraumAbleitungen(iDX: extended); inline;
+ procedure berechnePhasenraumAbleitungen; inline;
{$DEFINE LiKotRaumPunktHeader}
{$INCLUDE linearkombinationen.inc}
{$UNDEF LiKotRaumPunktHeader}
@@ -163,6 +162,7 @@ type
constructor create(groesse: longint; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aP: longint; dP: extended);
destructor destroy; override;
procedure berechneAbleitungen(dX,iDX: extended); inline;
+ procedure berechneGradienten(iDX,iDP: extended);
{$DEFINE LiKotFelderHeader}
{$INCLUDE linearkombinationen.inc}
{$UNDEF LiKotFelderHeader}
@@ -399,7 +399,7 @@ begin
if sT>=nNum+sDT/2 then
schreibeKopf;
- cnt:=floor((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)/sDX*gitter.dX+1);
+ cnt:=floor(((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)*gitter.dX+Min(gitter.dX,sDX)/2)/sDX+1);
cX:=gitter.xl;
sX:=cX+(cnt-1)*sDX;
@@ -587,8 +587,8 @@ end;
constructor tImpulsPunkt.create(anzTeilchen: longint);
var
- i: longint;
- b: boolean;
+ i,j: longint;
+ b: boolean;
begin
inherited create;
istNullPkt:=true; // als NullPunkt erstellen
@@ -598,6 +598,11 @@ begin
for i:=0 to length(werte)-1 do
for b:=false to true do
werte[i,b]:=0;
+ fillchar(grad,sizeof(grad),#0);
+ setlength(grad,anzTeilchen);
+ for i:=0 to length(grad)-1 do
+ for j:=0 to 1 do
+ grad[i,j]:=0;
fillchar(iMGamma,sizeof(iMGamma),#0);
setlength(iMGamma,length(werte));
for i:=0 to length(iMGamma)-1 do
@@ -610,9 +615,9 @@ end;
constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint);
var
- i: longint;
- b: boolean;
- np: tImpulsPunkt;
+ i,j: longint;
+ b: boolean;
+ np: tImpulsPunkt;
begin
inherited create;
istNullPkt:=false;
@@ -622,6 +627,11 @@ begin
for i:=0 to length(werte)-1 do
for b:=false to true do
werte[i,b]:=0;
+ fillchar(grad,sizeof(grad),#0);
+ setlength(grad,anzTeilchen);
+ for i:=0 to length(grad)-1 do
+ for j:=0 to 1 do
+ grad[i,j]:=0;
fillchar(iMGamma,sizeof(iMGamma),#0);
setlength(iMGamma,length(werte));
for i:=0 to length(iMGamma)-1 do
@@ -630,22 +640,19 @@ begin
nachbarn[0,true]:=nil;
nachbarn[1,false]:=linksImpuls;
nachbarn[1,true]:=nil;
- for i:=0 to 1 do
- if assigned(nachbarn[i,false]) and
- not nachbarn[i,false].istNullPunkt then
- nachbarn[i,false].nachbarn[i,true]:=self;
np:=nil;
for i:=0 to 1 do
- if assigned(nachbarn[i,false]) and
- not assigned(np) then
- np:=nachbarn[i,false].gibNullPunkt;
+ if assigned(nachbarn[i,false]) then begin
+ if nachbarn[i,false].nachbarn[i,true].istNullPunkt then
+ np:=nachbarn[i,false].nachbarn[i,true];
+ nachbarn[i,false].nachbarn[i,true]:=self;
+ end;
if not assigned(np) then
np:=tImpulsPunkt.create(anzTeilchen);
for i:=0 to 1 do
for b:=false to true do
- if not assigned(nachbarn[i,b]) then begin
+ if not assigned(nachbarn[i,b]) then
nachbarn[i,b]:=np;
- end;
end;
destructor tImpulsPunkt.destroy;
@@ -694,9 +701,9 @@ begin
raumpunkt.felder.gitter.abbrechen;
end; *)
end;
- if werte[i,false]>1E10 then begin
+ if werte[i,false]>1E100 then begin
pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'impulsPunkt.nichtnegativieren');
- pro.schreibe('An einer Stelle im Phasenraum sind bereits mehr als 10^100 Teilchen, ich breche ab!',true);
+ pro.schreibe('An einer Stelle im Phasenraum sind bereits mehr als 10^100 Teilchen, ich breche ab (t = '+floattostr(raumpunkt.felder.gitter.t)+' T)!',true);
pro.free;
raumpunkt.felder.gitter.abbrechen;
end;
@@ -755,7 +762,7 @@ begin
end;
end;
-procedure tImpulsPunkt.berechneAbleitungen(pX,iDX,iDPX: extended);
+procedure tImpulsPunkt.berechneAbleitungen(pX: extended);
var
i: longint;
tx,tp: extended;
@@ -770,73 +777,9 @@ begin
( raumpunkt.emFelder[efEX,false]
+ raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]);
werte[i,true]:=
- gradient(0,i,iDX,tX) +
- gradient(1,i,iDPX,tP);
-(* tX * (
- (nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false])
- + sign(tX)*raumpunkt.felder.teilchen[i].diffusion[0] * (
- nachbarn[0,true].werte[i,false] +
- nachbarn[0,false].werte[i,false] -
- 2*werte[i,false]
- )
- ) * iDX/2
- + tP * (
- (nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false])
- + sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * (
- nachbarn[1,true].werte[i,false] +
- nachbarn[1,false].werte[i,false] -
- 2*werte[i,false]
- )
- )* iDPX/2; *)
+ grad[i,0]*tX +
+ grad[i,1]*tP;
end;
-
- // Der jeweils zweite Summand entspricht (in etwa) einer thermischen
- // Verschmierung und soll die numerische Stabilität bis zu
- // d(ln n)/d[xp] * d[XP] = "diffusion" gewährleisten.
-end;
-
-function tImpulsPunkt.gradient(rtg,tlc: longint; iD,fak: extended): extended;
-begin
-
- {
- // Variante mit Diffusionsterm
- result:=
- fak * ((nachbarn[rtg,true].werte[tlc,false] - nachbarn[rtg,false].werte[tlc,false])
- + sign(fak)*raumpunkt.felder.teilchen[tlc].diffusion[rtg] * (
- nachbarn[rtg,true].werte[tlc,false] +
- nachbarn[rtg,false].werte[tlc,false] -
- 2*werte[tlc,false]
- )
- ) * iD/2;
- }
-
- result:=
- fak * ( - 1/60 * nachbarn[rtg,false].nachbarn[rtg,false].nachbarn[rtg,false].werte[tlc,false]
- + 3/20 * nachbarn[rtg,false].nachbarn[rtg,false].werte[tlc,false]
- - 3/4 * nachbarn[rtg,false].werte[tlc,false]
- + 3/4 * nachbarn[rtg,true].werte[tlc,false]
- - 3/20 * nachbarn[rtg,true].nachbarn[rtg,true].werte[tlc,false]
- + 1/60 * nachbarn[rtg,true].nachbarn[rtg,true].nachbarn[rtg,true].werte[tlc,false]
- ) * iD;
-
-
-// probier mal den Gradienten anders zu berechnen und zwar so:
-// https://en.wikipedia.org/wiki/Finite_difference_coefficient
-end;
-
-function tImpulsPunkt.gibNullPunkt: tImpulsPunkt;
-var
- i: longint;
-begin
- if istNullPunkt then begin
- result:=self;
- exit;
- end;
- result:=nil;
- for i:=0 to 1 do
- if assigned(nachbarn[i,false]) and
- not assigned(result) then
- result:=nachbarn[i,false].gibNullPunkt;
end;
// tRaumPunkt ******************************************************************
@@ -867,7 +810,7 @@ begin
rN:=nil;
if assigned(lN) then
lN.rN:=self;
- aP:=_aP+byte(not odd(_aP));
+ aP:=_aP;
dP:=_dP;
fillchar(phasenraum,sizeof(phasenraum),#0);
setlength(phasenraum,aP);
@@ -956,13 +899,13 @@ begin
emFelder[efAY,true]:=emFelder[efDAYDT,false];
end;
-procedure tRaumPunkt.berechnePhasenraumAbleitungen(iDX: extended);
+procedure tRaumPunkt.berechnePhasenraumAbleitungen;
var
i: longint;
begin
// df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f
for i:=0 to aP-1 do
- phasenraum[i].berechneAbleitungen((i-(aP div 2))*dP,iDX,1/dP);
+ phasenraum[i].berechneAbleitungen((i-(aP div 2))*dP);
end;
function tRaumPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen
@@ -1128,7 +1071,7 @@ begin
for i:=0 to length(teilchen)-1 do
teilchen[i]:=tTeilchenSpezies.create(_teilchen[i]);
fillchar(inhalt,sizeof(inhalt),#0);
- setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen
+ setlength(inhalt,groesse);
inhalt[0]:=tRaumPunkt.create(nil,self,length(teilchen),aP,dP);
for i:=1 to length(inhalt)-1 do
inhalt[i]:=tRaumPunkt.create(inhalt[i-1],self,length(teilchen),aP,dP);
@@ -1190,6 +1133,74 @@ begin
exprToFloat(false,lichters[rechts][i],gitter.kvs,nil);
end;
+procedure tFelder.berechneGradienten(iDX,iDP: extended);
+var
+ res,ims: array of extended;
+ i,j,tlc,len: longint;
+ tmp: extended;
+begin
+
+ // Ableitungen nach PX
+
+ write('.');
+
+ len:=length(inhalt[0].phasenraum);
+ write(len);
+ setlength(res,len);
+ setlength(ims,len);
+ for i:=0 to length(inhalt)-1 do
+ for tlc:=0 to length(teilchen)-1 do begin
+ if len<>length(inhalt[i].phasenraum) then
+ raise exception.create('Unterschiedliche Diskretisierung im Phasenraum an unterschiedlichen Orten!');
+ for j:=0 to length(inhalt[i].phasenraum)-1 do begin
+ res[j]:=inhalt[i].phasenraum[j].werte[tlc,false];
+ ims[j]:=0;
+ end;
+ fft(res,ims,false);
+ for j:=0 to (len div 2)-1 do begin
+ tmp:=res[j];
+ res[j]:=ims[j]*j;
+ ims[j]:=-tmp*j;
+
+ tmp:=res[len-1-j];
+ res[len-1-j]:=ims[len-1-j]*(-j-1);
+ ims[len-1-j]:=-tmp*(-j-1);
+ end;
+ fft(res,ims,true);
+ for j:=0 to length(inhalt[i].phasenraum)-1 do
+ inhalt[i].phasenraum[j].grad[tlc,0]:=res[j]*iDP*0;
+ end;
+
+ // Ableitungen nach X
+
+ write('.');
+ len:=length(inhalt);
+ write(len);
+ setlength(res,2*len);
+ setlength(ims,2*len);
+ for i:=0 to length(inhalt[0].phasenraum)-1 do
+ for tlc:=0 to length(teilchen)-1 do begin
+ for j:=0 to length(inhalt)-1 do begin
+ res[j]:=inhalt[j].phasenraum[i].werte[tlc,false];
+ ims[j]:=0;
+ end;
+ fft(res,ims,false);
+ for j:=0 to (len div 2)-1 do begin
+ tmp:=res[j];
+ res[j]:=ims[j]*j;
+ ims[j]:=-tmp*j;
+
+ tmp:=res[len-1-j];
+ res[len-1-j]:=ims[len-1-j]*(-j-1);
+ ims[len-1-j]:=-tmp*(-j-1);
+ end;
+ fft(res,ims,true);
+ for j:=0 to length(inhalt)-1 do
+ inhalt[j].phasenraum[i].grad[tlc,1]:=res[j]*iDX*0;
+ end;
+ write('>');
+end;
+
procedure tFelder.berechneAbleitungen(dX,iDX: extended);
var
i: longint;
@@ -1214,8 +1225,10 @@ begin
for i:=1 to length(inhalt)-2 do // sonstige Ableitungen des A-Feldes berechnen
inhalt[i].berechneEMAbleitungen(iDX);
+ berechneGradienten(iDX,1/inhalt[0].dP); // Gradienten der Phasenraumbesetzungsdichte berechnen
+
for i:=1 to length(inhalt)-2 do // Ableitungen der Phasenraumbesetzungsdichten berechnen (Rand wird ausgelassen!)
- inhalt[i].berechnePhasenraumAbleitungen(iDX);
+ inhalt[i].berechnePhasenraumAbleitungen;
end;
procedure tFelder.setzeNull;
@@ -1251,6 +1264,15 @@ var
i: longint;
begin
inherited create;
+
+ i:=aP;
+ aP:=round(power(2,ceil(ln(aP-0.5)/ln(2))));
+ dP:=dP*(i-1)/(aP-1);
+
+ i:=size;
+ size:=round(power(2,ceil(ln(size+4-0.5)/ln(2)))); // zwei Felder links und rechts extra für Randbedingungen
+ deltaX:=deltaX*(i-1)/(size-4-1);
+
abbruch:=false;
besitzer:=derBesitzer;
zeitverfahren:=zv;
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-11-21 23:18:15 +0000