summaryrefslogtreecommitdiff
path: root/Physikunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'Physikunit.pas')
-rw-r--r--Physikunit.pas116
1 files changed, 68 insertions, 48 deletions
diff --git a/Physikunit.pas b/Physikunit.pas
index 25a84a5..0a2ca5e 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -24,7 +24,6 @@ type
);
tMaterieFeldInhalt = (
mfN,mfDPsiDX, // Teilchendichte, Geschwindigkeitsdichte
- mfRho, // Ladungsdichte
mfP,mfPX,mfPY,mfPZ, // Impuls
mfGamma,mfIGamma // rel. Gamma-Faktor
);
@@ -43,6 +42,9 @@ type
pre,suf: string;
datei: file;
pro: tProtokollant;
+ buf: pointer;
+ bufLen,bufPos: longint;
+ procedure schreibeBuffer(minBufLen: longint);
public
zeitAnz: longint;
sDT: extended;
@@ -87,7 +89,6 @@ type
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;
@@ -200,8 +201,9 @@ const
'A','AX','AY','AZ','DAXDT','DAYDT','DAZDT','DPHIDX'
);
matFeldNamen: array[tMaterieFeldInhalt] of string = (
- 'N','DPSIDX','RHO','P','PX','PY','PZ','GAMMA','IGAMMA'
+ 'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA'
);
+ minAusgabeBuffer = 1024*1024;
var
sighupSimulationen: array of tSimulation;
@@ -268,12 +270,19 @@ begin
if ableitung then
pre:='d'+pre;
+ bufLen:=0;
+ buf:=nil;
+ bufPos:=0;
+
pre:=prefix+pre;
suf:=suffix;
end;
destructor tAusgabeDatei.destroy;
begin
+ if bufPos<>0 then
+ schreibeBuffer(0);
+ freeMem(buf,bufLen);
if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin
pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true);
pro.destroyall;
@@ -282,6 +291,20 @@ begin
inherited destroy;
end;
+procedure tAusgabeDatei.schreibeBuffer(minBufLen: longint);
+begin
+ if buf=nil then begin
+ bufLen:=max(minAusgabeBuffer,minBufLen);
+ getmem(buf,bufLen);
+ end
+ else begin
+ reset(datei,1);
+ seek(datei,fileSize(datei));
+ blockWrite(datei,buf^,bufPos);
+ end;
+ bufPos:=0;
+end;
+
function tAusgabeDatei.dump: string;
begin
if teilchen<0 then
@@ -300,6 +323,8 @@ begin
pro.destroyall;
halt(1);
end;
+ if bufPos<>0 then
+ schreibeBuffer(0);
tCnt:=0;
assignfile(datei,pre+'-'+inttostr(nNum)+suf);
inc(nNum);
@@ -331,17 +356,29 @@ begin
sX:=cX+(cnt-1)*sDX;
inc(tCnt);
- reset(datei,1);
- seek(datei,fileSize(datei));
- blockWrite(datei,cX,sizeof(extended));
- blockWrite(datei,sX,sizeof(extended));
- blockWrite(datei,cnt,sizeof(integer));
+ if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(extended) then
+ schreibeBuffer(sizeof(integer) + (2+cnt)*sizeof(extended));
+
+ if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(extended) then begin
+ pro.schreibe('Schreibbuffer ist zu klein ('+inttostr(bufLen)+' statt mindestens '+inttostr(sizeof(integer) + (2+cnt)*sizeof(extended))+' Bytes)!');
+ pro.destroyall;
+ halt(1);
+ end;
+
+ move(cX,(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
+ move(sX,(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
+ move(cnt,(buf+bufPos)^,sizeof(integer));
+ inc(bufPos,sizeof(integer));
+
sX:=cX-Min(gitter.dX,sDX)/2;
if teilchen<0 then begin
for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
while cX>=sX do begin
dec(cnt);
- blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(extended));
+ move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
sX:=sX+sDX;
end;
cX:=cX+gitter.dX;
@@ -351,13 +388,13 @@ begin
for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
while cX>=sX do begin
dec(cnt);
- blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(extended));
+ move(gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
sX:=sX+sDX;
end;
cX:=cX+gitter.dX;
end;
end;
- closeFile(datei);
if cnt<>0 then begin
pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true);
pro.destroyall;
@@ -535,13 +572,13 @@ begin
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 + lN.matWerte[i,mfPX,false])
- - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2;
+ ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false] * iDX)
+ + lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false] * iDX)
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) / 2;
// Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung
// 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
+ // alpha = pDNMax * dX^2
end;
procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean);
@@ -552,28 +589,18 @@ begin
for i:=0 to length(matWerte)-1 do
// rechter Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
matWerte[i,mfN,true]:=
- ( lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false])
- - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax - matWerte[i,mfPX,false])) * iDX/2;
+ ( lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false] * iDX)
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax - matWerte[i,mfPX,false] * iDX)) / 2;
end
else begin
for i:=0 to length(matWerte)-1 do
// linker Randterm von 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])
- - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax + matWerte[i,mfPX,false])) * iDX/2;
+ ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false] * iDX)
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax + matWerte[i,mfPX,false] * iDX)) / 2;
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;
@@ -591,7 +618,7 @@ begin
for i:=0 to length(matWerte)-1 do
emWerte[efDPhiDX,true]:=
emWerte[efDPhiDX,true]
- + (lN.matWerte[i,mfRho,true] + matWerte[i,mfRho,true])*dX/2;
+ + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2;
// d2A/dt2 = Laplace(A) - Nabla(phi) ...
emWerte[efDAXDT,true]:=
@@ -605,13 +632,13 @@ begin
for i:=0 to length(matWerte)-1 do begin
emWerte[efDAXDT,true]:=
emWerte[efDAXDT,true]
- - (matWerte[i,mfRho,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]);
+ - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]);
emWerte[efDAYDT,true]:=
emWerte[efDAYDT,true]
- - (matWerte[i,mfRho,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]);
+ - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]);
emWerte[efDAZDT,true]:=
emWerte[efDAZDT,true]
- - (matWerte[i,mfRho,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]);
+ - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]);
end;
// dA/dt = dA/dt
@@ -697,7 +724,6 @@ begin
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
@@ -769,9 +795,6 @@ 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
@@ -884,6 +907,7 @@ begin
// prot.destroyall;
// halt(1);
end;
+ prot.spuelen;
{$ENDIF}
exit;
end;
@@ -994,7 +1018,7 @@ begin
end{of case};
break;
- until false;
+ until abbruch;
aktuelleFelder:=1-aktuelleFelder;
{$IFDEF Dichteueberwachung}
@@ -1028,31 +1052,24 @@ end;
function tGitter.dumpErhaltungsgroessen: boolean;
var i,j: integer;
- ns,rhos: tExtendedArray;
+ ns: tExtendedArray;
pro: tProtokollant;
s: string;
begin
pro:=tProtokollant.create(prot,'dumpErhaltungsgroessen');
setlength(ns,felders[aktuelleFelder].matAnz);
- setlength(rhos,length(ns));
result:=true;
s:='';
- for i:=0 to length(ns)-1 do begin
+ for i:=0 to length(ns)-1 do
ns[i]:=0;
- rhos[i]:=0;
- end;
for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
- for j:=0 to length(ns)-1 do begin
+ for j:=0 to length(ns)-1 do
ns[j]:= ns[j] + felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false];
- rhos[j]:=rhos[j] + felders[aktuelleFelder].inhalt[i].matWerte[j,mfRho,false];
- end;
for i:=0 to length(ns)-1 do begin
ns[i]:=ns[i]*dX;
- rhos[i]:=rhos[i]*dX;
s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]);
- s:=s+ ' rho['+inttostr(i)+']='+floattostr(rhos[i]);
{$IFDEF Dichteueberwachung}
if ns[i]>1000 then begin
result:=false;
@@ -1425,6 +1442,9 @@ destructor tSimulation.destroy;
var
i,j: longint;
begin
+ for i:=0 to length(ausgabeDateien)-1 do
+ ausgabeDateien[i].free;
+ setlength(ausgabeDateien,0);
gitter.free;
prot.free;
for i:=length(sighupSimulationen)-1 downto 0 do