diff --git a/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m b/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m index 59df0a6..82ca429 100644 --- a/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m +++ b/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m @@ -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 diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 683d5d0..2bc0f81 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -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 diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index afc632c..227399a 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -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');