From 373ed9802b6d05e6c020553a68e6529d604cd5ef Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Wed, 29 Jul 2015 14:15:32 +0200 Subject: lauffaehig, bisher ohne Teilchen --- Physikunit.pas | 431 ++++++++++++++++++++++++++++++++------------------ Plasmapropagation | Bin 3387904 -> 3295712 bytes Plasmapropagation.lpi | 6 +- Plasmapropagation.lpr | 8 +- Plasmapropagation.lps | 171 ++++++++++---------- input.epost | 2 +- input.plap | 14 +- 7 files changed, 370 insertions(+), 262 deletions(-) diff --git a/Physikunit.pas b/Physikunit.pas index 492866c..1d50c5a 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -8,7 +8,7 @@ unit Physikunit; interface uses - Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, crt; + Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit; type tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier); @@ -24,27 +24,35 @@ type mfGamma,mfIGamma ); + tGitter = class; + { tAusgabeDatei } tAusgabeDatei = class - ableitung: boolean; - emInhalt: tEMFeldInhalt; - matInhalt: tMaterieFeldInhalt; - teilchen: longint; // wenn < 0, dann EM-Feld - name: string; - datei: file; - constructor create(feldName,prefix,suffix: string; prot: tProtokollant); + private + ableitung: boolean; + emInhalt: tEMFeldInhalt; + matInhalt: tMaterieFeldInhalt; + teilchen,nNum,tCnt: longint; // wenn teilchen < 0, dann EM-Feld + pre,suf: string; + datei: file; + pro: tProtokollant; + public + zeitAnz: longint; + sDT: extended; + constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); destructor destroy; override; + function dump: string; + procedure schreibeKopf; + procedure speichereWerte(gitter: tGitter; sDX: extended); end; - tGitter = class; - { tWertePunkt } tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen - matWerte: array of array[tMaterieFeldInhalt,boolean] of extended; + matWerte: array of array[tMaterieFeldInhalt,boolean] of double; // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen - emWerte: array[tEMFeldInhalt,boolean] of extended; + emWerte: array[tEMFeldInhalt,boolean] of double; // EM-Felder und deren Zeitableitungen // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! lN,rN: tWertePunkt; @@ -53,7 +61,7 @@ type destructor destroy; override; procedure berechneGammaUndP; procedure berechneNAbleitung(iDX,pDNMax: extended); - procedure berechneAbleitungen(dT,dX,iDX: extended); + procedure berechneAbleitungen(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; @@ -64,14 +72,14 @@ type private matAnz: longint; lichters: array[boolean] of tMyStringlist; - procedure setzeRaender(iDX: extended); + procedure setzeRaender(dT,iDX,iDT: extended); public inhalt: array of tWertePunkt; par: tGitter; constructor create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter); destructor destroy; override; - procedure berechneAbleitungen(dT,dX,iDX,pDNMax: extended); + procedure berechneAbleitungen(dT,dX,iDT,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; @@ -91,11 +99,11 @@ 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; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); + constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); destructor destroy; override; - procedure iteriereSchritt(dT: extended); - procedure macheAusgabe(aD: tAusgabedatei; sDX: extended); + procedure iteriereSchritt(dT,iDT: extended); function gibErhaltungsgroessen: string; + procedure berechne(was: char; teilchen: longint); end; { tSimulation } @@ -108,7 +116,7 @@ type fortschrittsAnzeige: boolean; ausgabeDateien: array of tAusgabeDatei; public - constructor create(inName: string; Protokollant: tProtokollant); + constructor create(inName: string; protokollant: tProtokollant; name: string); destructor destroy; override; function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? end; @@ -125,7 +133,7 @@ const { tAusgabeDatei } -constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant); +constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); var emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; @@ -133,6 +141,7 @@ var abl,gef: boolean; begin inherited create; + pro:=tProtokollant.create(prot,nam); num:=length(feldName); while (num>0) and (feldName[num] in ['0'..'9']) do dec(num); @@ -147,6 +156,9 @@ begin emInhalt:=efA; matInhalt:=mfN; feldName:=ansiUpperCase(feldName); + nNum:=0; // Header 0 wurde also noch nicht geschrieben + zeitAnz:=-1; + tCnt:=-1; gef:=false; for abl:=false to true do begin @@ -168,29 +180,108 @@ begin end; if not gef then begin - prot.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true); + pro.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true); + pro.destroyAll; halt(1); end; if teilchen>=0 then - name:=matFeldNamen[maF]+inttostr(teilchen) + pre:=matFeldNamen[maF]+inttostr(teilchen) else - name:=emFeldNamen[emF]; + pre:=emFeldNamen[emF]; - if abl then - name:='d'+name; - name:=prefix+name+suffix; + if ableitung then + pre:='d'+pre; - assignFile(datei,name); - rewrite(datei); + pre:=prefix+pre; + suf:=suffix; end; destructor tAusgabeDatei.destroy; begin - closeFile(datei); + 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); + halt(1); + end; inherited destroy; end; +function tAusgabeDatei.dump: string; +begin + if teilchen<0 then + result:=emFeldNamen[emInhalt] + else + result:=matFeldNamen[matInhalt]+inttostr(teilchen); + if ableitung then + result:=result+''''; + result:=result+' -> '+pre+'-$i'+suf; +end; + +procedure tAusgabeDatei.schreibeKopf; +begin + 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); + halt(1); + end; + tCnt:=0; + assignfile(datei,pre+'-'+inttostr(nNum)+suf); + inc(nNum); + rewrite(datei,1); + blockwrite(datei,nNum,sizeof(integer)); + blockwrite(datei,zeitAnz,sizeof(integer)); + closefile(datei); +end; + +procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sDX: extended); +var i,cnt: longint; + sX,cX: extended; +begin + if (teilchen<0) and (emInhalt=efA) then + gitter.berechne('A',teilchen); + if (teilchen>=0) and (matInhalt=mfP) then + gitter.berechne('P',teilchen); + + if gitter.t>=nNum+sDT/2 then + schreibeKopf; + + cnt:=floor((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)/sDX*gitter.dX+1); + cX:=gitter.xl; + 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)); + 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(double)); + sX:=sX+sDX; + end; + cX:=cX+gitter.dX; + end; + end + else 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(double)); + 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); + halt(1); + end; +end; + { tWertePunkt } constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint); @@ -275,7 +366,7 @@ begin // alpha = pDNMax * dX end; -procedure tWertePunkt.berechneAbleitungen(dT,dX,iDX: extended); // Zeitableitungen berechnen +procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen var i: longint; begin @@ -316,7 +407,9 @@ begin end; // dA/dt = dA/dt - emWerte[efAX,true]:=emWerte[efDAXDT,false] + emWerte[efDAXDT,true]*dT; + emWerte[efAX,true]:=emWerte[efDAXDT,false];// + emWerte[efDAXDT,true]*dT; + emWerte[efAY,true]:=emWerte[efDAYDT,false];// + emWerte[efDAYDT,true]*dT; + emWerte[efAZ,true]:=emWerte[efDAZDT,false];// + emWerte[efDAZDT,true]*dT; end; procedure tWertePunkt.liKo(in1,in2: tWertePunkt; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt @@ -399,7 +492,7 @@ begin end; parent.kvs.rem('x'); for rechts:=false to true do begin - lichters[rechts]:=tMyStringlist.create(nil); + lichters[rechts]:=tMyStringlist.create(nil,''); lichters[rechts].text:=lichter.text; end; lichters[false].grep('^links .*'); @@ -415,9 +508,12 @@ begin inherited destroy; end; -procedure tFelder.setzeRaender(iDX: extended); +procedure tFelder.setzeRaender(dT,iDX,iDT: extended); var - emF: tEMFeldInhalt; + emF: tEMFeldInhalt; + rechts: boolean; + i: longint; + oVal: extended; begin for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld inhalt[0].emWerte[emF,true]:= @@ -428,15 +524,18 @@ begin 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; *) + for rechts:=false to true do + for i:=0 to lichters[rechts].count-1 do begin + par.kvs.add('t',par.t-dT); + oVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil); + par.kvs.add('t',par.t); + inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true]:= + inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true] + + (exprToFloat(false,lichters[rechts][i],par.kvs,nil) - oVal)*iDT; + end; end; -procedure tFelder.berechneAbleitungen(dt,dx,iDX,pDNMax: extended); +procedure tFelder.berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended); var i: longint; begin @@ -449,12 +548,12 @@ begin -abs(inhalt[length(inhalt)-2].matWerte[i,mfPX,false]); end; - setzeRaender(iDX); + setzeRaender(dT,iDX,iDT); 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); + inhalt[i].berechneAbleitungen(dX,iDX); end; procedure tFelder.liKo(in1,in2: tFelder; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt @@ -475,14 +574,14 @@ end; { tGitter } -constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); +constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); var i: longint; begin inherited create; zeitverfahren:=ZV; kvs:=bekannteWerte; - Prot:=TProtokollant.create(Protokollant,'tGitter'); + prot:=tProtokollant.create(protokollant,name); dTMaximum:=deltaT; dX:=deltaX; iDX:=1/dX; @@ -502,6 +601,7 @@ begin aktuelleFelder:=0; t:=0; + kvs.add('t',t); end; destructor tGitter.destroy; @@ -514,28 +614,30 @@ begin inherited destroy; end; -procedure tGitter.iteriereSchritt(dT: extended); +procedure tGitter.iteriereSchritt(dT,iDT: extended); var i: longint; {$IFDEF Zeitschrittueberwachung} pro: tProtokollant; j: longint; {$ENDIF} begin - felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDX,pDNMax); // y' = y'(t,y(t)) + t:=t+dT; + kvs.add('t',t); + felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t)) case zeitverfahren of zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT); zfRungeKuttaVier: begin - 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[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2 + felders[2].berechneAbleitungen(dT/2,dX,iDT,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[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2 + felders[3].berechneAbleitungen(dT/2,dX,iDT,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[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt + felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc) felders[1-aktuelleFelder].liKo( felders[aktuelleFelder], @@ -551,7 +653,6 @@ begin end; end{of case}; - t:=t+dT; aktuelleFelder:=1-aktuelleFelder; {$IFDEF Dichteueberwachung} with felders[aktuelleFelder] do @@ -581,58 +682,10 @@ begin {$ENDIF} end; -procedure tGitter.macheAusgabe(aD: tAusgabedatei; sDX: extended); -var i: longint; - sX,cX,tmp: extended; - emF: tEMFeldInhalt; - maF: tMaterieFeldInhalt; -begin - if (aD.teilchen<0) and (aD.emInhalt=efA) then - for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin - tmp:=0; - 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; - 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; - 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: integer; ns: tExtendedArray; - Pro: TProtokollant; + pro: tProtokollant; begin pro:=tProtokollant.create(prot,'gibErhaltungsgroessen'); setlength(ns,felders[aktuelleFelder].matAnz); @@ -684,23 +737,53 @@ begin end; end; +procedure tGitter.berechne(was: char; teilchen: longint); +var + i: longint; + tmp: extended; + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; +begin + if was='A' then begin + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + tmp:=0; + 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; + exit; + end; + if was='P' then begin + 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[teilchen,maF,false]); + felders[aktuelleFelder].inhalt[i].matWerte[teilchen,maF,false]:=sqrt(tmp); + end; + exit; + end; + prot.schreibe('Kann '''+was+''' nicht berechnen, weil ich es nicht verstehe!',true); + halt; +end; + { tSimulation } -constructor tSimulation.create(inName: string; Protokollant: tProtokollant); +constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string); var ifile: tMyStringlist; zeitverfahren: tZeitverfahren; s,t,aSuffix,aPrefix: string; - deltaX,deltaT,breite,pDNMax: extended; - i,j: longint; + deltaX,breite,pDNMax: extended; + i: longint; kvs: tKnownValues; dichten,lichter: tMyStringlist; + pro: tProtokollant; begin inherited create; - prot:=tProtokollant.create(Protokollant,'tSimulation'); + prot:=tProtokollant.create(protokollant,name); kvs:=tKnownValues.create; - ifile:=tMyStringlist.create(prot); + ifile:=tMyStringlist.create(prot,'ifile'); ifile.loadfromfile(inName); if not ifile.unfoldMacros then begin prot.schreibe('Fehlerhafte Macros in Parameterdatei '''+inName+'''!',true); @@ -712,8 +795,9 @@ begin deltaX:=1e-2; kvs.add('λ',1); kvs.add('T',1); + kvs.add('ω',2*pi); kvs.add('dX',deltaX); - deltaT:=-1; + dT:=-1; sDT:=-1; sDX:=-1; tEnde:=100; @@ -727,18 +811,26 @@ begin setlength(ausgabeDateien,0); // Standardeinstellungen Breich 'lichtVon...' - lichter:=tMyStringlist.create(prot); + lichter:=tMyStringlist.create(prot,'lichter'); // Standardeinstellungen Breich 'teilchen...' - dichten:=tMyStringlist.create(prot); + dichten:=tMyStringlist.create(prot,'teilchen'); + + repeat + if not ifile.readln(s) then begin + prot.schreibe('Unerwartetes Dateiende in '''+inName+'''!',true); + halt(1); + end; + + if s='Dateiende' then + break; - 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; writeln(s); + end; if s='allgemeinEnde' then break; if s='runge-Kutta-4' then begin Zeitverfahren:=zfRungeKuttaVier; @@ -754,8 +846,8 @@ begin continue; end; if startetMit('zeitschritt',s) then begin - deltaT:=exprtofloat(false,s,kvs,nil); - kvs.add('dT',deltaT); + dT:=exprtofloat(false,s,kvs,nil); + kvs.add('dT',dT); continue; end; if startetMit('diffusionsterm',s) then begin @@ -812,7 +904,7 @@ begin while s<>'' do begin t:=erstesArgument(s,','); setlength(ausgabeDateien,length(ausgabeDateien)+1); - ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot); + ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot,'ausgabeDateien['+inttostr(length(ausgabeDateien)-1)+']'); end; continue; end; @@ -822,32 +914,73 @@ begin continue; end; + if startetMit('licht von ',s) then begin + if startetMit('links ',s) then begin + lichter.add('links '+s); + continue; + end; + if startetMit('rechts ',s) then begin + lichter.add('rechts '+s); + continue; + end; + prot.schreibe('Licht kann nicht von '''+erstesArgument(s)+''' kommen!',true); + halt(1); + end; + prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true); halt(1); - end; + until false; + + ifile.free; if length(ausgabedateien)=0 then begin prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true); halt(1); end; - if deltaT<0 then - deltaT:=deltaX/10; + if dT<0 then + dT:=deltaX/10; if sDT<0 then - sDT:=deltaT; + sDT:=dT; + sDT:=1/round(1/sDT); + sT:=sDT/2; 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(tEnde/sDT+1); - BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); + pro:=tProtokollant.create(prot,'create'); + case zeitverfahren of + zfEulerVorwaerts: + pro.schreibe('Iteration mittels Euler-Vorwärts'); + zfRungeKuttaVier: + pro.schreibe('Iteration mittels Runge-Kutta-4'); + else + pro.schreibe('Iteration mittels unbekanntem Verfahren'); + end{of case}; + pro.schreibe('dX = '+floattostr(deltaX)); + pro.schreibe('dT = '+floattostr(dT)); + pro.schreibe('breite = '+floattostr(breite)+' (= '+inttostr(round(breite/deltaX)+1)+' Punkte)'); + pro.schreibe('pDNMax = '+floattostr(pDNMax)); + pro.schreibe('bekannte Werte:'); + kvs.dump(pro,' '); + if dichten.count>0 then begin + pro.schreibe('dichten:'); + dichten.dump(pro,' '); end; - ifile.free; + if lichter.count>0 then begin + pro.schreibe('lichter:'); + lichter.dump(pro,' '); + end; + if length(ausgabeDateien)>0 then begin + pro.schreibe('ausgaben:'); + for i:=0 to length(ausgabeDateien)-1 do begin + ausgabeDateien[i].zeitAnz:=round(1/sDT); + ausgabeDateien[i].sDT:=sDT; + pro.schreibe(' '+ausgabeDateien[i].dump); + end; + end; + pro.free; - gitter:=tGitter.create(round(breite/deltaX),deltaT,deltaX,pDNMax,kvs,dichten,lichter,zeitverfahren,prot); + gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,dichten,lichter,zeitverfahren,prot,'gitter'); dichten.free; lichter.free; end; @@ -862,51 +995,33 @@ 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); + gitter.iteriereSchritt(dT,1/dT); zeitPhysik:=zeitPhysik+now; zeitDatei:=zeitDatei-now; - while gitter.t>=sT do begin + while (gitter.t>=sT) and (sT=tEnde; + result:=gitter.t - @@ -78,10 +77,7 @@ - - - - + diff --git a/Plasmapropagation.lpr b/Plasmapropagation.lpr index 5c914d5..5559f51 100644 --- a/Plasmapropagation.lpr +++ b/Plasmapropagation.lpr @@ -8,9 +8,7 @@ uses {$ENDIF}{$ENDIF} Classes, SysUtils, CustApp, { you can add units after this } - math, Physikunit, crt, protokollunit; - -var ErrorCode: longint; + math, Physikunit, protokollunit; type @@ -45,7 +43,7 @@ begin zeitDatei:=0; zeitPhysik:=0; - simulation:=tSimulation.create(paramstr(1),prot); + simulation:=tSimulation.create(paramstr(1),prot,'simulation'); while simulation.iteriereSchritt(start,zeitPhysik,zeitDatei) do ; @@ -66,7 +64,7 @@ begin prot.schreibe('Das hat '+s+' gedauert,',true); prot.schreibe(' davon '+t+' für Dateizugriffe',true); prot.schreibe('und nur '+u+' für die eigentliche Physik!',true); - prot.Free; + prot.free; // stop program loop Terminate; diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps index cd5e8e9..cbe3abd 100644 --- a/Plasmapropagation.lps +++ b/Plasmapropagation.lps @@ -7,74 +7,72 @@ - - - + + - - - + + + + - + - - - - + + + + + - - - - + + + - - - - + + + + + - - - - + + + - - - - - + + + + - - - - - - + + + + + @@ -107,121 +105,122 @@ + - + - - + + - - + + - - + + - - + + - - + + - - + + - + + - - + + - - + + - - + + - - + + - - + + - - + + - + - + - - + + - - + + - - + - - + + - - + + - - + + - + - + - - + + - - + + - - + + - - + + - - + + diff --git a/input.epost b/input.epost index 8f9dc9e..68a9f8f 100644 --- a/input.epost +++ b/input.epost @@ -18,7 +18,7 @@ Ende Threadanzahl: 10 # !Schleife: $Feld: dDPHIDX dDAYDT dAY dAX -!Schleife: $Feld: AX +!Schleife: $Feld: AY DAYDT Daten einlesen Genauigkeit: extended diff --git a/input.plap b/input.plap index d2debd1..673a4e2 100644 --- a/input.plap +++ b/input.plap @@ -3,9 +3,9 @@ allgemein runge-Kutta-4 ortsschritt 10^-2 * λ - zeitschritt 10^-3 * T - zeit 100 * T - breite 10 * λ + zeitschritt 10^-2 * T + zeit 20 * T + breite 5 * λ mit Fortschrittsanzeige allgemeinEnde @@ -15,12 +15,12 @@ ausgaben # felder AX,AY,dAYDT,N,dPhiDX,dPsiDX,Gamma ausgabenEnde -!setze $tFwhm: 5 T -!setze $tMitte: 20 T +!setze $tFwhm: (2.5 * T) +!setze $tMitte: (1 * T) -lichtVonLinks 0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * Sin(ω*t) +licht von links 0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * sin(ω*t) -dateiEnde +Dateiende teilchen1 ladung -q -- cgit v1.2.3-70-g09d2