diff --git a/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp index 1bf8ece4..a8a93032 100644 --- a/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp @@ -673,13 +673,14 @@ TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, con fOmegaMatrix = new double[stepsSq]; // |psi|^2 (x,y) fOmegaDiffMatrix = new fftw_complex[stepsSq]; // grad omega - fFFTin = new fftw_complex[(fSteps/2 + 1)*fSteps]; // aK matrix + fFFTin = new fftw_complex[stepsSq]; // aK matrix fBkMatrix = new fftw_complex[stepsSq]; // bK matrix + fRealSpaceMatrix = new fftw_complex[stepsSq]; fQMatrix = new fftw_complex[stepsSq]; fQMatrixA = new fftw_complex[stepsSq]; - fCheckAkConvergence = new double[fSteps/2 + 1]; + fCheckAkConvergence = new double[fSteps]; fCheckBkConvergence = new double[fSteps]; // Load wisdom from file if it exists and should be used @@ -704,16 +705,16 @@ TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, con if (fUseWisdom) { // use the first plan from the base class here - it will be destroyed by the base class destructor - fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fOmegaMatrix, FFTW_EXHAUSTIVE); + fFFTplan = fftw_plan_dft_2d(fSteps, fSteps, fFFTin, fRealSpaceMatrix, FFTW_BACKWARD, FFTW_EXHAUSTIVE); fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_BACKWARD, FFTW_EXHAUSTIVE); - fFFTplanOmegaToAk = fftw_plan_dft_r2c_2d(fSteps, fSteps, fOmegaMatrix, fFFTin, FFTW_EXHAUSTIVE); + fFFTplanOmegaToAk = fftw_plan_dft_2d(fSteps, fSteps, fRealSpaceMatrix, fFFTin, FFTW_FORWARD, FFTW_EXHAUSTIVE); fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_FORWARD, FFTW_EXHAUSTIVE); } else { // use the first plan from the base class here - it will be destroyed by the base class destructor - fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fOmegaMatrix, FFTW_ESTIMATE); + fFFTplan = fftw_plan_dft_2d(fSteps, fSteps, fFFTin, fRealSpaceMatrix, FFTW_BACKWARD, FFTW_ESTIMATE); fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_BACKWARD, FFTW_ESTIMATE); - fFFTplanOmegaToAk = fftw_plan_dft_r2c_2d(fSteps, fSteps, fOmegaMatrix, fFFTin, FFTW_ESTIMATE); + fFFTplanOmegaToAk = fftw_plan_dft_2d(fSteps, fSteps, fRealSpaceMatrix, fFFTin, FFTW_FORWARD, FFTW_ESTIMATE); fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_FORWARD, FFTW_ESTIMATE); } } @@ -730,6 +731,7 @@ TBulkTriVortexNGLFieldCalc::~TBulkTriVortexNGLFieldCalc() { delete[] fOmegaDiffMatrix; fOmegaDiffMatrix = 0; delete[] fBkMatrix; fBkMatrix = 0; + delete[] fRealSpaceMatrix; fRealSpaceMatrix = 0; delete[] fQMatrix; fQMatrix = 0; delete[] fQMatrixA; fQMatrixA = 0; @@ -745,6 +747,117 @@ void TBulkTriVortexNGLFieldCalc::CalculateGradient() const { const int NFFT_2(fSteps/2); const int NFFTsq(fSteps*fSteps); + int i, j, l, index; + + // Take the derivative of the Fourier sum of omega + + // First save a copy of the real aK-matrix in the imaginary part of the bK-matrix + #pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; ++l) { + fBkMatrix[l][1] = fFFTin[l][0]; + } + + // dw/dx = sum_K aK Kx sin(Kx*x + Ky*y) + // First multiply the aK with Kx, then call FFTW + + const double coeffKx(TWOPI/(sqrt3*fLatticeConstant)); + + // even rows + for (i = 0; i < NFFT; i += 2) { + // j = 0 + fFFTin[NFFT*i][0] = 0.0; + // j != 0 + for (j = 2; j < NFFT_2; j += 2) { + fFFTin[(j + NFFT*i)][0] *= coeffKx*static_cast(j); + } + for (j = NFFT_2; j < NFFT; j += 2) { + fFFTin[(j + NFFT*i)][0] *= coeffKx*static_cast(j - NFFT); + } + } + + // odd rows + for (i = 1; i < NFFT; i += 2) { + for (j = 0; j < NFFT_2; j += 2) { + fFFTin[(j + 1 + NFFT*i)][0] *= coeffKx*static_cast(j + 1); + } + for (j = NFFT_2; j < NFFT; j += 2) { + fFFTin[(j + 1 + NFFT*i)][0] *= coeffKx*static_cast(j + 1 - NFFT); + } + } + + fftw_execute(fFFTplan); + + // Copy the results to the gradient matrix and restore the original aK-matrix + #pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; ++l) { + fOmegaDiffMatrix[l][0] = fRealSpaceMatrix[l][1]; + fFFTin[l][0] = fBkMatrix[l][1]; + } + + // dw/dy = sum_K aK Ky sin(Kx*x + Ky*y) + // First multiply the aK with Ky, then call FFTW + + const double coeffKy(TWOPI/fLatticeConstant); + double ky; + + // even rows + // i = 0 + for (j = 0; j < NFFT; j += 2) { + fFFTin[j][0] = 0.0; + } + // i != 0 + for (i = 2; i < NFFT_2; i += 2) { + ky = coeffKy*static_cast(i); + for (j = 0; j < NFFT; j += 2) { + fFFTin[(j + NFFT*i)][0] *= ky; + } + } + for (i = NFFT_2; i < NFFT; i += 2) { + ky = coeffKy*static_cast(i - NFFT); + for (j = 0; j < NFFT; j += 2) { + fFFTin[(j + NFFT*i)][0] *= ky; + } + } + + // odd rows + for (i = 1; i < NFFT_2; i += 2) { + ky = coeffKy*static_cast(i); + for (j = 0; j < NFFT; j += 2) { + fFFTin[(j + 1 + NFFT*i)][0] *= ky; + } + } + for (i = NFFT_2 + 1; i < NFFT; i += 2) { + ky = coeffKy*static_cast(i - NFFT); + for (j = 0; j < NFFT; j += 2) { + fFFTin[(j + 1 + NFFT*i)][0] *= ky; + } + } + + fftw_execute(fFFTplan); + + // Copy the results to the gradient matrix and restore the original aK-matrix + #pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; ++l) { + fOmegaDiffMatrix[l][1] = fRealSpaceMatrix[l][1]; + fFFTin[l][0] = fBkMatrix[l][1]; + fBkMatrix[l][1] = 0.0; + } + + // Ensure that omega at the vortex-core positions is zero + fOmegaMatrix[0] = 0.0; + fOmegaMatrix[(NFFT+1)*NFFT_2] = 0.0; + + // Ensure that the derivatives at the vortex-core positions are zero + fOmegaDiffMatrix[0][0] = 0.0; + fOmegaDiffMatrix[(NFFT+1)*NFFT_2][0] = 0.0; + fOmegaDiffMatrix[0][1] = 0.0; + fOmegaDiffMatrix[(NFFT+1)*NFFT_2][1] = 0.0; + +/* + const int NFFT(fSteps); + const int NFFT_2(fSteps/2); + const int NFFTsq(fSteps*fSteps); + const double denomY(2.0*fLatticeConstant/static_cast(NFFT)); const double denomX(denomY*sqrt3); @@ -792,7 +905,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGradient() const { fOmegaDiffMatrix[(NFFT+1)*NFFT_2][0] = 0.0; fOmegaDiffMatrix[0][1] = 0.0; fOmegaDiffMatrix[(NFFT+1)*NFFT_2][1] = 0.0; - +*/ return; } @@ -808,9 +921,9 @@ void TBulkTriVortexNGLFieldCalc::FillAbrikosovCoefficients() const { } else { sign = -1.0; } - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast(l*l); - for (k = 0; k < NFFT_2 - 1; k += 2) { + for (k = 0; k < NFFT_2; k += 2) { sign = -sign; Gsq = static_cast(k*k) + ll; fFFTin[lNFFT_2 + k][0] = sign*exp(-pi_4sqrt3*Gsq); @@ -818,11 +931,14 @@ void TBulkTriVortexNGLFieldCalc::FillAbrikosovCoefficients() const { fFFTin[lNFFT_2 + k + 1][0] = 0.0; fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - sign = -sign; - Gsq = static_cast(k*k) + ll; - fFFTin[lNFFT_2 + k][0] = sign*exp(-pi_4sqrt3*Gsq); - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + sign = -sign; + Gsq = static_cast((k-NFFT)*(k-NFFT)) + ll; + fFFTin[lNFFT_2 + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } for (l = NFFT_2; l < NFFT; l += 2) { @@ -831,9 +947,9 @@ void TBulkTriVortexNGLFieldCalc::FillAbrikosovCoefficients() const { } else { sign = -1.0; } - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast((l-NFFT)*(l-NFFT)); - for (k = 0; k < NFFT_2 - 1; k += 2) { + for (k = 0; k < NFFT_2; k += 2) { sign = -sign; Gsq = static_cast(k*k) + ll; fFFTin[lNFFT_2 + k][0] = sign*exp(-pi_4sqrt3*Gsq); @@ -841,31 +957,38 @@ void TBulkTriVortexNGLFieldCalc::FillAbrikosovCoefficients() const { fFFTin[lNFFT_2 + k + 1][0] = 0.0; fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - sign = -sign; - Gsq = static_cast(k*k) + ll; - fFFTin[lNFFT_2 + k][0] = sign*exp(-pi_4sqrt3*Gsq); - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT_2; k += 2) { + sign = -sign; + Gsq = static_cast((k-NFFT)*(k-NFFT)) + ll; + fFFTin[lNFFT_2 + k][0] = sign*exp(-pi_4sqrt3*Gsq); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } // intermediate rows for (l = 1; l < NFFT_2; l += 2) { - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast(l*l); - for (k = 0; k < NFFT_2 - 1; k += 2) { + for (k = 0; k < NFFT_2; k += 2) { Gsq = static_cast((k + 1)*(k + 1)) + ll; fFFTin[lNFFT_2 + k][0] = 0.0; fFFTin[lNFFT_2 + k][1] = 0.0; fFFTin[lNFFT_2 + k + 1][0] = exp(-pi_4sqrt3*Gsq); fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - fFFTin[lNFFT_2 + k][0] = 0.0; - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + Gsq = static_cast((k + 1 - NFFT)*(k + 1 - NFFT)) + ll; + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = exp(-pi_4sqrt3*Gsq); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } for (l = NFFT_2 + 1; l < NFFT; l += 2) { - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast((l-NFFT)*(l-NFFT)); for (k = 0; k < NFFT_2 - 1; k += 2) { Gsq = static_cast((k+1)*(k+1)) + ll; @@ -874,9 +997,13 @@ void TBulkTriVortexNGLFieldCalc::FillAbrikosovCoefficients() const { fFFTin[lNFFT_2 + k + 1][0] = exp(-pi_4sqrt3*Gsq); fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - fFFTin[lNFFT_2 + k][0] = 0.0; - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + Gsq = static_cast((k+1 - NFFT)*(k+1 - NFFT)) + ll; + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = exp(-pi_4sqrt3*Gsq); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } fFFTin[0][0] = 0.0; @@ -896,23 +1023,26 @@ void TBulkTriVortexNGLFieldCalc::ManipulateFourierCoefficientsA() const { double Gsq, ll; for (l = 0; l < NFFT_2; l += 2) { - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast(l*l); - for (k = 0; k < NFFT_2 - 1; k += 2) { + for (k = 0; k < NFFT_2; k += 2) { Gsq = coeff1*(static_cast(k*k) + ll); fFFTin[lNFFT_2 + k][0] *= coeff2/(Gsq+coeff3); fFFTin[lNFFT_2 + k][1] = 0.0; fFFTin[lNFFT_2 + k + 1][0] = 0.0; fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - Gsq = coeff1*(static_cast(k*k) + ll); - fFFTin[lNFFT_2 + k][0] *= coeff2/(Gsq+coeff3); - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + Gsq = coeff1*(static_cast((k - NFFT)*(k - NFFT)) + ll); + fFFTin[lNFFT_2 + k][0] *= coeff2/(Gsq+coeff3); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } for (l = NFFT_2; l < NFFT; l += 2) { - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast((l-NFFT)*(l-NFFT)); for (k = 0; k < NFFT_2 - 1; k += 2) { Gsq = coeff1*(static_cast(k*k) + ll); @@ -921,42 +1051,53 @@ void TBulkTriVortexNGLFieldCalc::ManipulateFourierCoefficientsA() const { fFFTin[lNFFT_2 + k + 1][0] = 0.0; fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - Gsq = coeff1*(static_cast(k*k) + ll); - fFFTin[lNFFT_2 + k][0] *= coeff2/(Gsq+coeff3); - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + Gsq = coeff1*(static_cast((k-NFFT)*(k-NFFT)) + ll); + fFFTin[lNFFT_2 + k][0] *= coeff2/(Gsq+coeff3); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } //intermediate rows for (l = 1; l < NFFT_2; l += 2) { - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast(l*l); - for (k = 0; k < NFFT_2 - 1; k += 2) { + for (k = 0; k < NFFT_2; k += 2) { Gsq = coeff1*(static_cast((k+1)*(k+1)) + ll); fFFTin[lNFFT_2 + k][0] = 0.0; fFFTin[lNFFT_2 + k][1] = 0.0; fFFTin[lNFFT_2 + k + 1][0] *= coeff2/(Gsq+coeff3); fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - fFFTin[lNFFT_2 + k][0] = 0.0; - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + Gsq = coeff1*(static_cast((k+1-NFFT)*(k+1-NFFT)) + ll); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] *= coeff2/(Gsq+coeff3); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } for (l = NFFT_2 + 1; l < NFFT; l += 2) { - lNFFT_2 = l*(NFFT_2 + 1); + lNFFT_2 = l*(NFFT); ll = 3.0*static_cast((l-NFFT)*(l-NFFT)); - for (k = 0; k < NFFT_2 - 1; k += 2) { + for (k = 0; k < NFFT_2; k += 2) { Gsq = coeff1*(static_cast((k+1)*(k+1)) + ll); fFFTin[lNFFT_2 + k][0] = 0.0; fFFTin[lNFFT_2 + k][1] = 0.0; fFFTin[lNFFT_2 + k + 1][0] *= coeff2/(Gsq+coeff3); fFFTin[lNFFT_2 + k + 1][1] = 0.0; } - k = NFFT_2; - fFFTin[lNFFT_2 + k][0] = 0.0; - fFFTin[lNFFT_2 + k][1] = 0.0; + for (k = NFFT_2; k < NFFT; k += 2) { + Gsq = coeff1*(static_cast((k+1-NFFT)*(k+1-NFFT)) + ll); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] *= coeff2/(Gsq+coeff3); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } } fFFTin[0][0] = 0.0; @@ -1192,7 +1333,8 @@ void TBulkTriVortexNGLFieldCalc::ManipulateFourierCoefficientsForQy() const { void TBulkTriVortexNGLFieldCalc::CalculateSumAk() const { const int NFFT_2(fSteps/2); const int NFFTsq_2((fSteps/2 + 1)*fSteps); - + const int NFFTsq(fSteps*fSteps); +/* double SumFirstColumn(0.0); fSumAk = 0.0; @@ -1203,6 +1345,11 @@ void TBulkTriVortexNGLFieldCalc::CalculateSumAk() const { SumFirstColumn += fFFTin[l][0]; } fSumAk = 2.0*fSumAk - SumFirstColumn; +*/ + fSumAk = 0.0; + for (int l(0); l < NFFTsq; ++l) { + fSumAk += fFFTin[l][0]; + } return; } @@ -1248,7 +1395,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { FillAbrikosovCoefficients(); #pragma omp parallel for default(shared) private(l) schedule(dynamic) - for (l = 0; l < NFFT_2 + 1; l++) { + for (l = 0; l < NFFT; l++) { fCheckAkConvergence[l] = fFFTin[l][0]; } @@ -1262,7 +1409,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { #pragma omp parallel for default(shared) private(l) schedule(dynamic) for (l = 0; l < NFFTsq; l++) { - fOmegaMatrix[l] = fSumAk - fOmegaMatrix[l]; + fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0]; } // Calculate the gradient of omega @@ -1286,7 +1433,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { fQMatrix[l][0] = fQMatrixA[l][0]; fQMatrix[l][1] = fQMatrixA[l][1]; } - +/* fQMatrixA[0][0] = fQMatrixA[NFFT][0]; fQMatrixA[(NFFT+1)*NFFT_2][0] = fQMatrixA[0][0]; fQMatrix[0][0] = fQMatrixA[0][0]; @@ -1295,7 +1442,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { fQMatrixA[(NFFT+1)*NFFT_2][1] = fQMatrixA[0][1]; fQMatrix[0][1] = fQMatrixA[0][1]; fQMatrix[(NFFT+1)*NFFT_2][1] = fQMatrixA[0][1]; - +*/ // initialize B(x,y) with the mean field #pragma omp parallel for default(shared) private(l) schedule(dynamic) @@ -1306,6 +1453,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { bool akConverged(false), bkConverged(false), akInitiallyConverged(false), firstBkCalculation(true); double fourKappaSq(4.0*fKappa*fKappa); + while (!akConverged || !bkConverged) { // First iteration step for aK @@ -1313,17 +1461,18 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { #pragma omp parallel for default(shared) private(l) schedule(dynamic) for (l = 0; l < NFFTsq; l++) { if (fOmegaMatrix[l]) { - fOmegaMatrix[l] = fOmegaMatrix[l]*(fOmegaMatrix[l] + fQMatrix[l][0]*fQMatrix[l][0] + fQMatrix[l][1]*fQMatrix[l][1] - 2.0) + \ + fRealSpaceMatrix[l][0] = fOmegaMatrix[l]*(fOmegaMatrix[l] + fQMatrix[l][0]*fQMatrix[l][0] + fQMatrix[l][1]*fQMatrix[l][1] - 2.0) + \ (fOmegaDiffMatrix[l][0]*fOmegaDiffMatrix[l][0] + fOmegaDiffMatrix[l][1]*fOmegaDiffMatrix[l][1])/(fourKappaSq*fOmegaMatrix[l]); } else { - fOmegaMatrix[l] = 0.0; + fRealSpaceMatrix[l][0] = 0.0; } + fRealSpaceMatrix[l][1] = 0.0; } // At the two vortex cores g(r) is diverging; since all of this should be a smooth function anyway, I set the value of the next neighbour r // for the two vortex core positions in my matrix - fOmegaMatrix[0] = fOmegaMatrix[NFFT]; - fOmegaMatrix[(NFFT+1)*NFFT_2] = fOmegaMatrix[0]; + fRealSpaceMatrix[0][0] = fRealSpaceMatrix[NFFT][0]; + fRealSpaceMatrix[(NFFT+1)*NFFT_2][0] = fRealSpaceMatrix[0][0]; fftw_execute(fFFTplanOmegaToAk); @@ -1348,7 +1497,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { #pragma omp parallel for default(shared) private(l) schedule(dynamic) for (l = 0; l < NFFTsq; l++) { - fOmegaMatrix[l] = fSumAk - fOmegaMatrix[l]; + fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0]; } CalculateGradient(); @@ -1382,22 +1531,26 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { akConverged = true; - for (l = 0; l < NFFT_2 + 1; l++) { - if (((fabs(fFFTin[l][0]) > 1.0E-6) && (fabs(fCheckAkConvergence[l] - fFFTin[l][0])/fFFTin[l][0] > 1.0E-6)) || \ - (fCheckAkConvergence[l]/fFFTin[l][0] < 0.0)) { - // cout << "old: " << fCheckAkConvergence[l] << ", new: " << fFFTin[l][0] << endl; - akConverged = false; - break; + for (l = 0; l < NFFT; l++) { + if (fFFTin[l][0]){ + if (((fabs(fFFTin[l][0]) > 1.0E-6) && (fabs(fCheckAkConvergence[l] - fFFTin[l][0])/fFFTin[l][0] > 1.0E-6)) || \ + (fCheckAkConvergence[l]/fFFTin[l][0] < 0.0)) { + //cout << "old: " << fCheckAkConvergence[l] << ", new: " << fFFTin[l][0] << endl; + akConverged = false; + //cout << "index = " << l << endl; + break; + } } } if (!akConverged) { #pragma omp parallel for default(shared) private(l) schedule(dynamic) - for (l = 0; l < NFFT_2 + 1; l++) { + for (l = 0; l < NFFT; l++) { fCheckAkConvergence[l] = fFFTin[l][0]; } } else { akInitiallyConverged = true; + //break; } // cout << "Ak Convergence: " << akConverged << endl; @@ -1406,7 +1559,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { CalculateSumAk(); - // cout << "fSumAk = " << fSumAk << endl; + // cout << "fSumAk = " << fSumAk << " count = " << count << endl; // Do the Fourier transform to get omega(x,y) @@ -1414,13 +1567,13 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { #pragma omp parallel for default(shared) private(l) schedule(dynamic) for (l = 0; l < NFFTsq; l++) { - fOmegaMatrix[l] = fSumAk - fOmegaMatrix[l]; + fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0]; } CalculateGradient(); if (akInitiallyConverged) { // if the aK iterations converged, go on with the bK calculation - + //cout << "converged, count=" << count << endl; #pragma omp parallel for default(shared) private(l) schedule(dynamic) for (l = 0; l < NFFTsq; l++) { fBkMatrix[l][0] = fOmegaMatrix[l]*fFFTout[l] + fSumAk*(scaledB - fFFTout[l]) + \ @@ -1451,11 +1604,13 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { bkConverged = true; for (l = 0; l < NFFT; l++) { - if (((fabs(fBkMatrix[l][0]) > 1.0E-6) && (fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fabs(fBkMatrix[l][0]) > 1.0E-6)) || \ - (fCheckBkConvergence[l]/fBkMatrix[l][0] < 0.0)) { - // cout << "old: " << fCheckBkConvergence[l] << ", new: " << fBkMatrix[l][0] << endl; - bkConverged = false; - break; + if (fBkMatrix[l][0]) { + if (((fabs(fBkMatrix[l][0]) > 1.0E-6) && (fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fabs(fBkMatrix[l][0]) > 1.0E-6)) || \ + (fCheckBkConvergence[l]/fBkMatrix[l][0] < 0.0)) { + // cout << "old: " << fCheckBkConvergence[l] << ", new: " << fBkMatrix[l][0] << endl; + bkConverged = false; + break; + } } } diff --git a/src/external/TFitPofB-lib/classes/TFilmTriVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TFilmTriVortexFieldCalc.cpp index bf034bb1..28a31bd9 100644 --- a/src/external/TFitPofB-lib/classes/TFilmTriVortexFieldCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TFilmTriVortexFieldCalc.cpp @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2010/02/01 + 2010/02/21 ***************************************************************************/ @@ -834,12 +834,18 @@ void TFilmTriVortexNGLFieldCalc::CalculateGradient() const { fBkMatrix[l][1] = 0.0; } } -/* If the numerics is fine, this part is not needed +/* If the numerics is fine, this part is not needed */ // Ensure that omega and the gradient at the vortex-core positions are zero for (k = 0; k < NFFTz; ++k) { fOmegaMatrix[k] = 0.0; fOmegaMatrix[k + NFFTz*(NFFT+1)*NFFT_2] = 0.0; - + fOmegaDiffMatrix[0][k] = 0.0; + fOmegaDiffMatrix[0][k + NFFTz*(NFFT+1)*NFFT_2] = 0.0; + fOmegaDiffMatrix[1][k] = 0.0; + fOmegaDiffMatrix[1][k + NFFTz*(NFFT+1)*NFFT_2] = 0.0; + fOmegaDiffMatrix[2][k] = 0.0; + fOmegaDiffMatrix[2][k + NFFTz*(NFFT+1)*NFFT_2] = 0.0; + }/* for (i = 0; i < NFFT; ++i) { // j = 0 fOmegaDiffMatrix[0][k + NFFTz*NFFT*i] = 0.0; @@ -981,9 +987,10 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsA() const { const int NFFT(fSteps), NFFT_2(fSteps/2), NFFTz(fStepsZ), NFFTz_2(fStepsZ/2), NFFTsqStZ(fSteps*fSteps*fStepsZ), NFFTsq(fSteps*fSteps); // Divide EHB's coefficient no2 by two since we are considering "the full 3D reciprocal lattice", not only the half space! - const float symCorr(0.5f); + // Additionally treat all K the same (no difference between Kperp and K with Kz != 0) + const float symCorr(1.0f); const float coeff1(4.0f/3.0f*pow(PI/fLatticeConstant,2.0f)); - const float coeff3(4.0f*fKappa*fKappa); + const float coeff3(2.0f*fKappa*fKappa); const float coeff2(symCorr*coeff3/static_cast(NFFTsqStZ)); const float coeff4(4.0f*pow(PI/fThickness,2.0f)); @@ -1156,7 +1163,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsA() const { } fFFTin[k][0] = 0.0f; } - /* Comment out the negative Kz - setting them to zero instead to stay closer to EHB's solution */ + for (k = NFFTz_2; k < NFFTz; ++k) { kk = coeff4*static_cast((k - NFFTz)*(k - NFFTz)); for (i = 0; i < NFFT_2; i += 2) { @@ -1255,7 +1262,8 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsA() const { void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { const int NFFT(fSteps), NFFTsq(fSteps*fSteps), NFFT_2(fSteps/2), NFFTz(fStepsZ), NFFTz_2(fStepsZ/2); - // Divide EHB's PK and c by two since we are considering "the full 3D reciprocal lattice", not only the half space! + // Divide EHB's PK by two since we are considering "the full 3D reciprocal lattice", not only the half space! + // Additionally treat all K the same (no difference between Kperp and K with Kz != 0) const float coeffKsq(4.0f/3.0f*pow(PI/fLatticeConstant,2.0f)); const float coeffKy(TWOPI/fLatticeConstant); const float coeffKx(coeffKy/sqrt3); @@ -1277,7 +1285,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast(j*j) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1290,7 +1298,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast((j - NFFT)*(j - NFFT)) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1307,7 +1315,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast(j*j) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1320,7 +1328,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast((j - NFFT)*(j - NFFT)) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1341,7 +1349,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index2][1] = 0.0; } @@ -1354,7 +1362,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index2][1] = 0.0; } } @@ -1371,7 +1379,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index2][1] = 0.0; } @@ -1384,7 +1392,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(Gsq + 1.0f*fSumSum); fBkMatrix[index2][1] = 0.0; } } @@ -1407,7 +1415,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast(j*j) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1420,7 +1428,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast((j - NFFT)*(j - NFFT)) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1437,7 +1445,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast(j*j) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1450,7 +1458,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast((j - NFFT)*(j - NFFT)) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1471,7 +1479,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } @@ -1484,7 +1492,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } } @@ -1501,7 +1509,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } @@ -1514,7 +1522,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } } @@ -1534,7 +1542,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast(j*j) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1547,7 +1555,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast((j - NFFT)*(j - NFFT)) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1564,7 +1572,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast(j*j) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1577,7 +1585,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { Gsq = coeffKsq*(static_cast((j - NFFT)*(j - NFFT)) + ii); fBkMatrix[index][0] = \ (coeffPk*(ky*fQMatrix[index][1] + kx*fPkMatrix[index][1]) + \ - fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = 0.0; fBkMatrix[index2][1] = 0.0; @@ -1598,7 +1606,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } @@ -1611,7 +1619,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } } @@ -1628,7 +1636,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } @@ -1641,7 +1649,7 @@ void TFilmTriVortexNGLFieldCalc::ManipulateFourierCoefficientsB() const { fBkMatrix[index][1] = 0.0; fBkMatrix[index2][0] = \ (coeffPk*(ky*fQMatrix[index2][1] + kx*fPkMatrix[index2][1]) + \ - fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0*(Gsq + kk) + fSumSum); + 1.0f*fSumSum*fBkMatrix[index2][0] - sign*coeffBkS*sqrt(Gsq)*fBkS[j + 1 + NFFT*i][0])/(1.0f*(Gsq + kk) + fSumSum); fBkMatrix[index2][1] = 0.0; } } @@ -2541,24 +2549,7 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { // ... now fill in the Fourier components (Abrikosov) if everything was okay above - FillAbrikosovCoefficients(scaledB); - -/* try zeros - - for (k = 0; k < NFFTz_2; ++k) { - for (i = NFFT_2; i < NFFT; ++i) { - for (j = 0; j < NFFT; ++j) { - index = k + NFFTz*(j + NFFT*i); - fFFTin[index][0] = 0.0; - } - } - for (j = NFFT_2; j < NFFT; ++j) { - index = k + NFFTz*j; - fFFTin[index][0] = 0.0; - } - } - - end zeros */ + FillAbrikosovCoefficients(0.0); // save a few coefficients for the convergence check @@ -2598,7 +2589,6 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { CalculateGradient(); - // Calculate Q-Abrikosov float denomQAInv; @@ -2688,6 +2678,7 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { ManipulateFourierCoefficientsA(); + // Second iteration step for aK, first recalculate omega and its gradient CalculateSumAk(); @@ -2706,7 +2697,7 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { CalculateGradient(); -// CalculateGatVortexCore(); + //CalculateGatVortexCore(); // Get the spacial averages of the second iteration step for aK @@ -2723,8 +2714,8 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { fOmegaDiffMatrix[2][index]*fOmegaDiffMatrix[2][index])/(fourKappaSq*fOmegaMatrix[index]); } else { // cout << "! fOmegaMatrix at index " << index << endl; -// fSumSum -= fGstorage[k]; - index = k + fStepsZ*(j + fSteps*(i + 1)); + // fSumSum -= fGstorage[k]; + index = k + fStepsZ*(j + fSteps*(i + 1)); if (i < NFFT - 1 && fOmegaMatrix[index]) { fSumSum += fOmegaMatrix[index]*(1.0 - (fQMatrix[index][0]*fQMatrix[index][0] + fQMatrix[index][1]*fQMatrix[index][1])) - \ (fOmegaDiffMatrix[0][index]*fOmegaDiffMatrix[0][index] + \ @@ -2833,7 +2824,7 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { fBkS[index][1] = 0.f; } -// cout << "fC = " << fC << ", meanAk = " << meanAk << endl; + // cout << "fC = " << fC << ", meanAk = " << meanAk << endl; fSumSum = fC*meanAk; @@ -2890,7 +2881,7 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { } if (count == 50) { - cout << "3D iterations aborted after 50 steps" << endl; + cout << "3D iterations aborted after " << count << " steps" << endl; break; } @@ -2902,6 +2893,7 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { if (bkConverged && akConverged) { if (!fFind3dSolution) { //cout << "count = " << count << " 2D converged" << endl; + //cout << "2D iterations converged after " << count << " steps" << endl; //break; akConverged = false; bkConverged = false; @@ -3054,17 +3046,10 @@ void TFilmTriVortexNGLFieldCalc::CalculateGrid() const { fBkMatrix[l][1] = 0.0; } */ -/* - // If the iterations have converged, rescale the field from Brandt's units to Gauss - - #pragma omp parallel for default(shared) private(l) schedule(dynamic) - for (l = 0; l < NFFTsq; l++) { - fFFTout[l] *= Hc2_kappa; - } // Set the flag which shows that the calculation has been done fGridExists = true; return; -*/ + } diff --git a/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h index 0d5b1d51..08055544 100644 --- a/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h +++ b/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h @@ -148,6 +148,7 @@ private: mutable double *fOmegaMatrix; mutable fftw_complex *fOmegaDiffMatrix; mutable fftw_complex *fBkMatrix; + mutable fftw_complex *fRealSpaceMatrix; mutable fftw_complex *fQMatrix; mutable fftw_complex *fQMatrixA;