From 3012d4ee51fa18d9dae7d54d0ee9fab3f0bb847c Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Tue, 4 Aug 2015 12:57:51 +0200 Subject: Untergrenze fuer Zeitschritt konfigurierbar gemacht --- Physikunit.pas | 724 ++++++++------------------------------------------------- 1 file changed, 92 insertions(+), 632 deletions(-) (limited to 'Physikunit.pas') diff --git a/Physikunit.pas b/Physikunit.pas index 2c6b8ce..25a84a5 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -23,12 +23,14 @@ type efDPhiDX ); tMaterieFeldInhalt = ( - mfN,mfDPsiDX, - mfP,mfPX,mfPY,mfPZ, - mfGamma,mfIGamma + mfN,mfDPsiDX, // Teilchendichte, Geschwindigkeitsdichte + mfRho, // Ladungsdichte + mfP,mfPX,mfPY,mfPZ, // Impuls + mfGamma,mfIGamma // rel. Gamma-Faktor ); tGitter = class; + tFelder = class; { tAusgabeDatei } @@ -70,6 +72,9 @@ type { tWertePunkt } tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen + private + besitzer: tFelder; + public matWerte: array of array[tMaterieFeldInhalt,boolean] of extended; // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen emWerte: array[tEMFeldInhalt,boolean] of extended; @@ -77,11 +82,12 @@ type // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! lN,rN: tWertePunkt; - constructor create(linkerNachbar: tWertePunkt; materien: longint); + constructor create(linkerNachbar: tWertePunkt; derBesitzer: tFelder); destructor destroy; override; procedure berechneGammaUndP; procedure berechneNAbleitung(iDX,pDNMax: extended); overload; procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload; + procedure berechneRhoAbleitung; procedure berechneAbleitungen(dX,iDX: extended); 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; @@ -152,11 +158,11 @@ type private prot: tProtokollant; zeitverfahren: tZeitverfahren; - pDNMaxDynamisch: boolean; + pDNMaxDynamisch,abbruch: boolean; dX,iDX,pDNMax,dTMaximum,xl,t: extended; kvs: tKnownValues; procedure diffusionsTermAnpassen(pro: tProtokollant); - function pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert? + function pruefeMaxDT(aF: longint; var mDT,dT: extended; dTMin: extended): boolean; // wurde verkleinert? public aktuelleFelder: longint; @@ -164,7 +170,7 @@ type constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); destructor destroy; override; - procedure iteriereSchritt(var dT: extended); + procedure iteriereSchritt(var dT: extended; dTMin: extended); function dumpErhaltungsgroessen: boolean; procedure berechne(was: char; teilchen: longint); end; @@ -175,7 +181,7 @@ type private prot: tProtokollant; gitter: tGitter; - dT,tEnde,sT,sDT,sDX: extended; + dT,tEnde,sT,sDT,sDX,dTMin: extended; fortschrittsAnzeige: boolean; ausgabeDateien: array of tAusgabeDatei; public @@ -194,7 +200,7 @@ const 'A','AX','AY','AZ','DAXDT','DAYDT','DAZDT','DPHIDX' ); matFeldNamen: array[tMaterieFeldInhalt] of string = ( - 'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA' + 'N','DPSIDX','RHO','P','PX','PY','PZ','GAMMA','IGAMMA' ); var @@ -458,7 +464,7 @@ end; { tWertePunkt } -constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint); +constructor tWertePunkt.create(linkerNachbar: tWertePunkt; derBesitzer: tFelder); var emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; @@ -466,15 +472,13 @@ var i: longint; begin inherited create; + besitzer:=derBesitzer; fillchar(matWerte,sizeof(matWerte),#0); - setlength(matWerte,materien); - for i:=0 to length(matWerte)-1 do begin + setlength(matWerte,length(besitzer.spezLadungen)); + for i:=0 to length(matWerte)-1 do for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do for abl:=false to true do matWerte[i,maF,abl]:=0; - matWerte[i,mfGamma,false]:=0; - matWerte[i,mfIGamma,false]:=0; - end; for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do for abl:=false to true do emWerte[emF,abl]:=0; @@ -499,9 +503,9 @@ var i: longint; begin for i:=0 to length(matWerte)-1 do begin - matWerte[i,mfPX,false]:= emWerte[efAX,false] + matWerte[i,mfDPsiDX,false]; - matWerte[i,mfPY,false]:= emWerte[efAY,false]; - matWerte[i,mfPZ,false]:= emWerte[efAZ,false]; + matWerte[i,mfPX,false]:= besitzer.spezLadungen[i] * emWerte[efAX,false] + matWerte[i,mfDPsiDX,false]; + matWerte[i,mfPY,false]:= besitzer.spezLadungen[i] * emWerte[efAY,false]; + matWerte[i,mfPZ,false]:= besitzer.spezLadungen[i] * emWerte[efAZ,false]; matWerte[i,mfGamma,false]:= sqrt( @@ -535,7 +539,7 @@ begin + 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. + // und soll die numerische Stabilität bis zu p * d(ln n)/dx von pDNMax gewährleisten. // Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit // alpha = pDNMax * dX end; @@ -560,6 +564,16 @@ begin end end; +procedure tWertePunkt.berechneRhoAbleitung; +var + i: longint; +begin + for i:=0 to length(matWerte)-1 do + // dRho/dt = dN/dT * q_spez + matWerte[i,mfRho,true]:= + besitzer.spezLadungen[i]*matWerte[i,mfN,true]; +end; + procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen var i: longint; @@ -567,17 +581,17 @@ begin // d2Psi/dxdt = dPhi/dx - dGamma/dx for i:=0 to length(matWerte)-1 do matWerte[i,mfDPsiDX,true]:= - emWerte[efDPhiDX,false] + besitzer.spezLadungen[i] * emWerte[efDPhiDX,false] - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; - // d3Phi/dx2dt = dn/dt ... hier muss integriert werden: - // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + * deltax + // d3Phi/dx2dt = dRho/dt ... hier muss integriert werden: + // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + * deltax emWerte[efDPhiDX,true]:= lN.emWerte[efDPhiDX,true]; for i:=0 to length(matWerte)-1 do emWerte[efDPhiDX,true]:= emWerte[efDPhiDX,true] - + (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2; + + (lN.matWerte[i,mfRho,true] + matWerte[i,mfRho,true])*dX/2; // d2A/dt2 = Laplace(A) - Nabla(phi) ... emWerte[efDAXDT,true]:= @@ -587,17 +601,17 @@ begin ( rN.emWerte[efAY,false] - 2*emWerte[efAY,false] + lN.emWerte[efAY,false] )*sqr(iDX); emWerte[efDAZDT,true]:= ( rN.emWerte[efAZ,false] - 2*emWerte[efAZ,false] + lN.emWerte[efAZ,false] )*sqr(iDX); - // ... - n/gamma p + // ... - Rho/gamma p for i:=0 to length(matWerte)-1 do begin emWerte[efDAXDT,true]:= emWerte[efDAXDT,true] - - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]); + - (matWerte[i,mfRho,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]); emWerte[efDAYDT,true]:= emWerte[efDAYDT,true] - - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]); + - (matWerte[i,mfRho,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]); emWerte[efDAZDT,true]:= emWerte[efDAZDT,true] - - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]); + - (matWerte[i,mfRho,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]); end; // dA/dt = dA/dt @@ -617,7 +631,6 @@ begin result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]); end; - { linearkombinationsspezifische Methoden von tWertePunkt und tFelder } {$INCLUDE linearkombination.inc} @@ -678,12 +691,13 @@ begin 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); + inhalt[0]:=tWertePunkt.create(nil,self); for i:=1 to length(inhalt)-1 do - inhalt[i]:=tWertePunkt.create(inhalt[i-1],matAnz); + inhalt[i]:=tWertePunkt.create(inhalt[i-1],self); for i:=0 to length(inhalt)-1 do for j:=0 to length(teilchen)-1 do begin inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs); + inhalt[i].matWerte[j,mfRho,false]:=teilchen[j].spezifischeLadung * inhalt[i].matWerte[j,mfN,false]; inhalt[i].matWerte[j,mfDPsiDX,false]:=0; end; for rechts:=false to true do begin @@ -755,6 +769,9 @@ begin inhalt[1].berechneNAbleitung(iDX,pDNMax,false); inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true); + for i:=1 to length(inhalt)-1 do + inhalt[i].berechneRhoAbleitung; + setzeRaender(dT,iDX,iDT); for i:=1 to length(inhalt)-2 do @@ -777,7 +794,8 @@ var i: longint; begin inherited create; - zeitverfahren:=ZV; + abbruch:=false; + zeitverfahren:=zv; kvs:=bekannteWerte; prot:=tProtokollant.create(protokollant,name); dTMaximum:=deltaT; @@ -849,7 +867,7 @@ begin end; end; -function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert? +function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended; dTMin: extended): boolean; // wurde verkleinert? begin mDT:=min(mDT,felders[aF].maxDT); // das maximale dT, welches bei den Momentanen Ableitungen nicht zu @@ -860,17 +878,18 @@ begin 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); + if dT=sT) and (sT