From a63868369a70801f4da5c9606e84dcf9ae172786 Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Fri, 6 Nov 2009 14:09:27 +0000 Subject: [PATCH] Temporary upload in order to fix an error in my calculations - some parts are not working at all --- .../classes/TBulkVortexFieldCalc.cpp | 905 +++++++++++++++++- .../classes/TFitPofBStartupHandler.cpp | 75 +- .../include/TBulkVortexFieldCalc.h | 62 +- .../include/TFitPofBStartupHandler.h | 5 +- .../TFitPofB-lib/test/TFitPofB_startup.xml | 4 +- .../TFitPofB-lib/test/testVortex-v2.cpp | 10 +- 6 files changed, 983 insertions(+), 78 deletions(-) diff --git a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp index 31e34b5f..1d5efd8f 100644 --- a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp @@ -40,10 +40,14 @@ using namespace std; #define PI 3.14159265358979323846 #define TWOPI 6.28318530717958647692 -const double fluxQuantum(2.067833667e7); // 10e14 times CGS units %% in CGS units should be 10^-7 - // in this case this is Gauss times square nm +const double fluxQuantum(2.067833667e7); // in this case this is Gauss times square nm + const double sqrt3(sqrt(3.0)); +//const double pi_4sqrt3(0.25*PI/sqrt3); +const double pi_4sqrt3(0.25*PI/sqrt(3.0)); +//const double pi_4sqrt3(PI*sqrt(3.0)); + double getXi(const double hc2) { // get xi given Hc2 in Gauss if (hc2) return sqrt(fluxQuantum/(TWOPI*hc2)); @@ -80,6 +84,40 @@ TBulkVortexFieldCalc::~TBulkVortexFieldCalc() { //fftw_cleanup_threads(); } +double TBulkVortexFieldCalc::GetBmin() const { + if (fGridExists) { + double min(fFFTout[0]); + unsigned int minindex(0), counter(0); + for (unsigned int j(0); j < fSteps * fSteps / 2; j++) { + if (fFFTout[j] <= 0.0) { + return 0.0; + } + if (fFFTout[j] < min) { + min = fFFTout[j]; + minindex = j; + } + counter++; + if (counter == fSteps/2) { // check only the first quadrant of B(x,y) + counter = 0; + j += fSteps/2; + } + } + return fFFTout[minindex]; + } else { + CalculateGrid(); + return GetBmin(); + } +} + +double TBulkVortexFieldCalc::GetBmax() const { + if (fGridExists) + return fFFTout[0]; + else { + CalculateGrid(); + return GetBmax(); + } +} + TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const string& wisdom, const unsigned int steps) { fWisdom = wisdom; if (steps % 2) { @@ -200,7 +238,7 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const { double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2])); double Hc2(getHc2(xi)); - double latConstTr(2.0*sqrt(fluxQuantum/(field*sqrt3))); + double latConstTr(sqrt(2.0*fluxQuantum/(field*sqrt3))); double xisq_2_scaled(2.0/3.0*pow(xi*PI/latConstTr,2.0)), lambdasq_scaled(4.0/3.0*pow(lambda*PI/latConstTr,2.0)); const int NFFT(fSteps); @@ -317,37 +355,846 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const { } -double TBulkTriVortexLondonFieldCalc::GetBmin() const { - if (fGridExists) { - double min(fFFTout[0]); - unsigned int minindex(0), counter(0); - for (unsigned int j(0); j < fSteps * fSteps / 2; j++) { - if (fFFTout[j] <= 0.0) { - return 0.0; +TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, const unsigned int steps) : fLatticeConstant(0.0), fSumAk(0.0), fSumOmegaSq(0.0) +{ + fWisdom = wisdom; + if (steps % 2) { + fSteps = steps + 1; + } else { + fSteps = steps; + } + fParam.resize(3); + fGridExists = false; + + int init_threads(fftw_init_threads()); + if (init_threads) + fftw_plan_with_nthreads(2); + + unsigned int stepsSq(fSteps*fSteps); + + fFFTout = new double[stepsSq]; + + fAkMatrix = new fftw_complex[stepsSq]; + fBkMatrix = new fftw_complex[stepsSq]; + fRealSpaceMatrix = new fftw_complex[stepsSq]; + fBMatrix = new double[stepsSq]; + fOmegaMatrix = new double[stepsSq]; + fOmegaSqMatrix = new double[stepsSq]; + fOmegaDiffXMatrix = new double[stepsSq]; + fOmegaDiffYMatrix = new double[stepsSq]; + fGMatrix = new double[stepsSq]; + fQxMatrix = new double[stepsSq]; + fQyMatrix = new double[stepsSq]; + fQxMatrixA = new double[stepsSq]; + fQyMatrixA = new double[stepsSq]; + + fCheckAkConvergence = new double[fSteps]; + fCheckBkConvergence = new double[fSteps]; + +// cout << "Check for the FFT plans..." << endl; + +// Load wisdom from file if it exists and should be used + + fUseWisdom = true; + int wisdomLoaded(0); + + FILE *wordsOfWisdomR; + wordsOfWisdomR = fopen(fWisdom.c_str(), "r"); + if (wordsOfWisdomR == NULL) { + fUseWisdom = false; + } else { + wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); + fclose(wordsOfWisdomR); + } + + if (!wisdomLoaded) { + fUseWisdom = false; + } + +// create the FFT plans + + if (fUseWisdom) { + fFFTplanAkToOmega = fftw_plan_dft_2d(fSteps, fSteps, fAkMatrix, fRealSpaceMatrix, FFTW_BACKWARD, FFTW_EXHAUSTIVE); + fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fRealSpaceMatrix, FFTW_BACKWARD, FFTW_EXHAUSTIVE); + fFFTplanOmegaToAk = fftw_plan_dft_2d(fSteps, fSteps, fRealSpaceMatrix, fAkMatrix, FFTW_FORWARD, FFTW_EXHAUSTIVE); + fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fRealSpaceMatrix, fBkMatrix, FFTW_FORWARD, FFTW_EXHAUSTIVE); + } + else { + fFFTplanAkToOmega = fftw_plan_dft_2d(fSteps, fSteps, fAkMatrix, fRealSpaceMatrix, FFTW_BACKWARD, FFTW_ESTIMATE); + fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fRealSpaceMatrix, FFTW_BACKWARD, FFTW_ESTIMATE); + fFFTplanOmegaToAk = fftw_plan_dft_2d(fSteps, fSteps, fRealSpaceMatrix, fAkMatrix, FFTW_FORWARD, FFTW_ESTIMATE); + fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fRealSpaceMatrix, fBkMatrix, FFTW_FORWARD, FFTW_ESTIMATE); + } +} + +TBulkTriVortexNGLFieldCalc::~TBulkTriVortexNGLFieldCalc() { + + // clean up + + fftw_destroy_plan(fFFTplanAkToOmega); + fftw_destroy_plan(fFFTplanBkToBandQ); + fftw_destroy_plan(fFFTplanOmegaToAk); + fftw_destroy_plan(fFFTplanOmegaToBk); + delete[] fAkMatrix; fAkMatrix = 0; + delete[] fBkMatrix; fBkMatrix = 0; + delete[] fRealSpaceMatrix; fRealSpaceMatrix = 0; + delete[] fBMatrix; fBMatrix = 0; + delete[] fOmegaMatrix; fOmegaMatrix = 0; + delete[] fOmegaSqMatrix; fOmegaSqMatrix = 0; + delete[] fOmegaDiffXMatrix; fOmegaDiffXMatrix = 0; + delete[] fOmegaDiffYMatrix; fOmegaDiffYMatrix = 0; + delete[] fGMatrix; fGMatrix = 0; + delete[] fQxMatrix; fQxMatrix = 0; + delete[] fQyMatrix; fQyMatrix = 0; + delete[] fQxMatrixA; fQxMatrixA = 0; + delete[] fQyMatrixA; fQyMatrixA = 0; + delete[] fCheckAkConvergence; fCheckAkConvergence = 0; + delete[] fCheckBkConvergence; fCheckBkConvergence = 0; +} + +void TBulkTriVortexNGLFieldCalc::CalculateIntermediateMatrices() const { + const int NFFT(fSteps); + const int NFFTsq(fSteps*fSteps); + +// Derivatives of omega + +// const double denomY(2.0*fLatticeConstant/static_cast(NFFT)); + const double denomY(2.0/static_cast(NFFT)); +// const double denomY(2.0); + const double denomX(denomY*sqrt3); + const double kappa(fParam[1]/fParam[2]); + const double fourKappaSq(4.0*kappa*kappa); + int i; + +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i = 0; i < NFFTsq; i += 1) { + if (!(i % NFFT)) { // first column + fOmegaDiffXMatrix[i] = (fOmegaMatrix[i+1]-fOmegaMatrix[i+NFFT-1])/denomX; + } else if (!((i + 1) % NFFT)) { // last column + fOmegaDiffXMatrix[i] = (fOmegaMatrix[i-NFFT+1]-fOmegaMatrix[i-1])/denomX; + } else { + fOmegaDiffXMatrix[i] = (fOmegaMatrix[i+1]-fOmegaMatrix[i-1])/denomX; + } + } + + for (i = 0; i < NFFT; i++) { // first row + fOmegaDiffYMatrix[i] = (fOmegaMatrix[i+NFFT]-fOmegaMatrix[NFFTsq-NFFT+i])/denomY; + } + + for (i = NFFTsq - NFFT; i < NFFTsq; i++) { // last row + fOmegaDiffYMatrix[i] = (fOmegaMatrix[i-NFFTsq+NFFT]-fOmegaMatrix[i-NFFT])/denomY; + } + +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i = NFFT; i < NFFTsq - NFFT; i++) { // the rest + fOmegaDiffYMatrix[i] = (fOmegaMatrix[i+NFFT]-fOmegaMatrix[i-NFFT])/denomY; + } + +// g-Matrix + +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i = 0; i < NFFTsq; i++) { + if (!fOmegaMatrix[i]) { + fGMatrix[i] = 0.0; + } else { + fGMatrix[i] = (fOmegaDiffXMatrix[i]*fOmegaDiffXMatrix[i]+fOmegaDiffYMatrix[i]*fOmegaDiffYMatrix[i])/(fourKappaSq*fOmegaMatrix[i]); + } + } + + return; +} + +void TBulkTriVortexNGLFieldCalc::ManipulateFourierCoefficients(fftw_complex* F, const double &coeff2, const double &coeff3) const { + const int NFFT(fSteps), NFFT_2(fSteps/2); + int lNFFT, l, k; + double Gsq; + double coeff1(4.0/3.0*pow(PI/fLatticeConstant,2.0)); + + for (l = 0; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + Gsq = coeff1*(static_cast(k*k) + 3.0*static_cast(l*l)); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] = 0.0; + F[lNFFT + k + 1][1] = 0.0; } - if (fFFTout[j] < min) { - min = fFFTout[j]; - minindex = j; + if (NFFT_2 % 2) { + Gsq = coeff1*(static_cast(k*k) + 3.0*static_cast(l*l)); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; } - counter++; - if (counter == fSteps/2) { // check only the first quadrant of B(x,y) - counter = 0; - j += fSteps/2; + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + Gsq = coeff1*(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast(l*l)); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] = 0.0; + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + Gsq = coeff1*(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast(l*l)); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; } } - return fFFTout[minindex]; - } else { - CalculateGrid(); - return GetBmin(); - } + + for (l = NFFT_2; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + Gsq = coeff1*(static_cast(k*k) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] = 0.0; + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + Gsq = coeff1*(static_cast(k*k) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + Gsq = coeff1*(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] = 0.0; + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + Gsq = coeff1*(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + F[lNFFT + k][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k][1] = 0.0; + } + } + + //intermediate rows + + for (l = 1; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + Gsq = coeff1*(static_cast((k+1)*(k+1)) + 3.0*static_cast(l*l)); + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + Gsq = coeff1*(static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast(l*l)); + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + } + } + + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + Gsq = coeff1*(static_cast((k+1)*(k+1)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + Gsq = coeff1*(static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + F[lNFFT + k + 1][0] *= coeff2/(Gsq+coeff3); + F[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] = 0.0; + F[lNFFT + k][1] = 0.0; + } + } + + F[0][0] = 0.0; + + return; } -double TBulkTriVortexLondonFieldCalc::GetBmax() const { - if (fGridExists) - return fFFTout[0]; - else { - CalculateGrid(); - return GetBmax(); - } +void TBulkTriVortexNGLFieldCalc::ManipulateFourierCoefficientsForQx(fftw_complex* F) const { + const int NFFT(fSteps), NFFT_2(fSteps/2); + int lNFFT, l, k; + double coeff1(1.5*fLatticeConstant/PI); + + for (l = 0; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + if (!k && !l) + F[0][0] = 0.0; + else + F[lNFFT + k][0] *= coeff1*static_cast(l)/(static_cast(k*k) + 3.0*static_cast(l*l)); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l)/(static_cast(k*k) + 3.0*static_cast(l*l)); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l)/(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast(l*l)); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l)/(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast(l*l)); + } + } + + for (l = NFFT_2; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l-NFFT)/(static_cast(k*k) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l-NFFT)/(static_cast(k*k) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l-NFFT)/(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(l-NFFT)/(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + } + + //intermediate rows + + for (l = 1; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(l)/(static_cast((k+1)*(k+1)) + 3.0*static_cast(l*l)); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(l)/(static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast(l*l)); + } + } + + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(l-NFFT)/(static_cast((k+1)*(k+1)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(l-NFFT)/(static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + } + + return; +} + +void TBulkTriVortexNGLFieldCalc::ManipulateFourierCoefficientsForQy(fftw_complex* F) const { + const int NFFT(fSteps), NFFT_2(fSteps/2); + int lNFFT, l, k; + double coeff1(0.5*sqrt3*fLatticeConstant/PI); + + for (l = 0; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + if (!k && !l) + F[0][0] = 0.0; + else + F[lNFFT + k][0] *= coeff1*static_cast(k)/(static_cast(k*k) + 3.0*static_cast(l*l)); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k)/(static_cast(k*k) + 3.0*static_cast(l*l)); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k-NFFT)/(static_cast((k-NFFT)*(k-NFFT))+3.0*static_cast(l*l)); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k-NFFT)/(static_cast((k-NFFT)*(k-NFFT))+3.0*static_cast(l*l)); + } + } + + for (l = NFFT_2; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k)/(static_cast(k*k)+3.0*static_cast((l-NFFT)*(l-NFFT))); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k)/(static_cast(k*k)+3.0*static_cast((l-NFFT)*(l-NFFT))); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k-NFFT)/(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + if (NFFT_2 % 2) { + F[lNFFT + k][0] *= coeff1*static_cast(k-NFFT)/(static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + } + + //intermediate rows + + for (l = 1; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(k+1)/(static_cast((k+1)*(k+1)) + 3.0*static_cast(l*l)); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(k+1-NFFT)/(static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast(l*l)); + } + } + + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(k+1)/(static_cast((k+1)*(k+1)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + F[lNFFT + k + 1][0] *= coeff1*static_cast(k+1-NFFT)/(static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT))); + } + } + + return; +} + +void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { + // SetParameters - method has to be called from the user before the calculation!! + if (fParam.size() < 3) { + cout << endl << "The SetParameters-method has to be called before B(x,y) can be calculated!" << endl; + return; + } + if (!fParam[0] || !fParam[1] || !fParam[2]) { + cout << endl << "The field, penetration depth and coherence length have to have finite values in order to calculate B(x,y)!" << endl; + return; + } + + double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2])); + double Hc2(getHc2(xi)), kappa(lambda/xi), S(TWOPI/(kappa*field)); + double scaledB(field/fluxQuantum*2.0*PI*xi*lambda); + +// fLatticeConstant = sqrt(2.0*fluxQuantum/(scaledB*sqrt3)); + fLatticeConstant = sqrt(2.0*TWOPI/(kappa*scaledB*sqrt3)); +// fLatticeConstant = sqrt(4.0*PI/(kappa*field*sqrt3)); +// fLatticeConstant = sqrt(2.0*S/sqrt3); +// double xisq_2_scaled(2.0/3.0*pow(xi*PI/latConstTr,2.0)), lambdasq_scaled(4.0/3.0*pow(lambda*PI/latConstTr,2.0)); + + const int NFFT(fSteps); + const int NFFT_2(fSteps/2); + const int NFFTsq(fSteps*fSteps); + + // fill the field Fourier components in the matrix + + // ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC + if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) { + int m; +#pragma omp parallel for default(shared) private(m) schedule(dynamic) + for (m = 0; m < NFFTsq; m++) { + fFFTout[m] = field; + } + // Set the flag which shows that the calculation has been done + fGridExists = true; + return; + } + + // ... now fill in the Fourier components (Abrikosov) if everything was okay above + double Gsq, sign; + int k, l, lNFFT; + +// omp causes problems with the fftw_complex*... comment it out for the moment +//#pragma omp parallel default(shared) private(l,lNFFT_2,k,Gsq) + { +// #pragma omp sections nowait + { +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = 0; l < NFFT_2; l += 2) { + if (!(l % 4)) { + sign = 1.0; + } else { + sign = -1.0; + } + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + sign = -sign; + Gsq = static_cast(k*k) + 3.0*static_cast(l*l); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = 0.0; + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + sign = -sign; + Gsq = static_cast(k*k) + 3.0*static_cast(l*l); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + } + + // sign = -sign; + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + sign = -sign; + Gsq = static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast(l*l); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = 0.0; + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + sign = -sign; + Gsq = static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast(l*l); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + } + } + +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = NFFT_2; l < NFFT; l += 2) { + if (!(l % 4)) { + sign = 1.0; + } else { + sign = -1.0; + } + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + sign = -sign; + Gsq = static_cast(k*k) + 3.0*static_cast((l-NFFT)*(l-NFFT)); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = 0.0; + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + sign = -sign; + Gsq = static_cast(k*k) + 3.0*static_cast((l-NFFT)*(l-NFFT)); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + } + + // sign = -sign; + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + sign = -sign; + Gsq = static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT)); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = 0.0; + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2) { + sign = -sign; + Gsq = static_cast((k-NFFT)*(k-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT)); + fAkMatrix[lNFFT + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k][1] = 0.0; + } + } + + // intermediate rows +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = 1; l < NFFT_2; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + Gsq = static_cast((k + 1)*(k + 1)) + 3.0*static_cast(l*l); + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2){ + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + Gsq = static_cast((k + 1-NFFT)*(k + 1-NFFT)) + 3.0*static_cast(l*l); + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2){ + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + } + } + +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT = l*NFFT; + for (k = 0; k < NFFT_2 - 1; k += 2) { + Gsq = static_cast((k+1)*(k+1)) + 3.0*static_cast((l-NFFT)*(l-NFFT)); + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2){ + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + } + + for (k = NFFT_2; k < NFFT - 1; k += 2) { + Gsq = static_cast((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast((l-NFFT)*(l-NFFT)); + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + fAkMatrix[lNFFT + k + 1][0] = exp(-pi_4sqrt3*Gsq); + fAkMatrix[lNFFT + k + 1][1] = 0.0; + } + if (NFFT_2 % 2){ + fAkMatrix[lNFFT + k][0] = 0.0; + fAkMatrix[lNFFT + k][1] = 0.0; + } + } + } /* end of sections */ + + } /* end of parallel section */ + + fAkMatrix[0][0] = 0.0; + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFT; l++) { + fCheckAkConvergence[l] = fAkMatrix[NFFT_2/2*NFFT + l][0]; + } + + fSumAk = 0.0; + for (l = 1; l < NFFTsq; l++) { + fSumAk += fAkMatrix[l][0]; + } + + cout << "fSumAk = " << fSumAk << endl; + + // Do the Fourier transform to get omega(x,y) - Abrikosov + + fftw_execute(fFFTplanAkToOmega); + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0]; + } + + CalculateIntermediateMatrices(); + + double denomQA; + +#pragma omp parallel for default(shared) private(l,denomQA) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + if (!fOmegaMatrix[l]) { + fQxMatrixA[l] = 0.0; + fQyMatrixA[l] = 0.0; + } else { + denomQA = 2.0*kappa*fOmegaMatrix[l]; + fQxMatrixA[l] = fOmegaDiffYMatrix[l]/denomQA; + fQyMatrixA[l] = -fOmegaDiffXMatrix[l]/denomQA; + } + fQxMatrix[l] = fQxMatrixA[l]; + fQyMatrix[l] = fQyMatrixA[l]; + } + + // initial B + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fBMatrix[l] = scaledB; + } + + bool akConverged(false), bkConverged(false), akInitiallyConverged(false), firstBkCalculation(true); + double coeff1, coeff2; + + while (!akConverged || !bkConverged) { + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fOmegaSqMatrix[l] = fOmegaMatrix[l]*fOmegaMatrix[l]; + } + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fRealSpaceMatrix[l][0] = fOmegaSqMatrix[l] + fOmegaMatrix[l]*(fQxMatrix[l]*fQxMatrix[l] + fQyMatrix[l]*fQyMatrix[l] - 2.0) + fGMatrix[l]; + fRealSpaceMatrix[l][1] = 0.0; + } + + fftw_execute(fFFTplanOmegaToAk); + + coeff1 = (4.0*kappa*kappa/static_cast(NFFTsq)); + coeff2 = (2.0*kappa*kappa); + + ManipulateFourierCoefficients(fAkMatrix, coeff1, coeff2); +// + fSumAk = 0.0; + for (l = 1; l < NFFTsq; l++) { + fSumAk += fAkMatrix[l][0]; + } + + cout << "fSumAk = " << fSumAk << endl; + + // Do the Fourier transform to get omega(x,y) + + fftw_execute(fFFTplanAkToOmega); + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0]; + } + + CalculateIntermediateMatrices(); +// +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fOmegaSqMatrix[l] = fOmegaMatrix[l]*fOmegaMatrix[l]; + } + + fSumOmegaSq = 0.0; + for (l = 0; l < NFFTsq; l++) { + fSumOmegaSq += fOmegaSqMatrix[l]; + } + + double coeffSum(0.0); + for (l = 0; l < NFFTsq; l++) { + coeffSum += fOmegaMatrix[l]*(1.0 - (fQxMatrix[l]*fQxMatrix[l] + fQyMatrix[l]*fQyMatrix[l])) - fGMatrix[l]; + } + + for (l = 0; l < NFFTsq; l++) { + fAkMatrix[l][0] *= coeffSum/fSumOmegaSq; + } + + akConverged = true; + + for (l = 0; l < NFFT; l++) { + if (fabs(fCheckAkConvergence[l] - fAkMatrix[NFFT_2/2*NFFT + l][0]) > 1.0E-12) { + akConverged = false; + break; + } + } + + if (!akConverged) { +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFT; l++) { + fCheckAkConvergence[l] = fAkMatrix[NFFT_2/2*NFFT + l][0]; + } + } else { + akInitiallyConverged = true; + } + + cout << "Ak Convergence: " << akConverged << endl; + + fSumAk = 0.0; + for (l = 1; l < NFFTsq; l++) { + fSumAk += fAkMatrix[l][0]; + } + + cout << "fSumAk = " << fSumAk << endl; + + // Do the Fourier transform to get omega(x,y) + + fftw_execute(fFFTplanAkToOmega); + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0]; + } + + CalculateIntermediateMatrices(); + + if (akInitiallyConverged) {// break; + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fRealSpaceMatrix[l][0] = (fOmegaMatrix[l]-fSumAk)*fBMatrix[l] + fSumAk*scaledB - fQxMatrix[l]*fOmegaDiffYMatrix[l] + fQyMatrix[l]*fOmegaDiffXMatrix[l]; + fRealSpaceMatrix[l][1] = 0.0; + } + + fftw_execute(fFFTplanOmegaToBk); + + if (firstBkCalculation) { +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFT; l++) { + fCheckBkConvergence[l] = fBkMatrix[NFFT_2/2*NFFT + l][0]; + } + firstBkCalculation = false; + akConverged = false; + } + + coeff1 = -2.0/static_cast(NFFTsq); + + ManipulateFourierCoefficients(fBkMatrix, coeff1, fSumAk); + + fftw_execute(fFFTplanBkToBandQ); + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fBMatrix[l] = scaledB + fRealSpaceMatrix[l][0]; + } + + bkConverged = true; + + for (l = 0; l < NFFT; l++) { + if (fabs(fCheckBkConvergence[l] - fBkMatrix[NFFT_2/2*NFFT + l][0]) > 1.0E-12) { + bkConverged = false; + break; + } + } + + cout << "Bk Convergence: " << bkConverged << endl; + + if (!bkConverged) { +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFT; l++) { + fCheckBkConvergence[l] = fBkMatrix[NFFT_2/2*NFFT + l][0]; + } + } else if (akConverged) { + break; + } + +// I need a copy of Bk... store it to Ak since already tooooooooooooooooooo much memory is allocated +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fAkMatrix[l][0] = fBkMatrix[l][0]; + fAkMatrix[l][1] = fBkMatrix[l][1]; + } + + ManipulateFourierCoefficientsForQx(fBkMatrix); + + fftw_execute(fFFTplanBkToBandQ); + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fQxMatrix[l] = fQxMatrixA[l] - fRealSpaceMatrix[l][1]; + } + +// Restore Bk... +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fBkMatrix[l][0] = fAkMatrix[l][0]; + fBkMatrix[l][1] = fAkMatrix[l][1]; + } + + ManipulateFourierCoefficientsForQy(fBkMatrix); + + fftw_execute(fFFTplanBkToBandQ); + +#pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fQyMatrix[l] = fQyMatrixA[l] + fRealSpaceMatrix[l][1]; + } + } + } + + // Set the flag which shows that the calculation has been done + + fGridExists = true; + return; + } diff --git a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp index 63b04c85..3fc0a024 100644 --- a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp +++ b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp @@ -97,7 +97,9 @@ void TFitPofBStartupHandler::OnEndDocument() */ void TFitPofBStartupHandler::OnStartElement(const char *str, const TList *attributes) { - if (!strcmp(str, "LEM")) { + if (!strcmp(str, "debug")) { + fKey = eDebug; + } else if (!strcmp(str, "LEM")) { fKey = eLEM; } else if (!strcmp(str, "VortexLattice")) { fKey = eVortex; @@ -115,8 +117,6 @@ void TFitPofBStartupHandler::OnStartElement(const char *str, const TList *attrib fKey = eNSteps; } else if (!strcmp(str, "N_VortexGrid")) { fKey = eGridSteps; - } else if (!strcmp(str, "VortexSymmetry")) { - fKey = eVortexSymmetry; } } @@ -144,6 +144,12 @@ void TFitPofBStartupHandler::OnEndElement(const char *str) void TFitPofBStartupHandler::OnCharacters(const char *str) { switch (fKey) { + case eDebug: + if (!strcmp(str, "1")) + fDebug = true; + else + fDebug = false; + break; case eLEM: fLEM = true; break; @@ -178,13 +184,6 @@ void TFitPofBStartupHandler::OnCharacters(const char *str) // convert str to int and assign it to the GridSteps-member fGridSteps = atoi(str); break; - case eVortexSymmetry: - // convert str to int and assign it to the VortexSymmetry-member (square = 0, triangular = 1) - if (!strcmp(str, "square")) - fVortexSymmetry = 2; - else - fVortexSymmetry = 1; - break; default: break; } @@ -270,46 +269,55 @@ void TFitPofBStartupHandler::CheckLists() // check if anything was set, and if not set some default stuff // check if delta_t is given, if not set default - cout << endl << "TFitPofBStartupHandler::CheckLists: check specified time resolution ... " << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check specified time resolution ... " << endl; if(!fDeltat) { cout << "TFitPofBStartupHandler::CheckLists: You did not specify the time resolution. Setting the default (10 ns)." << endl; fDeltat = 0.01; } else { - cout << fDeltat << " us" << endl; + if(fDebug) + cout << fDeltat << " us" << endl; } // check if delta_B is given, if not set default - cout << endl << "TFitPofBStartupHandler::CheckLists: check specified field resolution ..." << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check specified field resolution ..." << endl; if(!fDeltaB) { cout << "TFitPofBStartupHandler::CheckLists: You did not specify the field resolution. Setting the default (0.1 G)." << endl; fDeltaB = 0.1; } else { - cout << fDeltaB << " G" << endl; + if(fDebug) + cout << fDeltaB << " G" << endl; } // check if any wisdom-file is specified - cout << endl << "TFitPofBStartupHandler::CheckLists: check wisdom-file ..." << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check wisdom-file ..." << endl; if (!fWisdomFile.size()) { cout << "TFitPofBStartupHandler::CheckLists: You did not specify a wisdom file. No FFTW plans will be loaded or saved." << endl; fWisdomFile = ""; } else { - cout << fWisdomFile << endl; + if(fDebug) + cout << fWisdomFile << endl; } if (fLEM) { // check if any data path is given - cout << endl << "TFitPofBStartupHandler::CheckLists: check data path ..." << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check data path ..." << endl; if (!fDataPath.size()) { cout << "TFitPofBStartupHandler::CheckLists: This is not going to work, you have to set a valid data path where to find the rge-files in the xml-file!" << endl; exit(-1); } else { - cout << fDataPath << endl; + if(fDebug) + cout << fDataPath << endl; } // check if any energies are given - cout << endl << "TFitPofBStartupHandler::CheckLists: check energy list ..." << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check energy list ..." << endl; if (!fEnergyList.size()) { cout << "TFitPofBStartupHandler::CheckLists: Energy list empty! Setting the default list ( 0.0:0.1:32.9 keV)." << endl; char eChar[5]; @@ -320,39 +328,36 @@ void TFitPofBStartupHandler::CheckLists() } } } else { - for (unsigned int i (0); i < fEnergyList.size(); i++) - cout << fEnergyList[i] << " "; - cout << endl; + if(fDebug) { + for (unsigned int i (0); i < fEnergyList.size(); i++) + cout << fEnergyList[i] << " "; + cout << endl; + } } // check if any number of steps for the theory function is specified - cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for theory ..." << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for theory ..." << endl; if (!fNSteps) { cout << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the theory. Setting the default (3000)." << endl; fNSteps = 3000; } else { - cout << fNSteps << endl; + if(fDebug) + cout << fNSteps << endl; } } if (fVortex) { // check if any number of steps for the theory function is specified - cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex grid ..." << endl; + if(fDebug) + cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex grid ..." << endl; if (!fGridSteps) { cout << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the grid. Setting the default (256)." << endl; fGridSteps = 256; } else { - cout << fGridSteps << endl; - } - - // check if any number of steps for the theory function is specified - cout << endl << "TFitPofBStartupHandler::CheckLists: check symmetry of the Vortex lattice ..." << endl; - if (!fVortexSymmetry) { - cout << "TFitPofBStartupHandler::CheckLists: You did not specify the symmetry of the vortex lattice. Setting the default (triangular)." << endl; - fVortexSymmetry = 1; - } else { - cout << (fVortexSymmetry == 2 ? "square" : "triangular") << endl; + if(fDebug) + cout << fGridSteps << endl; } } } diff --git a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h index c5db18ec..82b3d9df 100644 --- a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h +++ b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h @@ -53,8 +53,8 @@ public: virtual double* DataB() const {return fFFTout;} virtual void SetParameters(const vector& par) {fParam = par; fGridExists = false;} virtual void CalculateGrid() const = 0; - virtual double GetBmin() const = 0; - virtual double GetBmax() const = 0; + virtual double GetBmin() const; + virtual double GetBmax() const; virtual bool GridExists() const {return fGridExists;} virtual void UnsetGridExists() const {fGridExists = false;} virtual unsigned int GetNumberOfSteps() const {return fSteps;} @@ -82,8 +82,62 @@ public: ~TBulkTriVortexLondonFieldCalc() {} void CalculateGrid() const; - double GetBmin() const; - double GetBmax() const; + +}; + +//-------------------- +// Class for triangular vortex lattice, Minimisation of the GL free energy à la Brandt +//-------------------- + +class TBulkTriVortexNGLFieldCalc : public TBulkVortexFieldCalc { + +public: + + TBulkTriVortexNGLFieldCalc(const string&, const unsigned int steps = 256); + ~TBulkTriVortexNGLFieldCalc(); + + void CalculateGrid() const; + fftw_complex* GetRealSpaceMatrix() const {return fRealSpaceMatrix;} + fftw_complex* GetAkMatrix() const {return fAkMatrix;} + fftw_complex* GetBkMatrix() const {return fBkMatrix;} + double* GetOmegaMatrix() const {return fOmegaMatrix;} + double* GetBMatrix() const {return fBMatrix;} + double* GetOmegaDiffXMatrix() const {return fOmegaDiffXMatrix;} + double* GetOmegaDiffYMatrix() const {return fOmegaDiffYMatrix;} + double* GetQxMatrix() const {return fQxMatrix;} + double* GetQyMatrix() const {return fQyMatrix;} + +private: + + void CalculateIntermediateMatrices() const; + void ManipulateFourierCoefficients(fftw_complex*, const double&, const double&) const; + void ManipulateFourierCoefficientsForQx(fftw_complex*) const; + void ManipulateFourierCoefficientsForQy(fftw_complex*) const; + + fftw_complex *fAkMatrix; + fftw_complex *fBkMatrix; + mutable fftw_complex *fRealSpaceMatrix; + mutable double *fBMatrix; + mutable double *fOmegaMatrix; + mutable double *fOmegaSqMatrix; + mutable double *fOmegaDiffXMatrix; + mutable double *fOmegaDiffYMatrix; + mutable double *fQxMatrix; + mutable double *fQyMatrix; + mutable double *fQxMatrixA; + mutable double *fQyMatrixA; + mutable double *fGMatrix; + + mutable double *fCheckAkConvergence; + mutable double *fCheckBkConvergence; + + mutable double fLatticeConstant; + mutable double fSumAk; + mutable double fSumOmegaSq; + fftw_plan fFFTplanAkToOmega; + fftw_plan fFFTplanBkToBandQ; + fftw_plan fFFTplanOmegaToAk; + fftw_plan fFFTplanOmegaToBk; }; diff --git a/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h b/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h index f47ad380..dcb70de9 100644 --- a/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h +++ b/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h @@ -66,13 +66,13 @@ class TFitPofBStartupHandler : public TQObject { virtual const string GetWisdomFile() const { return fWisdomFile; } virtual const unsigned int GetNSteps() const { return fNSteps; } virtual const unsigned int GetGridSteps() const { return fGridSteps; } - virtual const unsigned int GetVortexSymmetry() const { return fVortexSymmetry; } private: - enum EKeyWords {eEmpty, eComment, eLEM, eVortex, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps, eGridSteps, eVortexSymmetry}; + enum EKeyWords {eEmpty, eComment, eDebug, eLEM, eVortex, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps, eGridSteps}; EKeyWords fKey; + bool fDebug; bool fLEM; bool fVortex; string fDataPath; @@ -82,7 +82,6 @@ class TFitPofBStartupHandler : public TQObject { string fWisdomFile; unsigned int fNSteps; unsigned int fGridSteps; - unsigned int fVortexSymmetry; ClassDef(TFitPofBStartupHandler, 1) }; diff --git a/src/external/TFitPofB-lib/test/TFitPofB_startup.xml b/src/external/TFitPofB-lib/test/TFitPofB_startup.xml index ec3dd5ec..aa60aa77 100644 --- a/src/external/TFitPofB-lib/test/TFitPofB_startup.xml +++ b/src/external/TFitPofB-lib/test/TFitPofB_startup.xml @@ -9,6 +9,7 @@ * Parameters for bulk uSR data analysis (currently only one model to calculate field distributions for a triangular vortex lattice) Common parameters: + - debug : set it to 1 in order to obtain some information what is read from the xml-file - wisdom : sets the path to an FFTW-wisdom file - if there is no valid wisdom at the given path, wisdom handling will be disabled - delta_t : time resolution of P(t) in microseconds - delta_B : field resolution of P(B) in Gauss @@ -19,14 +20,13 @@ bulk parameters: - N_VortexGrid : determines the number of points used for the calculation of the vortex lattice field distribution (the grid will be N*N) - - VortexSymmetry : specify the vortex lattice symmetry - either "square" or "triangular" + 0 /home/l_wojek/analysis/WordsOfWisdom.dat 0.01 0.1 256 - triangular /home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000- diff --git a/src/external/TFitPofB-lib/test/testVortex-v2.cpp b/src/external/TFitPofB-lib/test/testVortex-v2.cpp index 18dab7e7..1626b43a 100644 --- a/src/external/TFitPofB-lib/test/testVortex-v2.cpp +++ b/src/external/TFitPofB-lib/test/testVortex-v2.cpp @@ -9,7 +9,7 @@ using namespace std; int main(){ - unsigned int NFFT(512); + unsigned int NFFT(256); vector parForVortex; parForVortex.resize(3); @@ -20,17 +20,17 @@ int main(){ vector parForPofB; parForPofB.push_back(0.01); //dt - parForPofB.push_back(2.0); //dB + parForPofB.push_back(0.1); //dB vector parForPofT; parForPofT.push_back(0.0); //phase parForPofT.push_back(0.01); //dt - parForPofT.push_back(2.0); //dB + parForPofT.push_back(0.1); //dB TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc("/home/l_wojek/analysis/WordsOfWisdom.dat", NFFT); - parForVortex[0] = 1000.0; //app.field - parForVortex[1] = 1000.0; //lambda + parForVortex[0] = 10.0; //app.field + parForVortex[1] = 200.0; //lambda parForVortex[2] = 4.0; //xi vortexLattice->SetParameters(parForVortex);