diff options
Diffstat (limited to 'Physikunit.pas')
| -rw-r--r-- | Physikunit.pas | 233 |
1 files changed, 177 insertions, 56 deletions
diff --git a/Physikunit.pas b/Physikunit.pas index 142ee4b..ab99170 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -8,6 +8,7 @@ unit Physikunit; {$DEFINE Zeitschrittueberwachung} {$DEFINE Dichteueberwachung} +{ $DEFINE negativeDichteueberwachung} interface @@ -25,11 +26,20 @@ type tMaterieFeldInhalt = ( mfN,mfDPsiDX, // Teilchendichte, Geschwindigkeitsdichte mfP,mfPX,mfPY,mfPZ, // Impuls - mfGamma,mfIGamma // rel. Gamma-Faktor + mfGamma,mfIGamma // rel. Gamma-Faktor ); +const + erstesEMFMitAbleitung = efAX; + letztesEMFMitAbleitung = efDAZDT; + erstesMatFMitAbleitung = mfN; + letztesMatFMitAbleitung = mfDPsiDX; + +type + tGitter = class; tFelder = class; + tSimulation = class; { tAusgabeDatei } @@ -81,15 +91,16 @@ type // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen emWerte: array[tEMFeldInhalt,boolean] of extended; // EM-Felder und deren Zeitableitungen - // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! + // A, dPhiDX, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! lN,rN: tWertePunkt; constructor create(linkerNachbar: tWertePunkt; derBesitzer: tFelder); destructor destroy; override; procedure berechneGammaUndP; + procedure berechneDPhiDX(dX: extended); procedure berechneNAbleitung(iDX,pDNMax: extended); overload; procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload; - procedure berechneAbleitungen(dX,iDX: extended); + procedure berechneAbleitungen(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; procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload; @@ -112,6 +123,8 @@ type procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; function maxDT: extended; + function dump(mitNachbarn: boolean): string; + procedure nichtnegativieren; end; { TFelder } @@ -119,7 +132,8 @@ type tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private matAnz: longint; - spezLadungen: tExtendedArray; + spezLadungen, + massen: tExtendedArray; lichters: array[boolean] of tMyStringlist; procedure setzeRaender(iDX: extended); public @@ -157,6 +171,7 @@ type tGitter = class(tObject) private + besitzer: tSimulation; prot: tProtokollant; zeitverfahren: tZeitverfahren; pDNMaxDynamisch,abbruch: boolean; @@ -170,7 +185,7 @@ type aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern - constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); + constructor create(derBesitzer: tSimulation; 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; dTMin: extended); function dumpErhaltungsgroessen: boolean; @@ -429,7 +444,7 @@ var i: longint; begin prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc'); - prot.schreibe(prefix+'q/m = '+floattostr(spezifischeLadung)+' e/me'); + prot.schreibe(prefix+'q/m = '+floattostr(-spezifischeLadung)+' e/me'); prot.schreibe(prefix+'n(x):'); for i:=0 to length(dichteStuecke)-1 do begin prot.schreibe(prefix+' '+dichteStuecke[i]); @@ -557,6 +572,29 @@ begin end; end; +procedure tWertePunkt.berechneDPhiDX(dX: extended); +var + i: longint; +begin + // d2Phi/dx2 = rho ... hier muss integriert werden: + // dPhi/dx = dPhi/dx(x- deltax) + <rho> * deltax + + if assigned(lN) then begin + emWerte[efDPhiDX,false]:=lN.emWerte[efDPhiDX,false]; + for i:=0 to length(matWerte)-1 do + emWerte[efDPhiDX,false]:= + emWerte[efDPhiDX,false] + + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,false] + matWerte[i,mfN,false]) * dX/2; + end + else begin + emWerte[efDPhiDX,false]:=0; + for i:=0 to length(matWerte)-1 do + emWerte[efDPhiDX,false]:= + emWerte[efDPhiDX,false] + + besitzer.spezLadungen[i] * matWerte[i,mfN,false] * dX/2; + end; +end; + procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen var i: longint; @@ -602,7 +640,7 @@ begin end end; -procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen +procedure tWertePunkt.berechneAbleitungen(iDX: extended); // Zeitableitungen berechnen var i: longint; begin @@ -610,16 +648,7 @@ begin for i:=0 to length(matWerte)-1 do matWerte[i,mfDPsiDX,true]:= besitzer.spezLadungen[i] * emWerte[efDPhiDX,false] - - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; - - // d3Phi/dx2dt = dRho/dt ... hier muss integriert werden: - // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + <dRho/dt> * deltax - emWerte[efDPhiDX,true]:= - lN.emWerte[efDPhiDX,true]; - for i:=0 to length(matWerte)-1 do - emWerte[efDPhiDX,true]:= - emWerte[efDPhiDX,true] - + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2; +;// - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; // d2A/dt2 = Laplace(A) - Nabla(phi) ... emWerte[efDAXDT,true]:= @@ -659,6 +688,73 @@ begin result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]); end; +function tWertePunkt.dump(mitNachbarn: boolean): string; +var + i: longint; +begin + if mitNachbarn and assigned(lN) then result:=lN.dump(false)+' << ' + else result:=''; + + for i:=0 to length(matWerte)-1 do + result:= + result+ + 'N['+ inttostr(i)+'] = '+floattostr(matWerte[i,mfN,false])+' '+ + 'N''['+inttostr(i)+'] = '+floattostr(matWerte[i,mfN,true])+' '; + + result:=result+'maxDT: '+floattostr(maxDT); + + if mitNachbarn and assigned(rN) then result:=result+' >> '+rN.dump(false); +end; + +procedure tWertePunkt.nichtnegativieren; // Dichten nicht negativ machen +var + i: longint; + defizit: extended; + li,re: tWertePunkt; + pro: tProtokollant; +begin + for i:=0 to length(matWerte)-1 do + if matWerte[i,mfN,false]<0 then begin + defizit:=-matWerte[i,mfN,false]; + matWerte[i,mfN,false]:=0; + li:=self.lN; + re:=self.rN; + repeat + if assigned(li) then begin + if li.matWerte[i,mfN,false]>0 then begin + if defizit>=li.matWerte[i,mfN,false] then begin + defizit:=defizit-li.matWerte[i,mfN,false]; + li.matWerte[i,mfN,false]:=0; + end + else begin + li.matWerte[i,mfN,false]:=li.matWerte[i,mfN,false]-defizit; + defizit:=0;; + end; + end; + li:=li.lN; + end; + if assigned(re) then begin + if re.matWerte[i,mfN,false]>0 then begin + if defizit>=re.matWerte[i,mfN,false] then begin + defizit:=defizit-re.matWerte[i,mfN,false]; + re.matWerte[i,mfN,false]:=0; + end + else begin + re.matWerte[i,mfN,false]:=re.matWerte[i,mfN,false]-defizit; + defizit:=0;; + end; + end; + re:=re.rN; + end; + until (defizit<=0) or not (assigned(li) or assigned(re)); + if defizit>0 then begin + pro:=tProtokollant.create(besitzer.par.prot,'WertePunkt.nichtnegativieren'); + pro.schreibe('Ich konnte ein Teilchendefizit nicht auffüllen!',true); + besitzer.par.abbrechen; + end; + end; +end; + { linearkombinationsspezifische Methoden von tWertePunkt und tFelder } {$INCLUDE linearkombination.inc} @@ -722,10 +818,15 @@ begin inhalt[0]:=tWertePunkt.create(nil,self); for i:=1 to length(inhalt)-1 do inhalt[i]:=tWertePunkt.create(inhalt[i-1],self); + fillchar(massen,sizeof(massen),#0); + setlength(massen,length(teilchen)); + for i:=0 to length(massen)-1 do + massen[i]:=0; 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,mfDPsiDX,false]:=0; + massen[j]:=massen[j]+inhalt[i].matWerte[j,mfN,false]*par.dX; end; for rechts:=false to true do begin lichters[rechts]:=tMyStringlist.create(nil,''); @@ -778,6 +879,8 @@ var begin for i:=0 to length(inhalt)-1 do inhalt[i].berechneGammaUndP; + for i:=0 to length(inhalt)-1 do + inhalt[i].berechneDPhiDX(dX); for i:=0 to matAnz-1 do begin // Teilchen werden reflektiert inhalt[1].matWerte[i,mfPX,false]:= abs(inhalt[1].matWerte[i,mfPX,false]); @@ -794,7 +897,7 @@ begin setzeRaender(iDX); for i:=1 to length(inhalt)-2 do - inhalt[i].berechneAbleitungen(dX,iDX); + inhalt[i].berechneAbleitungen(iDX); end; function tFelder.maxDT: extended; @@ -808,12 +911,13 @@ end; { tGitter } -constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); +constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); var i: longint; begin inherited create; abbruch:=false; + besitzer:=derBesitzer; zeitverfahren:=zv; kvs:=bekannteWerte; prot:=tProtokollant.create(protokollant,name); @@ -859,17 +963,21 @@ end; procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant); var - i,j: longint; - curMax: extended; + i,j,k: longint; + curMax,tmp: extended; begin curMax:=0; + k:=0; for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do - for j:=0 to felders[aktuelleFelder].matAnz-1 do - 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-100))); + for j:=0 to felders[aktuelleFelder].matAnz-1 do begin + tmp:=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-100)); + if tmp>curMax then begin + curMax:=tmp; + k:=i; + end; + end; if curMax > pDNMax then begin if pDNMaxDynamisch then begin pDNMax:=2*curMax; @@ -878,10 +986,10 @@ begin 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); + 'T, x='+floattostr(xL+k*dX)+'), die numerische Stabilität ist nicht mehr gewährleistet, Simulation wird abgebrochen!',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 + abbrechen; end; end; end; @@ -892,18 +1000,17 @@ begin // das maximale dT, welches bei den momentanen Ableitungen nicht zu // unphysikalischen Effekten (z.B. negativen Dichten) führt - result:=dTMax<dT*2; + result:=(dTMax<dT) and (dT>1.1*dTMin); if result then begin // beabsichtigter Zeitschritt ist zu groß, dT:=dTMax/4; // dann machen wir ihn doch kleiner! + if dT<dTMin then + dT:=dTMin {$IFDEF Zeitschrittueberwachung} - 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); - abbrechen; - end; - prot.spuelen; - {$ENDIF} - exit; + else begin + prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')'); + prot.spuelen; + end + {$ENDIF}; end; end; @@ -928,11 +1035,8 @@ begin end; procedure tGitter.iteriereSchritt(var dT: extended; dTMin: extended); -var i: longint; - {$IFDEF Zeitschrittueberwachung} - j: longint; - {$ENDIF} - {$IFDEF Dichteueberwachung} +var {$IFDEF negativeDichteueberwachung} + i,j: longint; pro: tProtokollant; {$ENDIF} dTMax: extended; @@ -973,17 +1077,18 @@ begin break; until abbruch; - prot.schreibe('dT = '+floattostr(dT)); + aktuelleFelder:=1-aktuelleFelder; - {$IFDEF Dichteueberwachung} + {$IFDEF negativeDichteueberwachung} for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do for j:=0 to felders[aktuelleFelder].matAnz-1 do - if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then + if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then begin + pro:=tProtokollant.create(prot,'iteriereSchritt'); + pro.schreibe('n<0 bei:',true); + pro.schreibe(' t = ' +floattostr(t)+',',true); + pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); + pro.schreibe(' alte Werte:',true); with felders[1-aktuelleFelder] do begin - pro:=tProtokollant.create(prot,'iteriereSchritt'); - pro.schreibe('n<0 bei:',true); - pro.schreibe(' t = ' +floattostr(t)+',',true); - pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); @@ -996,9 +1101,25 @@ begin inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); - pro.destroyall; - halt(1); end; + pro.schreibe(' neue Werte:',true); + with felders[aktuelleFelder] do begin + pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* + inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); + pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); + pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true); + pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]* + inhalt[i].matWerte[j,mfIGamma,false])+',',true); + pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true); + pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true); + pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]* + inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); + pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); + pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); + end; + pro.free; + abbrechen; + end; {$ENDIF} t:=t+dT; dT:=max(dT,dTMax/100); @@ -1025,9 +1146,9 @@ begin ns[i]:=ns[i]*dX; s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]); {$IFDEF Dichteueberwachung} - if ns[i]>1000 then begin + if (not abbruch) and (abs(ns[i]-felders[aktuelleFelder].massen[i])>0.5) then begin result:=false; - pro.schreibe(' n['+inttostr(i)+'] > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true); + pro.schreibe(' n['+inttostr(i)+'] ('+floattostr(ns[i])+') unterscheidet sich stark vom Anfangswert ('+floattostr(felders[aktuelleFelder].massen[i])+'), es scheinen sehr viele Teilchen entstanden oder verloren gegangen zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true); end; {$ENDIF} end; @@ -1276,7 +1397,7 @@ begin end; if s='teilchen'+inttostr(length(teilchen))+'Ende' then break; if startetMit('spezifische Ladung ',s) then begin - teilchen[length(teilchen)-1].spezifischeLadung:=exprToFloat(false,s,kvs,nil); + teilchen[length(teilchen)-1].spezifischeLadung:=-exprToFloat(false,s,kvs,nil); kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung); continue; end; @@ -1385,7 +1506,7 @@ begin setlength(sighupSimulationen,length(sighupSimulationen)+1); sighupSimulationen[length(sighupSimulationen)-1]:=self; - gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); + gitter:=tGitter.create(self,round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); |