diff options
Diffstat (limited to 'Physikunit.pas')
| -rw-r--r-- | Physikunit.pas | 1016 |
1 files changed, 633 insertions, 383 deletions
diff --git a/Physikunit.pas b/Physikunit.pas index 71e6a79..492866c 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -2,75 +2,99 @@ unit Physikunit; {$mode objfpc}{$H+} -{ $DEFINE Zeitschrittueberwachung} -{ $DEFINE Dichteueberwachung} +{$DEFINE Zeitschrittueberwachung} +{$DEFINE Dichteueberwachung} interface uses - Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit; + Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, crt; 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 + tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier); + tVerteilungsfunktion = function(x: extended): extended; + tEMFeldInhalt = ( + efA,efAX,efAY,efAZ, + efDAXDT,efDAYDT,efDAZDT, + efDPhiDX + ); + tMaterieFeldInhalt = ( + mfN,mfDPsiDX, + mfP,mfPX,mfPY,mfPZ, + mfGamma,mfIGamma + ); + + { tAusgabeDatei } + + tAusgabeDatei = class ableitung: boolean; - inhalt: tFeldInhalt; + emInhalt: tEMFeldInhalt; + matInhalt: tMaterieFeldInhalt; + teilchen: longint; // wenn < 0, dann EM-Feld name: string; datei: file; + constructor create(feldName,prefix,suffix: string; prot: tProtokollant); + destructor destroy; override; end; - tFelder = class; + tGitter = class; + + { tWertePunkt } - { tWerteGitter } + tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen + matWerte: array of array[tMaterieFeldInhalt,boolean] of extended; + // 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! + lN,rN: tWertePunkt; - tWerteGitter = class(tObject) // repräsentiert ein Gitter von Werten und deren Zeitableitungen - werte: array[boolean] of array of extended; - par: tFelder; + constructor create(linkerNachbar: tWertePunkt; materien: longint); + destructor destroy; override; + procedure berechneGammaUndP; + procedure berechneNAbleitung(iDX,pDNMax: extended); + procedure berechneAbleitungen(dT,dX,iDX: extended); + procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt + procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt end; { TFelder } tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen + private + matAnz: longint; + lichters: array[boolean] of tMyStringlist; + procedure setzeRaender(iDX: extended); public - materieFelder: array of tWerteGitter; // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen - elektroMagnetischeFelder: tWerteGitter; // EM-Felder und deren Zeitableitungen - // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! - dX,iDX,x,pDNMax: extended; + inhalt: array of tWertePunkt; + par: tGitter; - constructor create(deltaX,pDNMa: extended); - procedure berechneGammaUndP; - procedure berechneNAbleitung; - procedure berechneAbleitungen(dT: extended); + constructor create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter); + destructor destroy; override; + procedure berechneAbleitungen(dT,dX,iDX,pDNMax: extended); + procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt + procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt 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); + prot: tProtokollant; + zeitverfahren: tZeitverfahren; + pDNMaxDynamisch: boolean; + dX,iDX,pDNMax,dTMaximum,xl,t: extended; + kvs: tKnownValues; + procedure diffusionsTermAnpassen(pro: tProtokollant); public + aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern - dTMaximum,dX,iDX,xl,xr,t,pDNMax: extended; - pDNMaxDynamisch: boolean; - Zeitverfahren: TZeitverfahren; - Al,Ar: TVerteilungsfunktion; - constructor create(size: longint; deltaT,deltaX,pDNMa: extended; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); + + constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); destructor destroy; override; - procedure iteriereSchritt(var dT: extended); - procedure macheAusgabe(AD: TAusgabedatei; sDX: extended); + procedure iteriereSchritt(dT: extended); + procedure macheAusgabe(aD: tAusgabedatei; sDX: extended); function gibErhaltungsgroessen: string; end; @@ -78,59 +102,157 @@ type tSimulation = class(tObject) private - prot: tProtokollant; - gitter: tGitter; - kvs: tKnownValues; + prot: tProtokollant; + gitter: tGitter; + dT,tEnde,sT,sDT,sDX: extended; + fortschrittsAnzeige: boolean; + ausgabeDateien: array of tAusgabeDatei; public constructor create(inName: string; Protokollant: tProtokollant); destructor destroy; override; + function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? end; +implementation + 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'; + emFeldNamen: array[tEMFeldInhalt] of string = ( + 'A','AX','AY','AZ','DAXDT','DAYDT','DAZDT','DPHIDX' + ); + matFeldNamen: array[tMaterieFeldInhalt] of string = ( + 'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA' + ); -implementation +{ tAusgabeDatei } + +constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant); +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + num: longint; + abl,gef: boolean; +begin + inherited create; + num:=length(feldName); + while (num>0) and (feldName[num] in ['0'..'9']) do + dec(num); + inc(num); + if num<=length(feldName) then begin + teilchen:=strtoint(copy(feldName,num,length(feldName))); + delete(feldName,num,length(feldName)); + end + else + teilchen:=0; + + emInhalt:=efA; + matInhalt:=mfN; + feldName:=ansiUpperCase(feldName); + + gef:=false; + for abl:=false to true do begin + if not gef then + for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do + if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin + emInhalt:=emF; + teilchen:=-1; + gef:=true; + break; + end; + if not gef then + for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do + if copy('D',1,byte(abl))+matFeldNamen[maF] = feldName then begin + matInhalt:=maF; + gef:=true; + break; + end; + end; + + if not gef then begin + prot.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true); + halt(1); + end; + + if teilchen>=0 then + name:=matFeldNamen[maF]+inttostr(teilchen) + else + name:=emFeldNamen[emF]; -{ TFeld } + if abl then + name:='d'+name; + name:=prefix+name+suffix; -constructor TFeld.create(deltaX,pDNMa: extended); -var FI: TFeldinhalt; - Abl: Boolean; + assignFile(datei,name); + rewrite(datei); +end; + +destructor tAusgabeDatei.destroy; +begin + closeFile(datei); + inherited destroy; +end; + +{ tWertePunkt } + +constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint); +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + abl: boolean; + i: longint; begin inherited create; - lN:=nil; + fillchar(matWerte,sizeof(matWerte),#0); + setlength(matWerte,materien); + for i:=0 to length(matWerte)-1 do begin + 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; + lN:=linkerNachbar; 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; + if assigned(lN) then + lN.rN:=self; +end; + +destructor tWertePunkt.destroy; +begin + setlength(matWerte,0); + if assigned(lN) then + lN.rN:=nil; + if assigned(rN) then + rN.lN:=nil; + inherited destroy; end; -procedure TFeld.berechneGammaUndP; // aus aktuellem dPsi/dX und A -var b: byte; - tmp: extended; +procedure tWertePunkt.berechneGammaUndP; // aus aktuellem dPsi/dX und A +var + i: longint; 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)]); + 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,mfGamma,false]:= + sqrt( + 1 + + sqr(matWerte[i,mfPX,false]) + + sqr(matWerte[i,mfPY,false]) + + sqr(matWerte[i,mfPZ,false]) + ); + matWerte[i,mfIGamma,false]:= + 1/matWerte[i,mfGamma,false]; end; - Groessen[false,fiGamma]:=sqrt(1+tmp); - Groessen[false,fiiGamma]:=1/Groessen[false,fiGamma]; end; -procedure TFeld.berechneNAbleitung; // Zeitableitung von n berechnen +procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen +var + i: longint; begin (* // dn/dt = -d(n/Gamma p_x)/dx Groessen[true,fiN]:= @@ -141,104 +263,226 @@ begin )*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 + for i:=0 to length(matWerte)-1 do + // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n + matWerte[i,mfN,true]:= + ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false]) + - ln.matWerte[i,mfN,false] * ln.matWerte[i,mfIGamma,false] * (-pDNMax - rN.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. + // 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; +procedure tWertePunkt.berechneAbleitungen(dT,dX,iDX: extended); // Zeitableitungen berechnen +var + i: longint; begin // d2Psi/dxdt = dPhi/dx - dGamma/dx - Groessen[true,fidPsiDX]:= - Groessen[false,fidPhiDX] - - (rN^.Groessen[false,fiGamma] - lN^.Groessen[false,fiGamma])*iDX/2; + for i:=0 to length(matWerte)-1 do + matWerte[i,mfDPsiDX,true]:= + 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) + <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) + 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; + + // d2A/dt2 = Laplace(A) - Nabla(phi) ... + emWerte[efDAXDT,true]:= + ( rN.emWerte[efAX,false] - 2*emWerte[efAX,false] + lN.emWerte[efAX,false] )*sqr(iDX) + - emWerte[efDPhiDX,false]; + emWerte[efDAYDT,true]:= + ( 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 + 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]); + emWerte[efDAYDT,true]:= + emWerte[efDAYDT,true] + - (matWerte[i,mfN,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]); + end; + // 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; + emWerte[efAX,true]:=emWerte[efDAXDT,false] + emWerte[efDAXDT,true]*dT; end; -{ TGitter } +procedure tWertePunkt.liKo(in1,in2: tWertePunkt; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + i: longint; +begin +(* tEMFeldInhalt = ( + efA,efAX,efAY,efAZ, + efDAXDT,efDAYDT,efDAZDT, + efDPhiDX + ); *) + for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde + emWerte[emF,false]:= in1.emWerte[emF,false] + fak2 * in2.emWerte[emF,true]; + +(* tMaterieFeldInhalt = ( + mfN,mfDPsiDX, + mfP,mfPX,mfPY,mfPZ, + mfGamma,mfIGamma + ); *) + for i:=0 to length(matWerte)-1 do // siehe oben + for maF:=mfN to mfDPsiDX do + matWerte[i,maF,false]:= in1.matWerte[i,maF,false] + fak2 * in2.matWerte[i,maF,true]; +end; -procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended; out dTMax: extended); -var i: longint; +procedure tWertePunkt.liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + 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; +(* tEMFeldInhalt = ( + efA,efAX,efAY,efAZ, + efDAXDT,efDAYDT,efDAZDT, + efDPhiDX + ); *) + for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde + emWerte[emF,false]:= + in1.emWerte[emF,false] + + fak2 * in2.emWerte[emF,true] + + fak3 * in3.emWerte[emF,true] + + fak4 * in4.emWerte[emF,true] + + fak5 * in5.emWerte[emF,true]; + +(* tMaterieFeldInhalt = ( + mfN,mfDPsiDX, + mfP,mfPX,mfPY,mfPZ, + mfGamma,mfIGamma + ); *) + for i:=0 to length(matWerte)-1 do // siehe oben + for maF:=mfN to mfDPsiDX do + matWerte[i,maF,false]:= + in1.matWerte[i,maF,false] + + fak2 * in2.matWerte[i,maF,true] + + fak3 * in3.matWerte[i,maF,true] + + fak4 * in4.matWerte[i,maF,true] + + fak5 * in5.matWerte[i,maF,true]; end; -procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended); -var i: longint; +{ tFelder } + +constructor tFelder.create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter); +var + i,j: longint; + rechts: boolean; 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; + inherited create; + par:=parent; + matAnz:=dichten.count; + fillchar(inhalt,sizeof(inhalt),#0); + setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen + inhalt[0]:=tWertePunkt.create(nil,matAnz); + for i:=1 to length(inhalt)-1 do + inhalt[i]:=tWertePunkt.create(inhalt[i-1],matAnz); + for i:=0 to length(inhalt)-1 do begin + parent.kvs.add('x',par.xl+(i-2)/groesse*par.dX); + for j:=0 to dichten.count-1 do + inhalt[i].matWerte[j,mfN,false]:=exprToFloat(false,dichten[j],parent.kvs,nil); + end; + parent.kvs.rem('x'); + for rechts:=false to true do begin + lichters[rechts]:=tMyStringlist.create(nil); + lichters[rechts].text:=lichter.text; 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); + lichters[false].grep('^links .*'); + lichters[false].subst('^links *',''); + lichters[true].grep('^rechts .*'); + lichters[true].subst('^rechts *',''); end; -procedure TGitter.setzeRaender(Felder: longint); -var FI: TFeldInhalt; +destructor tFelder.destroy; 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] + lichters[false].free; + lichters[true].free; + inherited destroy; +end; + +procedure tFelder.setzeRaender(iDX: extended); +var + emF: tEMFeldInhalt; +begin + for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld + inhalt[0].emWerte[emF,true]:= + (inhalt[1].emWerte[emF,false] - + inhalt[0].emWerte[emF,false])*iDX; + inhalt[length(inhalt)-1].emWerte[emF,true]:= + (inhalt[length(inhalt)-2].emWerte[emF,false] - + inhalt[length(inhalt)-1].emWerte[emF,false])*iDX; + end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler bzw. nicht beachtete numerische Dispersion) + + (* inhalt[0].emWerte[efAY,true]:= + inhalt[0].emWerte[fiAY,true] + (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; + + (Ar(t+dTMaximum)-Ar(t))/dTMaximum; *) end; -constructor TGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); -var i,j: longint; +procedure tFelder.berechneAbleitungen(dt,dx,iDX,pDNMax: extended); +var + i: longint; +begin + for i:=0 to length(inhalt)-1 do + inhalt[i].berechneGammaUndP; + for i:=0 to matAnz-1 do begin // Teilchen werden reflektiert + inhalt[1].matWerte[i,mfPX,false]:= + abs(inhalt[1].matWerte[i,mfPX,false]); + inhalt[length(inhalt)-2].matWerte[i,mfPX,false]:= + -abs(inhalt[length(inhalt)-2].matWerte[i,mfPX,false]); + end; + + setzeRaender(iDX); + + for i:=1 to length(inhalt)-2 do + inhalt[i].berechneNAbleitung(iDX,pDNMax); + for i:=1 to length(inhalt)-2 do + inhalt[i].berechneAbleitungen(dT,dX,iDX); +end; + +procedure tFelder.liKo(in1,in2: tFelder; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt +var + i: longint; +begin + for i:=0 to length(inhalt)-1 do + inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],fak2); +end; + +procedure tFelder.liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt +var + i: longint; +begin + for i:=0 to length(inhalt)-1 do + inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],in4.inhalt[i],in5.inhalt[i],fak2,fak3,fak4,fak5); +end; + +{ tGitter } + +constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); +var + i: longint; begin inherited create; - Ar:=@nullfunktion; - Al:=@nullfunktion; - Zeitverfahren:=ZV; - groesse:=size; - Prot:=TProtokollant.create(Protokollant,'TGitter'); + zeitverfahren:=ZV; + kvs:=bekannteWerte; + Prot:=TProtokollant.create(Protokollant,'tGitter'); dTMaximum:=deltaT; dX:=deltaX; iDX:=1/dX; @@ -252,252 +496,192 @@ begin 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; - knownValues.add('x',felders[i,j].x); - 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); + + for i:=0 to length(felders)-1 do + felders[i]:=tFelder.create(size,dichten,lichter,self); + aktuelleFelder:=0; t:=0; end; -destructor TGitter.destroy; -var i,j: longint; +destructor tGitter.destroy; +var + i: longint; begin - for i:=0 to length(Felders)-1 do begin - for j:=0 to length(Felders[i])-1 do - Felders[i,j].free; - setlength(Felders[i],0); - end; - setlength(Felders,0); + for i:=0 to length(Felders)-1 do + felders[i].free; + setlength(felders,0); inherited destroy; end; -procedure TGitter.iteriereSchritt(var dT: extended); +procedure tGitter.iteriereSchritt(dT: extended); var i: longint; - FI: TFeldInhalt; - dTMax: extended; {$IFDEF Zeitschrittueberwachung} - Pro: TProtokollant; + pro: tProtokollant; + j: longint; {$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} + felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDX,pDNMax); // y' = y'(t,y(t)) - 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; + case zeitverfahren of + zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt + felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT); 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; + felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2 + felders[2].berechneAbleitungen(dT/2,dX,iDX,pDNMax); // ya' = y'(t+dt/2,ya) + + felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2 + felders[3].berechneAbleitungen(dT/2,dX,iDX,pDNMax); // yb' = y'(t+dt/2,yb) + + felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt + felders[4].berechneAbleitungen(dT/2,dX,iDX,pDNMax); // yc' = y'(t+dt,yc) + + felders[1-aktuelleFelder].liKo( + felders[aktuelleFelder], + felders[aktuelleFelder], + felders[2], + felders[3], + felders[4], + dT/6, + dT/3, + dT/3, + dT/6 + ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') 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; + {$IFDEF Dichteueberwachung} + with felders[aktuelleFelder] do + for i:=0 to length(inhalt)-1 do + for j:=0 to matAnz-1 do + if inhalt[i].matWerte[j,mfN,false]<0 then begin + pro:=tProtokollant.create(prot,'iteriereSchritt'); + pro.schreibe( + 'n<=0 bei '+ + floattostr(t)+' '+ + inttostr(i)+' '+ + floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* + inhalt[i-1].matWerte[j,mfIGamma,false])+' '+ + floattostr(inhalt[i-1].matWerte[j,mfN,false])+' '+ + floattostr(inhalt[i-1].matWerte[j,mfN,true])+' '+ + floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]* + inhalt[i].matWerte[j,mfIGamma,false])+' '+ + floattostr(inhalt[i].matWerte[j,mfN,false])+' '+ + floattostr(inhalt[i].matWerte[j,mfN,true])+' '+ + floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]* + inhalt[i+1].matWerte[j,mfIGamma,false])+' '+ + floattostr(inhalt[i+1].matWerte[j,mfN,false])+' '+ + floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); + pro.destroyall; + halt(1); + end; {$ENDIF} end; -procedure TGitter.macheAusgabe(AD: TAusgabedatei; sDX: extended); -var i: longint; - tmp,sX: extended; - b: byte; +procedure tGitter.macheAusgabe(aD: tAusgabedatei; sDX: extended); +var i: longint; + sX,cX,tmp: extended; + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; begin - if AD.Inhalt in [fiA,fiP] then - for i:=0 to length(Felders[aktuelleFelder])-1 do begin + if (aD.teilchen<0) and (aD.emInhalt=efA) then + for i:=0 to length(felders[aktuelleFelder].inhalt)-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); + for emF:=efAX to efAZ do + tmp:=tmp+sqr(felders[aktuelleFelder].inhalt[i].emWerte[emF,false]); + felders[aktuelleFelder].inhalt[i].emWerte[efA,false]:=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)); + if (aD.teilchen>=0) and (aD.matInhalt=mfP) then + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + tmp:=0; + for maF:=mfPX to mfPZ do + tmp:=tmp+sqr(felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,maF,false]); + felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,maF,false]:=sqrt(tmp); + end; + + i:=floor((length(felders[aktuelleFelder].inhalt)-1)/sDX*dX+1); + sX:=xl+(i-1)*sDX; + blockWrite(aD.datei,xl,sizeof(extended)); + blockWrite(aD.datei,sX,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; + cX:=xl; + if aD.teilchen<0 then begin + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + blockWrite(aD.datei,felders[aktuelleFelder].inhalt[i].emWerte[aD.emInhalt,aD.ableitung],sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+dX; end; + end + else begin + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + blockWrite(aD.datei,felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,aD.matInhalt,aD.ableitung],sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+dX; + end; + end; end; -function TGitter.gibErhaltungsgroessen: string; -var i,j,k: integer; - n: extended; - Pro: TProtokollant; +function tGitter.gibErhaltungsgroessen: string; +var i,j: integer; + ns: tExtendedArray; + 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; + pro:=tProtokollant.create(prot,'gibErhaltungsgroessen'); + setlength(ns,felders[aktuelleFelder].matAnz); + result:=''; + for i:=0 to length(ns)-1 do + ns[i]:=0; + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do + for j:=0 to length(ns)-1 do + ns[j]:=ns[j] + felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]; + for i:=0 to length(ns)-1 do begin + result:=result + ' n['+inttostr(i)+']='+floattostr(ns[i]); + {$IFDEF Dichteueberwachung} + if ns[i]>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; - 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)+')'); + {$ENDIF} end; - {$ENDIF} - Pro.free; + delete(result,1,1); + pro.free; +end; + +procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant); +var + i,j: longint; +begin + for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do + for j:=0 to felders[aktuelleFelder].matAnz-1 do + if 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-10)) + > pDNMax then begin + if pDNMaxDynamisch then begin + pDNMax:= + 2* + 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-10)); + 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; end; { tSimulation } @@ -507,11 +691,10 @@ var ifile: tMyStringlist; zeitverfahren: tZeitverfahren; s,t,aSuffix,aPrefix: string; - deltaX,deltaT,endzeit,breite,sDT,pDNMax: extended; - ausgabeDateien: array of tAusgabeDatei; - abl: boolean; - fi: tFeldinhalt; + deltaX,deltaT,breite,pDNMax: extended; i,j: longint; + kvs: tKnownValues; + dichten,lichter: tMyStringlist; begin inherited create; prot:=tProtokollant.create(Protokollant,'tSimulation'); @@ -527,25 +710,35 @@ begin // Standardeinstellungen Bereich 'allgemein' zeitverfahren:=zfRungeKuttaVier; deltaX:=1e-2; - kvs.add('λ',1/deltaX); + kvs.add('λ',1); + kvs.add('T',1); + kvs.add('dX',deltaX); deltaT:=-1; sDT:=-1; - endzeit:=100; + sDX:=-1; + tEnde:=100; breite:=10.0; pDNMax:=-1; + fortschrittsAnzeige:=false; // Standardeinstellungen Bereich 'ausgaben' aPrefix:=extractfilepath(paramstr(0)); aSuffix:='.dat'; setlength(ausgabeDateien,0); - while ifile.readln(s) do begin + // Standardeinstellungen Breich 'lichtVon...' + lichter:=tMyStringlist.create(prot); + + // Standardeinstellungen Breich 'teilchen...' + dichten:=tMyStringlist.create(prot); + + while ifile.readln(s) do begin writeln(s); if s='allgemein' then begin repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true); halt(1); - end; + end; writeln(s); if s='allgemeinEnde' then break; if s='runge-Kutta-4' then begin Zeitverfahren:=zfRungeKuttaVier; @@ -557,7 +750,7 @@ begin end; if startetMit('ortsschritt',s) then begin deltaX:=exprtofloat(false,s,kvs,nil); - kvs.add('λ',1/deltaX); + kvs.add('dX',deltaX); continue; end; if startetMit('zeitschritt',s) then begin @@ -570,13 +763,21 @@ begin continue; end; if startetMit('zeit',s) then begin - endzeit:=exprtofloat(false,s,kvs,nil); + tEnde:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('breite',s) then begin breite:=exprtofloat(false,s,kvs,nil); continue; end; + if s='mit Fortschrittsanzeige' then begin + fortschrittsAnzeige:=true; + continue; + end; + if s='ohne Fortschrittsanzeige' then begin + fortschrittsAnzeige:=false; + continue; + end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true); halt(1); until false; @@ -602,20 +803,16 @@ begin sDT:=exprtofloat(false,s,kvs,nil); continue; end; + if startetMit('ortsschritt',s) then begin + sDX:=exprtofloat(false,s,kvs,nil); + continue; + end; if startetMit('felder',s) then begin s:=s+','; while s<>'' do begin t:=erstesArgument(s,','); - for abl:=false to true do - for fI:=low(tFeldInhalt) to high(tFeldInhalt) do - if not (abl and (fI in [fiGamma,fiP])) then - if t=copy('d',1,byte(abl))+optionNamen[fI] then begin - setlength(ausgabedateien,length(ausgabedateien)+1); - ausgabedateien[length(ausgabedateien)-1].inhalt:=fI; - ausgabedateien[length(ausgabedateien)-1].ableitung:=abl; - ausgabedateien[length(ausgabedateien)-1].name:=aPrefix+t+aSuffix; - assignFile(ausgabedateien[length(ausgabedateien)-1].datei,ausgabedateien[length(ausgabedateien)-1].name); - end; + setlength(ausgabeDateien,length(ausgabeDateien)+1); + ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot); end; continue; end; @@ -638,26 +835,79 @@ begin deltaT:=deltaX/10; if sDT<0 then sDT:=deltaT; + if sDX<0 then + sDX:=deltaX; for i:=0 to length(ausgabedateien)-1 do begin rewrite(ausgabedateien[i].datei,1); j:=0; BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); - j:=floor(Endzeit/sDT+1); + j:=floor(tEnde/sDT+1); BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); end; ifile.free; - gitter:=tGitter.create(Breite,deltaT,deltaX,pDNMax,@Anfangsdichte,Zeitverfahren,Prot); + gitter:=tGitter.create(round(breite/deltaX),deltaT,deltaX,pDNMax,kvs,dichten,lichter,zeitverfahren,prot); + dichten.free; + lichter.free; end; destructor tSimulation.destroy; begin gitter.free; - kvs.free; prot.free; inherited destroy; end; +function tSimulation.iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? +var + i: longint; + c: char; +begin + result:=false; + + zeitPhysik:=zeitPhysik-now; + if errorCode<2 then + gitter.iteriereSchritt(dT); + zeitPhysik:=zeitPhysik+now; + zeitDatei:=zeitDatei-now; + while gitter.t>=sT do begin + sT:=sT+sDT; + for i:=0 to length(Ausgabedateien)-1 do + gitter.macheAusgabe(ausgabeDateien[i],sDX); + end; + zeitDatei:=zeitDatei+now; + + gitter.iteriereSchritt(dT); + + if fortschrittsAnzeige then begin + if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) or keypressed then begin + if keypressed then c:=readkey + else c:=#0; + case c of + #27,'q': begin + errorCode:=3; + exit; + end; + ' ': begin + writeln(' ... Pause (beliebige Taste drücken um fortzufahren) ...'); + readkey; + writeln(' ... weiter geht''s ...'); + end; + else begin + Prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true); + Prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true); + Prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true); + Prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true); + if errorCode<2 then + Prot.schreibe(gitter.gibErhaltungsgroessen); + end; + end{of case}; + end; + end; + + result:=gitter.t>=tEnde; +end; + end. |