unit Physikunit; {$mode objfpc}{$H+} {$IFNDEF UNIX} {$ERROR This program can be compiled only on/for Unix/Linux based systems.} {$ENDIF} {$DEFINE Zeitschrittueberwachung} {$DEFINE Dichteueberwachung} { $DEFINE negativeDichteueberwachung} interface uses Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix; type tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn); tVerteilungsfunktion = function(x: extended): extended; tEMFeldInhalt = ( efA,efAX,efAY,efAZ, efDAXDT,efDAYDT,efDAZDT, efDPhiDX ); tMaterieFeldInhalt = ( mfN,mfDPsiDX, // Teilchendichte, Geschwindigkeitsdichte mfP,mfPX,mfPY,mfPZ, // Impuls mfGamma,mfIGamma // rel. Gamma-Faktor ); const erstesEMFMitAbleitung = efAX; letztesEMFMitAbleitung = efDAZDT; erstesMatFMitAbleitung = mfN; letztesMatFMitAbleitung = mfDPsiDX; type tGitter = class; tFelder = class; tSimulation = class; { tAusgabeDatei } tAusgabeDatei = class private ableitung: boolean; emInhalt: tEMFeldInhalt; matInhalt: tMaterieFeldInhalt; teilchen,nNum,tCnt: longint; // wenn teilchen < 0, dann EM-Feld pre,suf: string; datei: file; pro: tProtokollant; buf: pointer; bufLen,bufPos: longint; procedure schreibeBuffer(minBufLen: longint); 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; sT, sDX: extended); end; { tTeilchenSpezies } tTeilchenSpezies = class private ortsGrenzen: array of extended; dichteStuecke: array of string; public nMin,nMax, // minimale und maximale Massendichte in nc spezifischeLadung: extended; // q/m in e/me constructor create; destructor destroy; override; procedure dump(prot: tProtokollant; prefix: string); function gibDichte(x: extended; kvs: tKnownValues): extended; // Massendichte procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant); end; { tWertePunkt } tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen private besitzer: tFelder; public 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, dPhiDX, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! lN,rN: tWertePunkt; constructor create(linkerNachbar: tWertePunkt; derBesitzer: tFelder); destructor destroy; override; procedure berechneGammaUndP; procedure berechneDPhiDX(dX: extended); procedure berechneNAbleitung(iDX,pDNMax: extended); overload; procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload; procedure berechneAbleitungen(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; procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload; procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6: tWertePunkt; fak2,fak3,fak4,fak5,fak6: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; function maxDT: extended; function dump(mitNachbarn: boolean): string; procedure nichtnegativieren; end; { TFelder } tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private matAnz: longint; spezLadungen, massen: tExtendedArray; lichters: array[boolean] of tMyStringlist; procedure setzeRaender(iDX: extended); public inhalt: array of tWertePunkt; par: tGitter; constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); destructor destroy; override; procedure berechneAbleitungen(dX,iDX,pDNMax: extended); procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt procedure liKo(in1,in2,in3: tFelder; fak2,fak3: extended); overload; procedure liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); overload; procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6: tFelder; fak2,fak3,fak4,fak5,fak6: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7: tFelder; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; function maxDT: extended; end; { tGitter } tGitter = class(tObject) private besitzer: tSimulation; prot: tProtokollant; zeitverfahren: tZeitverfahren; pDNMaxDynamisch,abbruch: boolean; dX,iDX,pDNMax,dTMaximum,xl,t: extended; kvs: tKnownValues; procedure diffusionsTermAnpassen(pro: tProtokollant); function pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert? procedure abbrechen; public aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern constructor create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); destructor destroy; override; procedure iteriereSchritt(var dT: extended; dTMin: extended); function dumpErhaltungsgroessen: boolean; procedure berechne(was: char; teilchen: longint); end; { tSimulation } tSimulation = class(tObject) private prot: tProtokollant; gitter: tGitter; dT,tEnde,sT,sDT,sDX,dTMin: extended; fortschrittsAnzeige: boolean; ausgabeDateien: array of tAusgabeDatei; public gotSighup: boolean; 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; procedure SignalCapture(signal : longint); cdecl; implementation const 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' ); minAusgabeBuffer = 1024*1024; var sighupSimulationen: array of tSimulation; { tAusgabeDatei } constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); var emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; num: longint; 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); inc(num); if num<=length(feldName) then begin teilchen:=strtoint(copy(feldName,num,length(feldName)))-1; delete(feldName,num,length(feldName)); end else teilchen:=0; 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 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 pro.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true); pro.destroyAll; halt(1); end; if teilchen>=0 then pre:=matFeldNamen[maF]+inttostr(teilchen+1) else pre:=emFeldNamen[emF]; 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; halt(1); end; 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 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); pro.destroyall; halt(1); end; if bufPos<>0 then schreibeBuffer(0); 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; sT, sDX: extended); var i,cnt: longint; sX,cX: extended; begin if (teilchen>=gitter.felders[gitter.aktuelleFelder].matAnz) then begin pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true); pro.destroyall; halt(1); end; if (teilchen<0) and (emInhalt=efA) then gitter.berechne('A',teilchen); if (teilchen>=0) and (matInhalt=mfP) then gitter.berechne('P',teilchen); if sT>=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); 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); 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; end; end else begin for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin while cX>=sX do begin dec(cnt); 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; if cnt<>0 then begin pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true); pro.destroyall; halt(1); end; end; { tTeilchenSpezies } constructor tTeilchenSpezies.create; begin inherited create; fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0); setlength(ortsGrenzen,0); fillchar(dichteStuecke,sizeof(dichteStuecke),#0); setlength(dichteStuecke,0); nMin:=1E-3; nMax:=nMin; spezifischeLadung:=0; end; destructor tTeilchenSpezies.destroy; begin setlength(ortsGrenzen,0); setlength(dichteStuecke,0); inherited destroy; end; procedure tTeilchenSpezies.dump(prot: tProtokollant; prefix: string); var i: longint; begin prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc'); prot.schreibe(prefix+'q/m = '+floattostr(-spezifischeLadung)+' e/me'); prot.schreibe(prefix+'n(x):'); for i:=0 to length(dichteStuecke)-1 do begin prot.schreibe(prefix+' '+dichteStuecke[i]); if i verteilung stückweise!',true); prot.destroyall; halt(1); end; if length(dichteStuecke)=length(ortsGrenzen) then begin // neues Funktionsstück wird definiert setlength(dichteStuecke,length(dichteStuecke)+1); dichteStuecke[length(dichteStuecke)-1]:=s; continue; end; // kein neues Funktionsstück => if s='stückweiseEnde' then break; // Ende oder // neue Ortsgrenze setlength(ortsGrenzen,length(ortsGrenzen)+1); ortsGrenzen[length(ortsGrenzen)-1]:=exprToFloat(false,s,kvs,nil); until false; end else if startetMit('wie teilchen',rest) then begin ori:=strtoint(rest)-1; if (ori<0) or (ori>=length(teilchen)-1) then begin prot.schreibe('Kann mich nicht auf die Dichte von Teilchen '+inttostr(ori+1)+' beziehen, weil es noch nicht definiert wurde!',true); prot.destroyall; halt(1); end; setlength(ortsGrenzen,length(teilchen[ori].ortsGrenzen)); for i:=0 to length(ortsGrenzen)-1 do ortsGrenzen[i]:=teilchen[ori].ortsGrenzen[i]; setlength(dichteStuecke,length(teilchen[ori].dichteStuecke)); for i:=0 to length(dichteStuecke)-1 do dichteStuecke[i]:=teilchen[ori].dichteStuecke[i]; end else begin setlength(ortsGrenzen,0); setlength(dichteStuecke,1); dichteStuecke[0]:=rest; end; end; { tWertePunkt } constructor tWertePunkt.create(linkerNachbar: tWertePunkt; derBesitzer: tFelder); var emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; abl: boolean; i: longint; begin inherited create; besitzer:=derBesitzer; fillchar(matWerte,sizeof(matWerte),#0); setlength(matWerte,length(besitzer.spezLadungen)); for i:=0 to length(matWerte)-1 do for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do for abl:=false to true do matWerte[i,maF,abl]:=0; for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do for abl:=false to true do emWerte[emF,abl]:=0; lN:=linkerNachbar; rN:=nil; 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 tWertePunkt.berechneGammaUndP; // aus aktuellem dPsi/dX und A var i: longint; begin for i:=0 to length(matWerte)-1 do begin matWerte[i,mfPX,false]:= besitzer.spezLadungen[i] * emWerte[efAX,false] + matWerte[i,mfDPsiDX,false]; matWerte[i,mfPY,false]:= besitzer.spezLadungen[i] * emWerte[efAY,false]; matWerte[i,mfPZ,false]:= besitzer.spezLadungen[i] * 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; end; procedure tWertePunkt.berechneDPhiDX(dX: extended); var i: longint; begin // d2Phi/dx2 = rho ... hier muss integriert werden: // dPhi/dx = dPhi/dx(x- deltax) + * deltax if assigned(lN) then begin emWerte[efDPhiDX,false]:=lN.emWerte[efDPhiDX,false]; for i:=0 to length(matWerte)-1 do emWerte[efDPhiDX,false]:= emWerte[efDPhiDX,false] + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,false] + matWerte[i,mfN,false]) * dX/2; end else begin emWerte[efDPhiDX,false]:=0; for i:=0 to length(matWerte)-1 do emWerte[efDPhiDX,false]:= emWerte[efDPhiDX,false] + besitzer.spezLadungen[i] * matWerte[i,mfN,false] * dX/2; end; end; 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]:= -( rN^.rN^.Groessen[false,fiN]*rN^.rN^.Groessen[false,fiiGamma]*rN^.rN^.Groessen[false,fiPX] + rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*rN^.Groessen[false,fiPX] - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*lN^.Groessen[false,fiPX] - lN^.lN^.Groessen[false,fiN]*lN^.lN^.Groessen[false,fiiGamma]*lN^.lN^.Groessen[false,fiPX] )*iDX/6; // es wird über die beiden nächsten linken bzw. rechten Nachbarn gemittelt *) 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] * 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^2 end; procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); var i: longint; begin if rechts then 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] * 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] * iDX) - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax + matWerte[i,mfPX,false] * iDX)) / 2; end end; procedure tWertePunkt.berechneAbleitungen(iDX: extended); // Zeitableitungen berechnen var i: longint; begin // d2Psi/dxdt = dPhi/dx - dGamma/dx for i:=0 to length(matWerte)-1 do matWerte[i,mfDPsiDX,true]:= besitzer.spezLadungen[i] * emWerte[efDPhiDX,false] - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; // d2A/dt2 = Laplace(A) ... ( - Nabla(dPhi/dt) wird auf dA/dt aufgeschlagen !!! ) emWerte[efDAXDT,true]:= ( rN.emWerte[efAX,false] - 2*emWerte[efAX,false] + lN.emWerte[efAX,false] )*sqr(iDX); 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); // ... - Rho/gamma p for i:=0 to length(matWerte)-1 do begin emWerte[efDAXDT,true]:= emWerte[efDAXDT,true] - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]); emWerte[efDAYDT,true]:= emWerte[efDAYDT,true] - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]); emWerte[efDAZDT,true]:= emWerte[efDAZDT,true] - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]); end; // dA/dt = dA/dt - Nabla(Phi) emWerte[efAX,true]:=emWerte[efDAXDT,false] - emWerte[efDPhiDX,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; function tWertePunkt.maxDT: extended; var i: longint; begin result:=infinity; for i:=0 to length(matWerte)-1 do if (matWerte[i,mfN,true]<0) and (matWerte[i,mfN,false]>0) then result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]); end; function tWertePunkt.dump(mitNachbarn: boolean): string; var i: longint; begin if mitNachbarn and assigned(lN) then result:=lN.dump(false)+' << ' else result:=''; for i:=0 to length(matWerte)-1 do result:= result+ 'N['+ inttostr(i)+'] = '+floattostr(matWerte[i,mfN,false])+' '+ 'N''['+inttostr(i)+'] = '+floattostr(matWerte[i,mfN,true])+' '; result:=result+'maxDT: '+floattostr(maxDT); if mitNachbarn and assigned(rN) then result:=result+' >> '+rN.dump(false); end; procedure tWertePunkt.nichtnegativieren; // Dichten nicht negativ machen var i: longint; defizit: extended; li,re: tWertePunkt; pro: tProtokollant; begin for i:=0 to length(matWerte)-1 do if matWerte[i,mfN,false]<0 then begin defizit:=-matWerte[i,mfN,false]; matWerte[i,mfN,false]:=0; li:=self.lN; re:=self.rN; repeat if assigned(li) then begin if li.matWerte[i,mfN,false]>0 then begin if defizit>=li.matWerte[i,mfN,false] then begin defizit:=defizit-li.matWerte[i,mfN,false]; li.matWerte[i,mfN,false]:=0; end else begin li.matWerte[i,mfN,false]:=li.matWerte[i,mfN,false]-defizit; defizit:=0;; end; end; li:=li.lN; end; if assigned(re) then begin if re.matWerte[i,mfN,false]>0 then begin if defizit>=re.matWerte[i,mfN,false] then begin defizit:=defizit-re.matWerte[i,mfN,false]; re.matWerte[i,mfN,false]:=0; end else begin re.matWerte[i,mfN,false]:=re.matWerte[i,mfN,false]-defizit; defizit:=0;; end; end; re:=re.rN; end; until (defizit<=0) or not (assigned(li) or assigned(re)); if defizit>0 then begin pro:=tProtokollant.create(besitzer.par.prot,'WertePunkt.nichtnegativieren'); pro.schreibe('Ich konnte ein Teilchendefizit nicht auffüllen!',true); besitzer.par.abbrechen; end; end; end; { linearkombinationsspezifische Methoden von tWertePunkt und tFelder } {$INCLUDE linearkombination.inc} {$DEFINE lkA3} {$INCLUDE linearkombination.inc} {$DEFINE lkA4} {$INCLUDE linearkombination.inc} {$DEFINE lkA5} {$INCLUDE linearkombination.inc} {$DEFINE lkA6} {$INCLUDE linearkombination.inc} {$DEFINE lkA7} {$INCLUDE linearkombination.inc} {$DEFINE lkA8} {$INCLUDE linearkombination.inc} {$DEFINE lkA9} {$INCLUDE linearkombination.inc} {$DEFINE lkA10} {$INCLUDE linearkombination.inc} {$DEFINE lkA11} {$INCLUDE linearkombination.inc} {$DEFINE lkA12} {$INCLUDE linearkombination.inc} {$DEFINE lkA14} {$INCLUDE linearkombination.inc} {$DEFINE lkA15} {$INCLUDE linearkombination.inc} {$DEFINE lkA16} {$INCLUDE linearkombination.inc} {$DEFINE lkA17} {$INCLUDE linearkombination.inc} {$DEFINE lkA18} {$INCLUDE linearkombination.inc} {$DEFINE lkA19} {$INCLUDE linearkombination.inc} {$DEFINE lkA22} {$INCLUDE linearkombination.inc} {$DEFINE lkA23} {$INCLUDE linearkombination.inc} {$DEFINE lkA31} {$INCLUDE linearkombination.inc} {$DEFINE lkA32} {$INCLUDE linearkombination.inc} { tFelder } constructor tFelder.create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); var i,j: longint; rechts: boolean; begin inherited create; par:=parent; matAnz:=length(teilchen); fillchar(spezLadungen,sizeof(spezLadungen),#0); setlength(spezLadungen,matAnz); for i:=0 to length(spezLadungen)-1 do spezLadungen[i]:=teilchen[i].spezifischeLadung; fillchar(inhalt,sizeof(inhalt),#0); setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen inhalt[0]:=tWertePunkt.create(nil,self); for i:=1 to length(inhalt)-1 do inhalt[i]:=tWertePunkt.create(inhalt[i-1],self); fillchar(massen,sizeof(massen),#0); setlength(massen,length(teilchen)); for i:=0 to length(massen)-1 do massen[i]:=0; 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,mfDPsiDX,false]:=0; massen[j]:=massen[j]+inhalt[i].matWerte[j,mfN,false]*par.dX; end; for rechts:=false to true do begin lichters[rechts]:=tMyStringlist.create(nil,''); lichters[rechts].text:=lichter.text; end; lichters[false].grep('^links .*'); lichters[false].subst('^links *',''); lichters[true].grep('^rechts .*'); lichters[true].subst('^rechts *',''); end; destructor tFelder.destroy; var i: longint; begin setlength(spezLadungen,0); for i:=0 to length(inhalt)-1 do inhalt[i].free; setlength(inhalt,0); lichters[false].free; lichters[true].free; inherited destroy; end; procedure tFelder.setzeRaender(iDX: extended); var emF: tEMFeldInhalt; rechts: boolean; i: longint; 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) for rechts:=false to true do for i:=0 to lichters[rechts].count-1 do 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); end; procedure tFelder.berechneAbleitungen(dX,iDX,pDNMax: extended); var i: longint; begin for i:=0 to length(inhalt)-1 do inhalt[i].berechneGammaUndP; for i:=0 to length(inhalt)-1 do inhalt[i].berechneDPhiDX(dX); 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; for i:=2 to length(inhalt)-3 do inhalt[i].berechneNAbleitung(iDX,pDNMax); inhalt[1].berechneNAbleitung(iDX,pDNMax,false); inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true); setzeRaender(iDX); for i:=1 to length(inhalt)-2 do inhalt[i].berechneAbleitungen(iDX); end; function tFelder.maxDT: extended; var i: longint; begin result:=inhalt[0].maxDT; for i:=1 to length(inhalt)-1 do result:=min(result,inhalt[i].maxDT); end; { tGitter } constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); var i: longint; begin inherited create; abbruch:=false; besitzer:=derBesitzer; zeitverfahren:=zv; kvs:=bekannteWerte; prot:=tProtokollant.create(protokollant,name); dTMaximum:=deltaT; dX:=deltaX; iDX:=1/dX; if pDNMa < 0 then begin pDNMax:=0; pDNMaxDynamisch:=true; end else pDNMax:=pDNMa; case Zeitverfahren of zfEulerVorwaerts: Setlength(Felders,2); zfRungeKuttaDreiAchtel,zfRungeKuttaVier: Setlength(Felders,5); zfRungeKuttaZehn: Setlength(Felders,18); zfRungeKuttaZwoelf: Setlength(Felders,26); zfRungeKuttaVierzehn: Setlength(Felders,36); end{of Case}; xl:=dX/2; for i:=0 to length(felders)-1 do felders[i]:=tFelder.create(size,teilchen,lichter,self); aktuelleFelder:=0; t:=0; kvs.add('t',t); end; destructor tGitter.destroy; var i: longint; begin for i:=0 to length(Felders)-1 do felders[i].free; setlength(felders,0); inherited destroy; end; procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant); var i,j,k: longint; curMax,tmp: extended; begin curMax:=0; k:=0; for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do for j:=0 to felders[aktuelleFelder].matAnz-1 do begin tmp:=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-100)); if tmp>curMax then begin curMax:=tmp; k:=i; end; end; if curMax > pDNMax then begin if pDNMaxDynamisch then begin pDNMax:=2*curMax; 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, x='+floattostr(xL+k*dX)+'), die numerische Stabilität ist nicht mehr gewährleistet, Simulation wird abgebrochen!',true); pro.schreibe(' Lösung: größeren Diffusionsterm wählen',true); pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern',true); abbrechen; end; end; end; function tGitter.pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert? begin dTMax:=min(dTMax,felders[aF].maxDT); // das maximale dT, welches bei den momentanen Ableitungen nicht zu // unphysikalischen Effekten (z.B. negativen Dichten) führt result:=(dTMax1.1*dTMin); if result then begin // beabsichtigter Zeitschritt ist zu groß, dT:=dTMax/4; // dann machen wir ihn doch kleiner! if dT '+floattostr(dT)+' (t = '+floattostr(t)+')'); prot.spuelen; end {$ENDIF}; end; end; procedure tGitter.abbrechen; var i,j,k: longint; maF: tMaterieFeldInhalt; emF: tEMFeldInhalt; abl: boolean; begin abbruch:=true; for i:=0 to length(felders)-1 do for j:=0 to length(felders[i].inhalt)-1 do begin for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do for abl:=false to true do felders[i].inhalt[j].emWerte[emF,false]:=0; for k:=0 to felders[i].matAnz-1 do for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do for abl:=false to true do felders[i].inhalt[j].matWerte[k,maF,false]:=0; end; end; procedure tGitter.iteriereSchritt(var dT: extended; dTMin: extended); var {$IFDEF negativeDichteueberwachung} i,j: longint; pro: tProtokollant; {$ENDIF} dTMax: extended; begin if abbruch then begin t:=t+dT; exit; end; kvs.add('t',t); diffusionsTermAnpassen(prot); dTMax:=dX*100; repeat felders[aktuelleFelder].berechneAbleitungen(dX,iDX,pDNMax); // y' = y'(t,y(t)) if pruefeMaxDT(aktuelleFelder,dTMax,dT,dTMin) then continue; case zeitverfahren of zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT); zfRungeKuttaDreiAchtel: begin {$INCLUDE rk3_8.inc} end; zfRungeKuttaVier: begin {$INCLUDE rk4.inc} end; zfRungeKuttaZehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk108.txt {$INCLUDE rk108.inc} end; zfRungeKuttaZwoelf: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1210.txt {$INCLUDE rk1210.inc} end; zfRungeKuttaVierzehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1412.txt {$INCLUDE rk1412.inc} end; end{of case}; break; until abbruch; aktuelleFelder:=1-aktuelleFelder; {$IFDEF negativeDichteueberwachung} for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do for j:=0 to felders[aktuelleFelder].matAnz-1 do if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then begin pro:=tProtokollant.create(prot,'iteriereSchritt'); pro.schreibe('n<0 bei:',true); pro.schreibe(' t = ' +floattostr(t)+',',true); pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); pro.schreibe(' alte Werte:',true); with felders[1-aktuelleFelder] do begin pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true); pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]* inhalt[i].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true); pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true); pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]* inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); end; pro.schreibe(' neue Werte:',true); with felders[aktuelleFelder] do begin pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true); pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]* inhalt[i].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true); pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true); pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]* inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); end; pro.free; abbrechen; end; {$ENDIF} t:=t+dT; dT:=max(dT,dTMax/100); end; function tGitter.dumpErhaltungsgroessen: boolean; var i,j: integer; ns: tExtendedArray; pro: tProtokollant; s: string; begin pro:=tProtokollant.create(prot,'dumpErhaltungsgroessen'); setlength(ns,felders[aktuelleFelder].matAnz); result:=true; s:=''; 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 ns[i]:=ns[i]*dX; s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]); {$IFDEF Dichteueberwachung} if (not abbruch) and (abs(ns[i]-felders[aktuelleFelder].massen[i])>0.5) then begin result:=false; pro.schreibe(' n['+inttostr(i)+'] ('+floattostr(ns[i])+') unterscheidet sich stark vom Anfangswert ('+floattostr(felders[aktuelleFelder].massen[i])+'), es scheinen sehr viele Teilchen entstanden oder verloren gegangen zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true); end; {$ENDIF} end; delete(s,1,1); prot.schreibe(s); pro.free; 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); prot.destroyall; halt; end; { tSimulation } constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string); var ifile: tMyStringlist; zeitverfahren: tZeitverfahren; s,t,aSuffix,aPrefix: string; deltaX,breite,pDNMax: extended; i: longint; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; pro: tProtokollant; na: pSigActionRec; begin inherited create; prot:=tProtokollant.create(protokollant,name); kvs:=tKnownValues.create; ifile:=tMyStringlist.create(prot,'ifile'); ifile.loadfromfile(inName); if not ifile.unfoldMacros then begin prot.schreibe('Fehlerhafte Macros in Parameterdatei '''+inName+'''!',true); prot.destroyall; halt(1); end; // Standardeinstellungen Bereich 'allgemein' zeitverfahren:=zfRungeKuttaVier; deltaX:=1e-2; kvs.add('λ',1); kvs.add('T',1); kvs.add('ω',2*pi); kvs.add('dX',deltaX); kvs.add('q',1); kvs.add('me',1); dT:=-1; sDT:=-1; sDX:=-1; tEnde:=100; breite:=10.0; pDNMax:=-1; fortschrittsAnzeige:=false; dTMin:=-1; gotSighup:=false; // Standardeinstellungen Bereich 'ausgaben' aPrefix:=extractfilepath(paramstr(0)); aSuffix:='.dat'; setlength(ausgabeDateien,0); // Standardeinstellungen Breich 'lichtVon...' lichter:=tMyStringlist.create(prot,'lichter'); // Standardeinstellungen Breich 'teilchen...' setlength(teilchen,0); repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in '''+inName+'''!',true); prot.destroyall; halt(1); end; if s='Dateiende' then break; if s='allgemein' then begin repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true); prot.destroyall; halt(1); end; if s='allgemeinEnde' then break; if s='runge-Kutta-3/8' then begin Zeitverfahren:=zfRungeKuttaDreiAchtel; continue; end; if s='runge-Kutta-4' then begin Zeitverfahren:=zfRungeKuttaVier; continue; end; if s='runge-Kutta-10' then begin Zeitverfahren:=zfRungeKuttaZehn; continue; end; if s='runge-Kutta-12' then begin Zeitverfahren:=zfRungeKuttaZwoelf; continue; end; if s='runge-Kutta-14' then begin Zeitverfahren:=zfRungeKuttaVierzehn; continue; end; if s='euler-Vorwärts' then begin Zeitverfahren:=zfEulerVorwaerts; continue; end; if startetMit('ortsschritt ',s) then begin deltaX:=exprtofloat(false,s,kvs,nil); kvs.add('dX',deltaX); continue; end; if startetMit('zeitschritt ',s) then begin dT:=exprtofloat(false,s,kvs,nil); kvs.add('dT',dT); continue; end; if startetMit('minZeitschritt ',s) then begin dTMin:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('diffusionsterm ',s) then begin pDNMax:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('zeit ',s) then begin 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); prot.destroyall; halt(1); until false; continue; end; if s='ausgaben' then begin repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); prot.destroyall; halt(1); end; if s='ausgabenEnde' then break; if startetMit('suffix ',s) then begin aSuffix:=s; continue; end; if startetMit('prefix ',s) then begin aPrefix:=s; continue; end; if startetMit('zeitschritt ',s) then 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,','); setlength(ausgabeDateien,length(ausgabeDateien)+1); ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot,'ausgabeDateien['+inttostr(length(ausgabeDateien)-1)+']'); end; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); prot.destroyall; halt(1); until false; 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); prot.destroyall; halt(1); end; if startetMit('teilchen',s) then begin if (s<>'') and (strtoint(s)<>length(teilchen)+1) then begin prot.schreibe('Ich erwarte die Teilchen beginnend bei 1 der Reihe nach in Parameterdatei '''+inName+'''!',true); prot.destroyall; halt(1); end; setlength(teilchen,length(teilchen)+1); teilchen[length(teilchen)-1]:=tTeilchenSpezies.create; repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich teilchen!',true); prot.destroyall; halt(1); end; if s='teilchen'+inttostr(length(teilchen))+'Ende' then break; if startetMit('spezifische Ladung ',s) then begin teilchen[length(teilchen)-1].spezifischeLadung:=-exprToFloat(false,s,kvs,nil); kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung); continue; end; if startetMit('maximaldichte ',s) then begin teilchen[length(teilchen)-1].nMax:=exprToFloat(false,s,kvs,nil); kvs.add('nmax'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMax); continue; end; if startetMit('minimaldichte ',s) then begin teilchen[length(teilchen)-1].nMin:=exprToFloat(false,s,kvs,nil); kvs.add('nmin'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMin); continue; end; if startetMit('verteilung ',s) then begin teilchen[length(teilchen)-1].liesDichteFunktion(s,ifile,kvs,teilchen,prot); continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich teilchen!',true); prot.destroyall; halt(1); until false; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true); prot.destroyall; halt(1); until false; ifile.free; if length(ausgabedateien)=0 then begin prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true); prot.destroyall; halt(1); end; if dT<0 then dT:=deltaX/10; if sDT<0 then sDT:=dT; sDT:=1/round(1/sDT); sT:=sDT/2; if sDX<0 then sDX:=deltaX; if dTMin<0 then dTMin:=min(1E-30,dT*1E-5); pro:=tProtokollant.create(prot,'create'); case zeitverfahren of zfEulerVorwaerts: pro.schreibe('Iteration mittels Euler-Vorwärts'); zfRungeKuttaDreiAchtel: pro.schreibe('Iteration mittels Runge-Kutta-3/8'); zfRungeKuttaVier: pro.schreibe('Iteration mittels Runge-Kutta-4'); zfRungeKuttaZehn: pro.schreibe('Iteration mittels Runge-Kutta-10'); zfRungeKuttaZwoelf: pro.schreibe('Iteration mittels Runge-Kutta-12'); zfRungeKuttaVierzehn: pro.schreibe('Iteration mittels Runge-Kutta-14'); 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 length(teilchen)>0 then begin pro.schreibe('teilchen:'); for i:=0 to length(teilchen)-1 do teilchen[i].dump(pro,' '); end; 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; if length(sighupSimulationen)=0 then begin new(na); na^.sa_Handler := sigActionHandler(@signalCapture); fillchar(na^.sa_Mask,sizeof(na^.sa_mask),#0); na^.sa_Flags := 0; {$ifdef Linux} // Linux specific na^.sa_Restorer := Nil; {$endif} if (fPSigaction(SIGHUP, na, nil) <> 0) then begin writeln(stderr, 'Error: ', fpgeterrno, '.'); halt(1); end; dispose(na); end; setlength(sighupSimulationen,length(sighupSimulationen)+1); sighupSimulationen[length(sighupSimulationen)-1]:=self; gitter:=tGitter.create(self,round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); lichter.free; end; 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 if sighupSimulationen[i]=self then begin for j:=i+1 to length(sighupSimulationen)-1 do sighupSimulationen[j-1]:=sighupSimulationen[j]; setlength(sighupSimulationen,length(sighupSimulationen)-1); end; inherited destroy; end; function tSimulation.iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? var i: longint; begin result:=false; if gitter.abbruch or (dT>sDT) then dT:=sDT; zeitPhysik:=zeitPhysik-now; if errorCode<2 then gitter.iteriereSchritt(dT,dTMin); zeitPhysik:=zeitPhysik+now; zeitDatei:=zeitDatei-now; while (gitter.t>=sT) and (sT