From 9a20545e798a68ebd73f92ca0f19a2cfbebf5983 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Tue, 4 Aug 2015 15:55:25 +0200 Subject: alle Runge-Kuttas in .inc ausgelagert, minimales dT auf dX gesetzt (sonst Bug z.B. in AY) --- Physikunit.pas | 202 +++++++++++++++++++------------------------------- Plasmapropagation.lps | 114 ++++++++++++++-------------- input.plap | 10 +-- rk108.inc | 64 ++++++++-------- rk1210.inc | 96 ++++++++++++------------ rk1412.inc | 136 ++++++++++++++++----------------- rk3_8.inc | 47 ++++++++++++ rk3_8.txt | 30 ++++++++ rk4.inc | 41 ++++++++++ rk4.txt | 30 ++++++++ rktopas | 18 ++--- 11 files changed, 448 insertions(+), 340 deletions(-) create mode 100644 rk3_8.inc create mode 100644 rk3_8.txt create mode 100644 rk4.inc create mode 100644 rk4.txt diff --git a/Physikunit.pas b/Physikunit.pas index 0a2ca5e..142ee4b 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -52,7 +52,7 @@ type destructor destroy; override; function dump: string; procedure schreibeKopf; - procedure speichereWerte(gitter: tGitter; sDX: extended); + procedure speichereWerte(gitter: tGitter; sT, sDX: extended); end; { tTeilchenSpezies } @@ -121,14 +121,14 @@ type matAnz: longint; spezLadungen: tExtendedArray; lichters: array[boolean] of tMyStringlist; - procedure setzeRaender(dT,iDX,iDT: extended); + 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(dT,dX,iDT,iDX,pDNMax: extended); + 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; @@ -163,7 +163,8 @@ type dX,iDX,pDNMax,dTMaximum,xl,t: extended; kvs: tKnownValues; procedure diffusionsTermAnpassen(pro: tProtokollant); - function pruefeMaxDT(aF: longint; var mDT,dT: extended; dTMin: extended): boolean; // wurde verkleinert? + function pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert? + procedure abbrechen; public aktuelleFelder: longint; @@ -334,7 +335,7 @@ begin closefile(datei); end; -procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sDX: extended); +procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: extended); var i,cnt: longint; sX,cX: extended; begin @@ -348,7 +349,7 @@ begin if (teilchen>=0) and (matInhalt=mfP) then gitter.berechne('P',teilchen); - if gitter.t>=nNum+sDT/2 then + if sT>=nNum+sDT/2 then schreibeKopf; cnt:=floor((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)/sDX*gitter.dX+1); @@ -749,12 +750,11 @@ begin inherited destroy; end; -procedure tFelder.setzeRaender(dT,iDX,iDT: extended); +procedure tFelder.setzeRaender(iDX: extended); var emF: tEMFeldInhalt; rechts: boolean; i: longint; - nVal: extended; begin for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld inhalt[0].emWerte[emF,true]:= @@ -766,17 +766,13 @@ begin 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 begin - par.kvs.add('t',par.t+dT); - nVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil); - par.kvs.add('t',par.t); + 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] + - (nVal - exprToFloat(false,lichters[rechts][i],par.kvs,nil))*iDT; - end; + exprToFloat(false,lichters[rechts][i],par.kvs,nil); end; -procedure tFelder.berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended); +procedure tFelder.berechneAbleitungen(dX,iDX,pDNMax: extended); var i: longint; begin @@ -795,7 +791,7 @@ begin inhalt[1].berechneNAbleitung(iDX,pDNMax,false); inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true); - setzeRaender(dT,iDX,iDT); + setzeRaender(iDX); for i:=1 to length(inhalt)-2 do inhalt[i].berechneAbleitungen(dX,iDX); @@ -863,49 +859,47 @@ end; procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant); var - i,j: longint; + i,j: longint; + curMax: extended; begin + curMax:=0; for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do for j:=0 to felders[aktuelleFelder].matAnz-1 do - if 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-10)) - > pDNMax then begin - if pDNMaxDynamisch then begin - pDNMax:= - 2* - abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]* + curMax:= + max(curMax, + 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-10)); - 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), die numerische Stabilität ist nicht mehr gewährleistet!',true); - pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true); - pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true); - pDNMax:=-1; - end; - end; + max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-100))); + 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), die numerische Stabilität ist nicht mehr gewährleistet!',true); + pro.schreibe(' Lösung: größeren Diffusionsterm wählen',true); + pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern',true); + pDNMax:=0; // Term komplett abschalten + end; + end; end; -function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended; dTMin: extended): boolean; // wurde verkleinert? +function tGitter.pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert? begin - mDT:=min(mDT,felders[aF].maxDT); - // das maximale dT, welches bei den Momentanen Ableitungen nicht zu + 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:=mDT '+floattostr(dT)); + prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')'); if dT> ${outputFile} - echo -ne "${lines}" | grep "^${i} \+" | sed "s/^[0-9]\+ \+\([0-9]\+\) \+[-0-9.*^]\+$/,\n felders[\1]/" | tr "\n" ";" | sed "s/;,/,/g" | sed "s/;$//" | tr ";" "\n" | sed "s/felders\[1\]/felders[aktuelleFelder]/g" >> ${outputFile} - echo -ne "${lines}" | grep "^${i} \+" | sed "s/^[0-9]\+ \+[0-9]\+ \+\([-0-9.*^]\+\)$/,\n \1 * dT/" | tr "\n" ";" | sed "s/,;/,\n/g" | sed "s/;,/,/g" | tr ";" " " | sed "s/ $//" | sed "s/\*10^/E/g" >> ${outputFile} + echo -ne "${lines}" | grep "^${i} \+" | sed "s/^[0-9]\+ \+\([0-9]\+\) \+[-0-9.*^/]\+$/,\n felders[\1]/" | tr "\n" ";" | sed "s/;,/,/g" | sed "s/;$//" | tr ";" "\n" | sed "s/felders\[1\]/felders[aktuelleFelder]/g" >> ${outputFile} + echo -ne "${lines}" | grep "^${i} \+" | sed "s/^[0-9]\+ \+[0-9]\+ \+\([-0-9.*^/]\+\)$/,\n \1 * dT/" | tr "\n" ";" | sed "s/,;/,\n/g" | sed "s/;,/,/g" | tr ";" " " | sed "s/ $//" | sed "s/\*10^/E/g" >> ${outputFile} echo -e "\n );" >> ${outputFile} - echo " felders[${i}].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax);" >> ${outputFile} + echo " felders[${i}].berechneAbleitungen(dX,iDX,pDNMax);" >> ${outputFile} echo >> ${outputFile} - echo " if pruefeMaxDT(${i},mDT,dT,dTMin) then" >> ${outputFile} + echo " if pruefeMaxDT(${i},dTMax,dT,dTMin) then" >> ${outputFile} echo " continue;" >> ${outputFile} echo >> ${outputFile} done - rest=$[$(wc -l ${inputFile} | cut -d " " -f 1)-$(grep --text -n "^ *k\s\+c\[k].$" ${inputFile} | cut -d ":" -f 1)] - lines=$(tail -n${rest} ${inputFile} | sed "s/^ \+//" | sed "s/[^0-9]$//" | grep --text "^[0-9]\+ \+[-0-9.^*]\+$" | grep -v "^[0-9]\+ \+0\(\.0\+\)\?$" | sed "s/0\+$//" | awk '{print $1+1" "$2}') + rest=$[$(wc -l ${inputFile} | cut -d " " -f 1)-$(grep --text -n "^ *k\s\+c\[k].\?$" ${inputFile} | cut -d ":" -f 1)] + lines=$(tail -n${rest} ${inputFile} | sed "s/^ \+//" | sed "s/[^0-9]$//" | grep --text "^[0-9]\+ \+[-0-9.^*/]\+$" | grep -v "^[0-9]\+ \+0\(\.0\+\)\?$" | sed "s/0\+$//" | awk '{print $1+1" "$2}') echo -e " felders[1-aktuelleFelder].liKo(\n felders[aktuelleFelder]," >> ${outputFile} - echo -e "${lines}" | sed "s/^\([0-9]\+\) \+[-0-9.*^]\+$/ felders[\1],/" | sed "s/felders\[1\]/felders[aktuelleFelder]/g" >> ${outputFile} - echo -ne "${lines}" | sed "s/^[0-9]\+ \+\([-0-9.*^]\+\)$/ \1 * dT,/" | tr "\n" "," | sed "s/\([^,]\),$/\1\n );\n/" | sed "s/,,/,\n/g" | sed "s/\*10^/E/g" >> ${outputFile} + echo -e "${lines}" | sed "s/^\([0-9]\+\) \+[-0-9.*^/]\+$/ felders[\1],/" | sed "s/felders\[1\]/felders[aktuelleFelder]/g" >> ${outputFile} + echo -ne "${lines}" | sed "s/^[0-9]\+ \+\([-0-9.*^/]\+\)$/ \1 * dT,/" | tr "\n" "," | sed "s/\([^,]\),$/\1\n );\n/" | sed "s/,,/,\n/g" | sed "s/\*10^/E/g" >> ${outputFile} done -- cgit v1.2.3-70-g09d2