From 196663e367b80b1286039d9ba99aca73580eff2e Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 22 Nov 2024 16:18:27 +0100 Subject: [PATCH] Added and tested 2D BdG solver - works until Matlab starts solving for the eigen-values and vectors. --- .../+BdGSolver2D/BdGMatrix.m | 35 ++++++++++ .../solveBogoliubovdeGennesIn2D.m | 38 ++-------- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 37 ++++++++-- .../+Scripts/run_on_cluster.m | 69 +++++++++++++++++++ 4 files changed, 141 insertions(+), 38 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+BdGSolver2D/BdGMatrix.m diff --git a/Dipolar-Gas-Simulator/+BdGSolver2D/BdGMatrix.m b/Dipolar-Gas-Simulator/+BdGSolver2D/BdGMatrix.m new file mode 100644 index 0000000..1c7bb80 --- /dev/null +++ b/Dipolar-Gas-Simulator/+BdGSolver2D/BdGMatrix.m @@ -0,0 +1,35 @@ +function BdGVec = BdGMatrix(gVec, psi, Params, VDk, VParams, Transf, muchem) + + 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; + + g = reshape(gVec, size(psi)); + + % Operate in order on g + BdG = muHC(muHC(g) + (2.*X(g))); + + BdGVec = reshape(BdG, [], 1); + +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m b/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m index 44a1430..404d771 100644 --- a/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m +++ b/Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m @@ -1,44 +1,16 @@ function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem) -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) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w; - -% Operate in order on g -BdG = @(g) muHC(muHC(g) + (2.*X(g))); - -syssize = size(psi); -opts.v0 = psi(:); +Size = length(psi(:)); +Neigs = length(psi(:)); opts.tol = 1e-16; opts.disp = 1; opts.issym = 0; opts.isreal = 1; opts.maxit = 1e4; -Neigs = syssize; -[g,D] = eigs(BdG,syssize,Neigs,'sr',opts); +BdGVec = @(g) BdGSolver2D.BdGMatrix(g, psi, Params, VDk, VParams, Transf, muchem); + +[g,D] = eigs(BdGVec,Size,Neigs,'sr',opts); evals = diag(D); clear D; diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index d6de98e..4f08fb3 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -88,7 +88,7 @@ OptionsStruct.WidthUpperBound = 12; OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.PlotLive = true; -OptionsStruct.JobNumber = 2; +OptionsStruct.JobNumber = 1; OptionsStruct.RunOnGPU = false; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Data_TriangularPhase'; @@ -100,7 +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; + +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); %% - Create Variational2D and Calculator object with specified options @@ -130,7 +156,7 @@ OptionsStruct.WidthUpperBound = 12; OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.PlotLive = true; -OptionsStruct.JobNumber = 2; +OptionsStruct.JobNumber = 1; OptionsStruct.RunOnGPU = false; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Data_StripePhase'; @@ -172,7 +198,7 @@ OptionsStruct.WidthUpperBound = 12; OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.PlotLive = true; -OptionsStruct.JobNumber = 2; +OptionsStruct.JobNumber = 1; OptionsStruct.RunOnGPU = false; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Data_HoneycombPhase'; @@ -211,4 +237,5 @@ Scripts.analyzeRun2D(SaveDirectory, JobNumber) SaveDirectory = './Data_HoneycombPhase'; JobNumber = 3; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) -Scripts.analyzeRun2D(SaveDirectory, JobNumber) \ No newline at end of file +Scripts.analyzeRun2D(SaveDirectory, JobNumber) + diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index b257e2a..ff0cf5e 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -10,6 +10,7 @@ % as = ((as/add)*Params.add)/Params.a0 % Critical point: 102.5133; Triangular phase: 98.0676; Stripe phase: 100.0289; Honeycomb phase: 101.9903 +%{ %% - Create Variational2D and Calculator object with specified options OptionsStruct = struct; @@ -135,3 +136,71 @@ solver.Potential = pot.trap(); %-% Run Solver %-% [Params, Transf, psi, V, VDk] = solver.run(); +%} +%% - Create Variational2D and Calculator object with specified options + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 2.0030e+07; +OptionsStruct.DipolarPolarAngle = 0; +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 98.0676; + +OptionsStruct.TrapFrequencies = [10, 10, 72.4]; +OptionsStruct.TrapPotentialType = 'None'; + +OptionsStruct.NumberOfGridPoints = [128, 128]; +OptionsStruct.Dimensions = [100, 100]; +OptionsStruct.TimeStepSize = 500E-6; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E6; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-04; +OptionsStruct.NoiseScaleFactor = 4; + +OptionsStruct.MaxIterations = 1; +OptionsStruct.VariationalWidth = 5.0; +OptionsStruct.WidthLowerBound = 0.2; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 1e-2; + +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 1; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Data_TriangularPhase'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + +%-% Run Solver %-% +% [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; + +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);