diff options
-rw-r--r-- | Physikunit.pas | 116 | ||||
-rw-r--r-- | Plasmapropagation.lps | 91 | ||||
-rw-r--r-- | 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 @@ <Unit0> <Filename Value="Plasmapropagation.lpr"/> <IsPartOfProject Value="True"/> - <CursorPos X="34" Y="11"/> - <UsageCount Value="183"/> + <TopLine Value="54"/> + <CursorPos X="37" Y="34"/> + <UsageCount Value="184"/> <Loaded Value="True"/> </Unit0> <Unit1> <Filename Value="Physikunit.pas"/> <IsPartOfProject Value="True"/> + <IsVisibleTab Value="True"/> <EditorIndex Value="1"/> - <TopLine Value="467"/> - <CursorPos X="43" Y="544"/> - <FoldState Value=" T3nG083 piYj90A5 pjBjZ08114]KalH0V[K4KkM[c5lAlQ0b113 T0*3Q/0F17121Z"/> - <UsageCount Value="124"/> + <TopLine Value="205"/> + <CursorPos X="23" Y="327"/> + <FoldState Value=" T3le0,312 piXlY054 pj6jN08123]Keld0B[K6T08]IlBlR0U]lp6q10G2 T0*KQ/0N1211"/> + <UsageCount Value="125"/> <Bookmarks Count="1"> - <Item0 Y="996"/> + <Item0 Y="1020"/> </Bookmarks> <Loaded Value="True"/> </Unit1> @@ -28,17 +30,15 @@ <Filename Value="../units/protokollunit.pas"/> <IsPartOfProject Value="True"/> <EditorIndex Value="-1"/> - <TopLine Value="20"/> - <CursorPos X="22" Y="31"/> - <UsageCount Value="86"/> + <CursorPos X="15" Y="46"/> + <UsageCount Value="87"/> </Unit2> <Unit3> <Filename Value="input.plap"/> <IsPartOfProject Value="True"/> - <IsVisibleTab Value="True"/> <EditorIndex Value="3"/> - <CursorPos X="23" Y="12"/> - <UsageCount Value="85"/> + <CursorPos X="42" Y="7"/> + <UsageCount Value="86"/> <Loaded Value="True"/> <DefaultSyntaxHighlighter Value="None"/> </Unit3> @@ -47,14 +47,14 @@ <IsPartOfProject Value="True"/> <EditorIndex Value="2"/> <CursorPos X="48" Y="220"/> - <UsageCount Value="33"/> + <UsageCount Value="34"/> <Loaded Value="True"/> </Unit4> <Unit5> <Filename Value="input.epost"/> <EditorIndex Value="4"/> <CursorPos X="33" Y="21"/> - <UsageCount Value="75"/> + <UsageCount Value="76"/> <Loaded Value="True"/> <DefaultSyntaxHighlighter Value="None"/> </Unit5> @@ -137,122 +137,123 @@ <JumpHistory Count="30" HistoryIndex="29"> <Position1> <Filename Value="Physikunit.pas"/> - <Caret Line="600" Column="53" TopLine="567"/> + <Caret Line="89" Column="44" TopLine="68"/> </Position1> <Position2> <Filename Value="Physikunit.pas"/> - <Caret Line="603" Column="53" TopLine="571"/> + <Caret Line="130" Column="55" TopLine="97"/> </Position2> <Position3> <Filename Value="Physikunit.pas"/> - <Caret Line="606" Column="53" TopLine="574"/> + <Caret Line="162" Column="18" TopLine="129"/> </Position3> <Position4> - <Filename Value="Physikunit.pas"/> - <Caret Line="999" Column="80" TopLine="967"/> + <Filename Value="Plasmapropagation.lpr"/> + <Caret Line="34" Column="37"/> </Position4> <Position5> <Filename Value="Physikunit.pas"/> - <Caret Line="1003" Column="78" TopLine="971"/> + <Caret Line="27" TopLine="8"/> </Position5> <Position6> <Filename Value="Physikunit.pas"/> - <Caret Line="881" TopLine="648"/> + <Caret Line="565" TopLine="545"/> </Position6> <Position7> <Filename Value="Physikunit.pas"/> - <Caret TopLine="664"/> + <Caret Line="582" Column="36" TopLine="545"/> </Position7> <Position8> <Filename Value="Physikunit.pas"/> - <Caret Line="87" Column="38" TopLine="80"/> + <Caret Line="596" Column="51" TopLine="577"/> </Position8> <Position9> <Filename Value="Physikunit.pas"/> - <Caret Line="503" Column="46" TopLine="466"/> + <Caret Line="668" TopLine="668"/> </Position9> <Position10> <Filename Value="Physikunit.pas"/> - <Caret Line="569" TopLine="387"/> + <Caret Line="740" TopLine="712"/> </Position10> <Position11> <Filename Value="Physikunit.pas"/> - <Caret Line="90" Column="15" TopLine="79"/> + <Caret Line="1026" Column="14" TopLine="1013"/> </Position11> <Position12> <Filename Value="Physikunit.pas"/> - <Caret Line="504" Column="23" TopLine="467"/> + <Caret Line="1475" TopLine="1029"/> </Position12> <Position13> <Filename Value="Physikunit.pas"/> - <Caret Line="87" Column="32" TopLine="73"/> + <Caret Line="1025" Column="31" TopLine="876"/> </Position13> <Position14> <Filename Value="Physikunit.pas"/> - <Caret Line="572" Column="14" TopLine="550"/> + <Caret Line="284" Column="24" TopLine="201"/> </Position14> <Position15> <Filename Value="Physikunit.pas"/> - <Caret Line="684" TopLine="675"/> + <Caret Line="354" Column="32" TopLine="331"/> </Position15> <Position16> <Filename Value="Physikunit.pas"/> + <Caret Line="212" Column="81" TopLine="191"/> </Position16> <Position17> <Filename Value="Physikunit.pas"/> - <Caret Line="26" Column="8"/> + <Caret Line="51" Column="71" TopLine="27"/> </Position17> <Position18> <Filename Value="Physikunit.pas"/> - <Caret Line="232" Column="17" TopLine="205"/> + <Caret Line="272" Column="51" TopLine="242"/> </Position18> <Position19> <Filename Value="Physikunit.pas"/> - <Caret Line="579" Column="4" TopLine="464"/> + <Caret Line="1276" Column="156" TopLine="1262"/> </Position19> <Position20> <Filename Value="Physikunit.pas"/> - <Caret Line="637" TopLine="464"/> + <Caret Line="212" Column="66" TopLine="192"/> </Position20> <Position21> <Filename Value="Physikunit.pas"/> - <Caret Line="715" Column="16" TopLine="678"/> + <Caret Line="51" Column="60" TopLine="47"/> </Position21> <Position22> <Filename Value="Physikunit.pas"/> - <Caret Line="872" Column="8" TopLine="728"/> + <Caret Line="351" Column="20" TopLine="301"/> </Position22> <Position23> <Filename Value="Physikunit.pas"/> - <Caret Line="1031" Column="18" TopLine="716"/> + <Caret Line="2"/> </Position23> <Position24> <Filename Value="Physikunit.pas"/> - <Caret Line="1070" Column="7" TopLine="728"/> + <Caret Line="351" Column="73" TopLine="331"/> </Position24> <Position25> <Filename Value="Physikunit.pas"/> - <Caret Line="577" TopLine="491"/> + <Caret Line="47" Column="29" TopLine="31"/> </Position25> <Position26> <Filename Value="Physikunit.pas"/> - <Caret Line="1139" Column="6" TopLine="1106"/> + <Caret Line="283" Column="21" TopLine="212"/> </Position26> <Position27> <Filename Value="Physikunit.pas"/> - <Caret Line="88" Column="44" TopLine="68"/> + <Caret Line="334" TopLine="212"/> </Position27> <Position28> <Filename Value="Physikunit.pas"/> - <Caret Line="89" Column="44" TopLine="68"/> + <Caret Line="1452" Column="31" TopLine="1420"/> </Position28> <Position29> <Filename Value="Physikunit.pas"/> - <Caret Line="130" Column="55" TopLine="97"/> + <Caret TopLine="178"/> </Position29> <Position30> <Filename Value="Physikunit.pas"/> - <Caret Line="162" Column="18" TopLine="129"/> + <Caret Line="1049" Column="41" TopLine="748"/> </Position30> </JumpHistory> </ProjectSession> @@ -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 * λ) |