From 30567a80fd9df1b6e65fce56abc9cedc26f0217c Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Fri, 13 Nov 2015 11:12:17 +0100 Subject: Phasenraum Schnappschüsse sind nun korrekt orientiert MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Physikunit.pas | 121 ++++++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 99 insertions(+), 22 deletions(-) diff --git a/Physikunit.pas b/Physikunit.pas index b5000e2..e72c48a 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -127,6 +127,8 @@ type _dX,_dP,iDX: double; ffts: array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen fftTmp,lFftTmp: Pcomplex_double; // Zwischenspeicher für ffts + transps: array of fftw_plan_double; // Algorithmen zum Transponieren des Phasenraumes + transpTmp: pDouble; // Zwischenspeicher dafür procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline; {$IFDEF exzessiveArrayBereichsTests} procedure pruefeArrayEnden(fehler: string); @@ -157,6 +159,7 @@ type {$INCLUDE linearkombinationen.inc} {$UNDEF LiKoInterface} procedure schreibePhasenraum(tlc: longint; datNamPre, datNamSuf: string); + function transponiereFuer(tlc: longint): pDouble; end; { tGitter } @@ -191,6 +194,7 @@ type fortschrittsAnzeige: boolean; feldAusgabeDateien: array of tFeldAusgabeDatei; impulsRaumAusgabeDateien: array of tImpulsRaumAusgabeDatei; + transpTmp: array of double; public gotSigusr1,gotSigterm, gotSigint: boolean; @@ -557,21 +561,22 @@ begin assignFile(datei,pre+'-'+inttostr(nNum-1)+'-'+inttostr(tCnt)+suf); rewrite(datei,1); + inc(tCnt); sF:=gitter.felders[gitter.aktuelleFelder]; - blockWrite(datei,gitter.xl,sizeof(double)); - tmpD:=gitter.xl + (sF.aX-1)*sF.dX; - blockWrite(datei,tmpD,sizeof(double)); - tmpI:=sF.aX; - blockWrite(datei,tmpI,sizeof(longint)); tmpD:=-(sF.aP-1)*sF.dP/2; blockWrite(datei,tmpD,sizeof(double)); tmpD:=-tmpD; blockWrite(datei,tmpD,sizeof(double)); tmpI:=sF.aP; blockWrite(datei,tmpI,sizeof(longint)); - blockWrite(datei,sF.impulsraum[teilchen,false]^,sizeof(double)*sF.aX*sF.aP); + blockWrite(datei,gitter.xl,sizeof(double)); + tmpD:=gitter.xl + (sF.aX-1)*sF.dX; + blockWrite(datei,tmpD,sizeof(double)); + tmpI:=sF.aX; + blockWrite(datei,tmpI,sizeof(longint)); + blockWrite(datei,sF.transponiereFuer(teilchen)^,sizeof(double)*sF.aX*sF.aP); closeFile(datei); end; @@ -742,6 +747,7 @@ var emF: tEMFeldGroesze; emQ: tEMQuellGroesze; maF: tMaterieFeldGroesze; + dims: Pfftw_iodim; begin inherited create; gitter:=parent; @@ -871,25 +877,60 @@ begin ); end; + (* + fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, + double *in, const int *inembed, + int istride, int idist, + fftw_complex *out, const int *onembed, + int ostride, int odist, + unsigned flags); + fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany, + fftw_complex *in, const int *inembed, + int istride, int idist, + double *out, const int *onembed, + int ostride, int odist, + unsigned flags); + *) + + fillchar(transps,sizeof(transps),#0); + setlength(transps,length(teilchen)); + fftw_getmem(transpTmp,(aX*aP+1)*sizeof(double)); + (transpTmp+aX*aP)^:=pi; + getmem(dims,2*sizeof(fftw_iodim)); + dims^.n:=aX; // aX Spalten + dims^.ins:=aP; // Werte im Input um jeweils aP versetzt + dims^.ous:=1; // und im Output um jeweils 1 + (dims+1)^.n:=aP; // aP Zeilen + (dims+1)^.ins:=1; // Werte im Input um jeweils 1 versetzt + (dims+1)^.ous:=aX; // und im Output um jeweils aX + for i:=0 to length(transps)-1 do begin + transps[i]:= + fftw_plan_guru_r2r( + 0, // null Dimensionen für die FFT, also nur kopieren! + nil, // Bockwurst + 2, // 2 Dimensionen + dims, // @(aX;aP) + impulsraum[i,false], // von hier + transpTmp, // nach dort + nil, // Bockwurst + [fftw_estimate] // einfach irgendwas nehmen + ); + end; + freemem(dims); + + (* + fftw_plan fftw_plan_guru_r2r(int rank, const fftw_iodim *dims, + int howmany_rank, + const fftw_iodim *howmany_dims, + double *in, double *out, + const fftw_r2r_kind *kind, + unsigned flags); + *) + {$IFDEF exzessiveArrayBereichsTests} - pruefeArrayEnden('FFT schreibt beim Planen über das Ende hinaus!'); + pruefeArrayEnden('FFT/Transponieren schreibt beim Planen über das Ende hinaus!'); {$ENDIF} - (* - fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, - double *in, const int *inembed, - int istride, int idist, - fftw_complex *out, const int *onembed, - int ostride, int odist, - unsigned flags); - fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany, - fftw_complex *in, const int *inembed, - int istride, int idist, - double *out, const int *onembed, - int ostride, int odist, - unsigned flags); - *) - fillchar(massen,sizeof(massen),#0); setlength(massen,length(teilchen)); for i:=0 to length(massen)-1 do @@ -950,6 +991,10 @@ begin fftw_destroy_plan(ffts[i,j,abl]); setlength(ffts,0); fftw_freemem(fftTmp); + for i:=0 to length(transps)-1 do + fftw_destroy_plan(transps[i]); + setlength(transps,0); + fftw_freemem(transpTmp); inherited destroy; end; @@ -974,6 +1019,11 @@ var i,j: longint; abl: boolean; begin + if abs((transpTmp+aP*aX)^-pi)>1e-10 then begin + gitter.prot.schreibe(fehler); + gitter.prot.schreibe(floattostr((transpTmp+aP*aX)^)+' vs. '+floattostr(pi)); + raise exception.create('Fehler in tFelder.pruefeArrayEnden!'); + end; if (abs(lFftTmp^.re-pi)>1e-10) or (abs(lFftTmp^.im-sqr(pi))>1e-10) then begin gitter.prot.schreibe(fehler); gitter.prot.schreibe(floattostr(lFftTmp^.re)+' vs. '+floattostr(pi)+' und '+floattostr(lFftTmp^.im)+' vs. '+floattostr(sqr(pi))); @@ -1423,6 +1473,29 @@ begin end; end; +function tFelder.transponiereFuer(tlc: longint): pDouble; +{$IFDEF exzessiveArrayBereichsTests} +var + i: longint; +{$ENDIF} +begin + {$IFDEF exzessiveArrayBereichsTests} + for i:=0 to aP*aX-1 do + (transpTmp+i)^:=sqrt(2); + {$ENDIF} + fftw_execute(transps[tlc]); + {$IFDEF exzessiveArrayBereichsTests} + pruefeArrayEnden('TransporniereFuer('+inttostr(tlc)+') schreibt über das Ende hinaus!'); + for i:=0 to aP*aX-1 do + if (transpTmp+i)^=sqrt(2) then begin + gitter.prot.schreibe('TransporniereFuer('+inttostr(tlc)+') schreibt nicht alle Werte!'); + gitter.prot.schreibe('Position '+inttostr(i)+' wurde ausgelassen!'); + raise exception.create('Fehler in tFelder.transponiereFuer!'); + end; + {$ENDIF} + result:=transpTmp; +end; + // tGitter ********************************************************************* constructor tGitter.create(besitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string); @@ -1608,6 +1681,9 @@ begin setlength(feldAusgabeDateien,0); setlength(impulsRaumAusgabeDateien,0); + fillchar(transpTmp,sizeof(transpTmp),#0); // Transpositions-Puffer initialisieren + setlength(transpTmp,0); + // Standardeinstellungen Bereich 'lichtVon...' lichter:=tMyStringlist.create(prot,'lichter'); @@ -1941,6 +2017,7 @@ begin simulationen[j-1]:=simulationen[j]; setlength(simulationen,length(simulationen)-1); end; + setlength(transpTmp,0); inherited destroy; end; -- cgit v1.2.3-54-g00ecf