From 6c4fc54c735407ef70a3de495a58600af730d07b Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Mon, 24 Aug 2015 09:40:58 +0200 Subject: Der Rand ist jetzt ein regulaerer Punkt. Gradient in Funktion gekapselt, ist jetzt O(6). --- Physikunit.pas | 108 +++++++++++++++++++++++++++++++++++++++++++++----- Plasmapropagation.lps | 103 +++++++++++++++++++++++++++++++++-------------- 2 files changed, 172 insertions(+), 39 deletions(-) diff --git a/Physikunit.pas b/Physikunit.pas index 5d8881a..27cca80 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -91,12 +91,14 @@ type tImpulsPunkt = class(tObject) // Besetzungsdichten verschiedener Teilchenspezies für ein Paar (r;p) private - raumpunkt: tRaumPunkt; - nachbarn: array[0..1,boolean] of tImpulsPunkt; // x- und px-Nachbarn (jeweil in Richtung - und +) + raumpunkt: tRaumPunkt; + nachbarn: array[0..1,boolean] of tImpulsPunkt; // x- und px-Nachbarn (jeweil in Richtung - und +) + istNullPkt: boolean; 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) - constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); + constructor create(anzTeilchen: longint); overload; // NullPunkt erstellen + constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); overload; destructor destroy; override; {$DEFINE LiKotImpulsPunktHeader} {$INCLUDE linearkombinationen.inc} @@ -106,6 +108,9 @@ type 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; + property istNullPunkt: boolean read istNullPkt; + function gibNullPunkt: tImpulsPunkt; end; { tRaumPunkt } @@ -580,12 +585,37 @@ end; // tImpulsPunkt **************************************************************** -constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); +constructor tImpulsPunkt.create(anzTeilchen: longint); var i: longint; b: boolean; begin inherited create; + istNullPkt:=true; // als NullPunkt erstellen + raumPunkt:=nil; + fillchar(werte,sizeof(werte),#0); + setlength(werte,anzTeilchen); + for i:=0 to length(werte)-1 do + for b:=false to true do + werte[i,b]:=0; + fillchar(iMGamma,sizeof(iMGamma),#0); + setlength(iMGamma,length(werte)); + for i:=0 to length(iMGamma)-1 do + iMGamma[i]:=1; + nachbarn[0,false]:=self; + nachbarn[0,true]:=self; + nachbarn[1,false]:=self; + nachbarn[1,true]:=self; +end; + +constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); +var + i: longint; + b: boolean; + np: tImpulsPunkt; +begin + inherited create; + istNullPkt:=false; raumPunkt:=rp; fillchar(werte,sizeof(werte),#0); setlength(werte,anzTeilchen); @@ -601,8 +631,21 @@ begin nachbarn[1,false]:=linksImpuls; nachbarn[1,true]:=nil; for i:=0 to 1 do - if assigned(nachbarn[i,false]) then + 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 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 + nachbarn[i,b]:=np; + end; end; destructor tImpulsPunkt.destroy; @@ -612,7 +655,8 @@ var begin for i:=0 to 1 do for b:=false to true do - if assigned(nachbarn[i,b]) then + if assigned(nachbarn[i,b]) and + not nachbarn[i,b].istNullPunkt then nachbarn[i,b].nachbarn[i,not b]:=nil; setlength(werte,0); setlength(iMGamma,0); @@ -726,7 +770,9 @@ begin ( raumpunkt.emFelder[efEX,false] + raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]); werte[i,true]:= - tX * ( + 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] + @@ -741,7 +787,7 @@ begin nachbarn[1,false].werte[i,false] - 2*werte[i,false] ) - )* iDPX/2; + )* iDPX/2; *) end; // Der jeweils zweite Summand entspricht (in etwa) einer thermischen @@ -749,6 +795,50 @@ begin // 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 ****************************************************************** constructor tRaumPunkt.create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aP: longint; _dP: extended); @@ -871,7 +961,7 @@ var i: longint; begin // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f - for i:=1 to aP-2 do // Rand wird ausgelassen! + for i:=0 to aP-1 do phasenraum[i].berechneAbleitungen((i-(aP div 2))*dP,iDX,1/dP); end; diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps index 9a7b886..9330883 100644 --- a/Plasmapropagation.lps +++ b/Plasmapropagation.lps @@ -15,13 +15,14 @@ + - - - - + + + + - + @@ -30,16 +31,15 @@ - + - - + @@ -49,14 +49,14 @@ - + - + @@ -155,83 +155,126 @@ - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -- cgit v1.2.3-70-g09d2