karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Solver runs till BdG equations are solved and results saved.

Browse Source
This commit is contained in:
Karthik 2024-11-23 22:53:58 +01:00
parent 5b1428e91a
commit 5049970246
3 changed files with 81 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

96
Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m
View File

@ -1,41 +1,67 @@
function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem)
% 2-D matrix will be unravelled to a single column vector and the corresponding BdG matrix of (N^2)^2 elements solved for.
Size = length(psi(:));
Neigs = length(psi(:));
opts.tol = 1e-16;
opts.disp = 1;
opts.issym = 0;
opts.isreal = 1;
opts.maxit = 1e4;
BdGVec = @(g) BdGSolver2D.BdGMatrix(g, psi, Params, VDk, VParams, Transf, muchem); % This function takes a column vector as input and returns a
% matrix-vector product which is also a column vector
[g,D] = eigs(BdGVec,Size,Neigs,'sr',opts);
evals = diag(D);
clear D;
gs = Params.gs;
gdd = Params.gdd;
gammaQF = Params.gammaQF;
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2);
g_pf_2D = 1/(sqrt(2*pi)*VParams.ell);
gQF_pf_2D = sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2));
Ez = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
muchem_tilde = muchem - Ez;
% eigs only works with column vectors
psi = psi.';
KEop = KEop.';
VDk = VDk.';
% Interaction Potential
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
% Operators
H = @(w) real(ifft(KEop.*fft(w)));
C = @(w) (((g_pf_2D*gs*abs(psi).^2) + (g_pf_2D*gdd*Phi)).*w) + (gQF_pf_2D*gammaQF*abs(psi).^3.*w);
muHC = @(w) (-muchem_tilde * w) + H(w) + C(w);
X = @(w) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w;
% Eigenvalues
evals = sqrt(evals);
% Obtain f from g
f = zeros(size(g));
for ii = 1:Neigs
f(:,:,ii) = (1/evals(ii)) * (muHC(g(:,:,ii)) + (2.*X(g(:,:,ii))));
end
% Obtain u and v from f and g
u = (f + g)/2;
v = (f - g)/2;
% Renormalize to \int |u|^2 - |v|^2 = 1
for ii=1:Neigs
normalization = sum(abs(u(:,:,ii)).^2 - abs(v(:,:,ii)).^2);
u(:,:,ii) = u(:,:,ii)/sqrt(normalization);
v(:,:,ii) = v(:,:,ii)/sqrt(normalization);
end
modes.u = u'; modes.v = v';
modes.g = g'; modes.f = f';
% 2-D matrix will be unravelled to a single column vector and the corresponding BdG matrix of (N^2)^2 elements solved for.
Size = length(psi(:));
Neigs = length(psi(:));
opts.tol = 1e-16;
opts.disp = 1;
opts.issym = 0;
opts.isreal = 1;
opts.maxit = 1e4;
BdGVec = @(g) BdGSolver2D.BdGMatrix(g, psi, Params, VDk, VParams, Transf, muchem); % This function takes a column vector as input and returns a
% matrix-vector product which is also a column vector
[g,D] = eigs(BdGVec,Size,Neigs,'sr',opts);
evals = diag(D);
clear D;
% Eigenvalues
evals = sqrt(evals);
% Obtain f from g
for ii = 1:Neigs
gres = reshape(g(:,ii), size(psi));
f(:,ii) = reshape((1/evals(ii)) * (muHC(gres) + (2.*X(gres))), [], 1);
end
% Obtain u and v from f and g
u = (f + g)/2;
v = (f - g)/2;
% Renormalize to \int |u|^2 - |v|^2 = 1
for ii=1:Neigs
normalization = sum(abs(u(:,ii)).^2 - abs(v(:,ii)).^2);
u(:,ii) = u(:,ii)/sqrt(normalization);
v(:,ii) = v(:,ii)/sqrt(normalization);
end
modes.u = u'; modes.v = v';
modes.g = g'; modes.f = f';
end

21
Dipolar-Gas-Simulator/+Scripts/run_locally.m
View File

@ -100,36 +100,33 @@ pot = VariationalSolver2D.Potentials(options{:
solver.Potential = pot.trap();
%-% Run Solver %-%
% [Params, Transf, psi, V, VDk] = solver.run();
[Params, Transf, psi, V, VDk] = solver.run();
% Solve BdG equations
% Load data
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Params','Transf');
Params = Data.Params;
Transf = Data.Transf;
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V');
Transf = Data.Transf;
Params = Data.Params;
VParams = Data.VParams;
V = Data.V;
if isgpuarray(Data.psi)
psi = gather(Data.psi);
else
psi = Data.psi;
end
VParams.ell = Params.ell;
% == DDI potential == %
VDk = solver.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);
% == Trap potential == %
X = Transf.X; Y = Transf.Y;
V = 0.0*(Params.gx.*X.^2+Params.gy.*Y.^2);
% == Chemical potential == %
muchem = solver.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
[evals, modes] = BdGSolver2D.solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem);
% Save the eigenvalues and eigenvectors to a .mat file
save(sprintf(strcat(solver.SaveDirectory, '/Run_%03i/bdg_eigen_data.mat'),solver.JobNumber), 'evals', 'modes');
save(sprintf(strcat(solver.SaveDirectory, '/Run_%03i/bdg_eigen_data.mat'),solver.JobNumber), 'evals', 'modes', '-v7.3');
%% - Create Variational2D and Calculator object with specified options

25
Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
View File

@ -164,9 +164,9 @@ OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = true;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = false;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data_TriangularPhase';
options = Helper.convertstruct2cell(OptionsStruct);
@ -177,33 +177,30 @@ pot = VariationalSolver2D.Potentials(options{:
solver.Potential = pot.trap();
%-% Run Solver %-%
% [Params, Transf, psi, V, VDk] = solver.run();
[Params, Transf, psi, V, VDk] = solver.run();
% Solve BdG equations
% Load data
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Params','Transf');
Params = Data.Params;
Transf = Data.Transf;
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V');
Transf = Data.Transf;
Params = Data.Params;
VParams = Data.VParams;
V = Data.V;
if isgpuarray(Data.psi)
psi = gather(Data.psi);
else
psi = Data.psi;
end
VParams.ell = Params.ell;
% == DDI potential == %
VDk = solver.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);
% == Trap potential == %
X = Transf.X; Y = Transf.Y;
V = 0.0*(Params.gx.*X.^2+Params.gy.*Y.^2);
% == Chemical potential == %
muchem = solver.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
[evals, modes] = BdGSolver2D.solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem);
% Save the eigenvalues and eigenvectors to a .mat file
save(sprintf(strcat(solver.SaveDirectory, '/Run_%03i/bdg_eigen_data.mat'),Params.njob), 'evals', 'modes');
save(sprintf(strcat(solver.SaveDirectory, '/Run_%03i/bdg_eigen_data.mat'),solver.JobNumber), 'evals', 'modes', '-v7.3');