From e425856940fd5237241d315c26d61b5f41e23ad4 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Tue, 4 Aug 2015 13:53:01 +0200 Subject: mfRho entfernt, Ausgabe gepuffert --- Physikunit.pas | 116 +++++++++++++++++++++++++++++--------------------- Plasmapropagation.lps | 91 +++++++++++++++++++-------------------- input.plap | 6 +-- 3 files changed, 117 insertions(+), 96 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 diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps index 98cc3b8..e959ac2 100644 --- a/Plasmapropagation.lps +++ b/Plasmapropagation.lps @@ -7,20 +7,22 @@ - - + + + + - - - - + + + + - + @@ -28,17 +30,15 @@ - - - + + - - - + + @@ -47,14 +47,14 @@ - + - + @@ -137,122 +137,123 @@ - + - + - + - - + + - + - + - + - + - + - + - + - + - + - + - + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + diff --git a/input.plap b/input.plap index ab2e416..33082aa 100644 --- a/input.plap +++ b/input.plap @@ -7,9 +7,9 @@ allgemein # runge-Kutta-3/8 # runge-Kutta-4 # euler-Vorwärts - ortsschritt 10^-4 * λ - zeitschritt 10^-4 * T - minZeitschritt 10^-7 * T + ortsschritt 10^-3 * λ + zeitschritt 10^-3 * T + minZeitschritt 10^-6 * T zeit 20 * T diffusionsterm 1 # max. p/Lp !setze $breite: (5 * λ) -- cgit v1.2.3-54-g00ecf