Temporary upload in order to fix an error in my calculations - some parts are not working at all
This commit is contained in:
Bastian M. Wojek
parent
43ec66ad7c
commit
a63868369a
@ -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<double>(NFFT));
|
||||
const double denomY(2.0/static_cast<double>(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<double>(k*k) + 3.0*static_cast<double>(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<double>(k*k) + 3.0*static_cast<double>(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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>(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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>(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<double>(k*k) + 3.0*static_cast<double>((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<double>(k*k) + 3.0*static_cast<double>((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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((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<double>((k+1)*(k+1)) + 3.0*static_cast<double>(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<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>(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<double>((k+1)*(k+1)) + 3.0*static_cast<double>((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<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>((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<double>(l)/(static_cast<double>(k*k) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(l)/(static_cast<double>(k*k) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(l)/(static_cast<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(l)/(static_cast<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>(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<double>(l-NFFT)/(static_cast<double>(k*k) + 3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(l-NFFT)/(static_cast<double>(k*k) + 3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(l-NFFT)/(static_cast<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(l-NFFT)/(static_cast<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((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<double>(l)/(static_cast<double>((k+1)*(k+1)) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k + 1][0] *= coeff1*static_cast<double>(l)/(static_cast<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>(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<double>(l-NFFT)/(static_cast<double>((k+1)*(k+1)) + 3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k + 1][0] *= coeff1*static_cast<double>(l-NFFT)/(static_cast<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>((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<double>(k)/(static_cast<double>(k*k) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(k)/(static_cast<double>(k*k) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(k-NFFT)/(static_cast<double>((k-NFFT)*(k-NFFT))+3.0*static_cast<double>(l*l));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(k-NFFT)/(static_cast<double>((k-NFFT)*(k-NFFT))+3.0*static_cast<double>(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<double>(k)/(static_cast<double>(k*k)+3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(k)/(static_cast<double>(k*k)+3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(k-NFFT)/(static_cast<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
if (NFFT_2 % 2) {
|
||||
F[lNFFT + k][0] *= coeff1*static_cast<double>(k-NFFT)/(static_cast<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((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<double>(k+1)/(static_cast<double>((k+1)*(k+1)) + 3.0*static_cast<double>(l*l));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k + 1][0] *= coeff1*static_cast<double>(k+1-NFFT)/(static_cast<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>(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<double>(k+1)/(static_cast<double>((k+1)*(k+1)) + 3.0*static_cast<double>((l-NFFT)*(l-NFFT)));
|
||||
}
|
||||
|
||||
for (k = NFFT_2; k < NFFT - 1; k += 2) {
|
||||
F[lNFFT + k + 1][0] *= coeff1*static_cast<double>(k+1-NFFT)/(static_cast<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>((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<double>(k*k) + 3.0*static_cast<double>(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<double>(k*k) + 3.0*static_cast<double>(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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>(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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>(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<double>(k*k) + 3.0*static_cast<double>((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<double>(k*k) + 3.0*static_cast<double>((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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((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<double>((k-NFFT)*(k-NFFT)) + 3.0*static_cast<double>((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<double>((k + 1)*(k + 1)) + 3.0*static_cast<double>(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<double>((k + 1-NFFT)*(k + 1-NFFT)) + 3.0*static_cast<double>(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<double>((k+1)*(k+1)) + 3.0*static_cast<double>((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<double>((k+1-NFFT)*(k+1-NFFT)) + 3.0*static_cast<double>((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<double>(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<double>(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;
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -53,8 +53,8 @@ public:
|
||||
virtual double* DataB() const {return fFFTout;}
|
||||
virtual void SetParameters(const vector<double>& 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;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -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)
|
||||
};
|
||||
|
||||
@ -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"
|
||||
</comment>
|
||||
<debug>0</debug>
|
||||
<wisdom>/home/l_wojek/analysis/WordsOfWisdom.dat</wisdom>
|
||||
<delta_t>0.01</delta_t>
|
||||
<delta_B>0.1</delta_B>
|
||||
<VortexLattice>
|
||||
<N_VortexGrid>256</N_VortexGrid>
|
||||
<VortexSymmetry>triangular</VortexSymmetry>
|
||||
</VortexLattice>
|
||||
<LEM>
|
||||
<data_path>/home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000-</data_path>
|
||||
|
||||
10
src/external/TFitPofB-lib/test/testVortex-v2.cpp
vendored
10
src/external/TFitPofB-lib/test/testVortex-v2.cpp
vendored
@ -9,7 +9,7 @@ using namespace std;
|
||||
|
||||
int main(){
|
||||
|
||||
unsigned int NFFT(512);
|
||||
unsigned int NFFT(256);
|
||||
|
||||
vector<double> parForVortex;
|
||||
parForVortex.resize(3);
|
||||
@ -20,17 +20,17 @@ int main(){
|
||||
|
||||
vector<double> parForPofB;
|
||||
parForPofB.push_back(0.01); //dt
|
||||
parForPofB.push_back(2.0); //dB
|
||||
parForPofB.push_back(0.1); //dB
|
||||
|
||||
vector<double> 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);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user