diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 40e7774..45c9e76 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -20,7 +20,7 @@ OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' -OptionsStruct.GradientDescentMethod = 'HeavyBall'; % 'HeavyBall' | 'NonLinearCGD' +OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; % 'HeavyBall' | 'NonLinearCGD' OptionsStruct.MaxIterationsForGD = 2E5; OptionsStruct.TimeStepSize = 1E-4; % in s OptionsStruct.MinimumTimeStepSize = 2E-10; % in s @@ -50,35 +50,6 @@ sim.Potential = pot.trap(); % + pot.repulsive_chopstick( Plotter.visualizeTrapPotential(sim.Potential,Params,Transf) %% - Plot initial wavefunction Plotter.visualizeWavefunction(psi,Params,Transf) -%% - Plot GS wavefunction -% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8'; -% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7'; -SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_8'; -JobNumber = 0; -Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) -%% -% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_SSD'; -% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Labyrinth'; -SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Honeycomb'; -JobNumber = 0; -Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) -%% - Plot GS wavefunction -SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio2_8'; -% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_7'; -% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_8'; -% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_9'; -JobNumber = 3; -Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) -%% -% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD'; -% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth'; -SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb'; -JobNumber = 0; -Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) -%% -SaveDirectory = './Results/Data_3D/GradientDescent'; -JobNumber = 0; -Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% % To reproduce results from the Blair Blakie paper: @@ -438,4 +409,40 @@ JobNumber = 2; % 82 % JobNumber = 3; % 83 Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber); - +%% - Plot GS wavefunction +% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8'; +% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7'; +SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_8'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_SSD'; +% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Labyrinth'; +SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Honeycomb'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% - Plot GS wavefunction +SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio2_8'; +% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_7'; +% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_8'; +% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_9'; +JobNumber = 3; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD'; +% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth'; +SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +SaveDirectory = './Results/Data_3D/TiltedDipoles15'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +SaveDirectory = './Results/Data_3D/TiltedDipoles30'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +SaveDirectory = './Results/Data_3D/TiltedDipoles45'; +JobNumber = 1; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 5280b0e..1608c1e 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -1,18 +1,20 @@ -%% Scaled parameters +theta = 0; +phi = 0; -ScalingFactor = (4/5)^2; +% - SSD: N = 1E5, as = 86ao OptionsStruct = struct; -OptionsStruct.NumberOfAtoms = sqrt(ScalingFactor) * 5E5; -OptionsStruct.DipolarPolarAngle = deg2rad(0); -OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 75; +OptionsStruct.NumberOfAtoms = 1E5; +OptionsStruct.DipolarPolarAngle = deg2rad(theta); +OptionsStruct.DipolarAzimuthAngle = deg2rad(phi); +OptionsStruct.ScatteringLength = 86; -AspectRatio = 2.8; -HorizontalTrapFrequency = 125/ScalingFactor; -VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; -OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; +% AspectRatio = 2.0; +% HorizontalTrapFrequency = 125; +% VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; +% OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; +OptionsStruct.TrapFrequencies = [150, 150, 300]; OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.NumberOfGridPoints = [128, 128, 64]; @@ -21,18 +23,18 @@ OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 1E-4; % in s -OptionsStruct.MinimumTimeStepSize = 2E-6; % in s +OptionsStruct.TimeStepSize = 1E-3; % in s +OptionsStruct.MinimumTimeStepSize = 2E-6; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; -OptionsStruct.PlotLive = false; -OptionsStruct.JobNumber = 2; -OptionsStruct.RunOnGPU = true; +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = false; OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/ApproximateLHY/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_')); +OptionsStruct.SaveDirectory = './Results/Data_3D/TiltedDipoles0'; options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct @@ -43,46 +45,3 @@ sim.Potential = pot.trap(); %-% Run Simulation %-% [Params, Transf, psi, V, VDk] = sim.run(); - -%{ -%% - Gradient Descent Test - -OptionsStruct = struct; - -OptionsStruct.NumberOfAtoms = 8E4; -OptionsStruct.DipolarPolarAngle = deg2rad(0); -OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 95; - -OptionsStruct.TrapFrequencies = [30, 60, 90]; -OptionsStruct.TrapPotentialType = 'Harmonic'; - -OptionsStruct.NumberOfGridPoints = [256, 128, 128]; -OptionsStruct.Dimensions = [30, 20, 20]; -OptionsStruct.UseApproximationForLHY = true; -OptionsStruct.IncludeDDICutOff = true; -OptionsStruct.CutoffType = 'Cylindrical'; -OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' -OptionsStruct.MaxIterationsForGD = 2E4; -% OptionsStruct.TimeStepSize = 1E-4; % in s -% OptionsStruct.MinimumTimeStepSize = 2E-10; % in s -% OptionsStruct.TimeCutOff = 2E6; % in s -% OptionsStruct.EnergyTolerance = 5E-10; -% OptionsStruct.ResidualTolerance = 1E-08; -OptionsStruct.NoiseScaleFactor = 0.01; - -OptionsStruct.PlotLive = true; -OptionsStruct.JobNumber = 0; -OptionsStruct.RunOnGPU = false; -OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_3D/GradientDescent'; -options = Helper.convertstruct2cell(OptionsStruct); -clear OptionsStruct - -sim = Simulator.DipolarGas(options{:}); -pot = Simulator.Potentials(options{:}); -sim.Potential = pot.trap(); % + pot.repulsive_chopstick(); - -%-% Run Simulation %-% -[Params, Transf, psi, V, VDk] = sim.run(); -%} \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m index 898d19b..0d16c30 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m @@ -41,8 +41,19 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad) end % == Setting up the initial wavefunction == % - psi = this.setupWavefunction(Params,Transf); + WavefunctionFile = fullfile(this.SaveDirectory, 'psi_init.mat'); + if isfile(WavefunctionFile) + loadedData = load(WavefunctionFile); + if isgpuarray(loadedData.psi) + psi = gather(loadedData.psi); + else + psi = loadedData.psi; + end + else + psi = this.setupWavefunction(Params,Transf); + end + if this.RunOnGPU psi = gpuArray(psi); end diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m index 8ac935c..727c48a 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m @@ -12,6 +12,7 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) Observ.res = 1; psi_old = psi; % Previous psi value (for heavy-ball method) + % Live Plotter if this.PlotLive Plotter.plotLive(psi,Params,Transf,Observ) drawnow @@ -24,7 +25,6 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) J = compute_gradient(psi, Params, Transf, VDk, V); % Calculate chemical potential and norm - % Can also be calculated as --> muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); muchem = sum(real(conj(psi(:)) .* J(:))) / sum(abs(psi(:)).^2); % Calculate residual and check convergence @@ -83,12 +83,11 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) disp('Saving data...'); save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V'); disp('Save complete!'); - case 'NonLinearCGD' - % Define the function handle - f = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N; + case 'NonLinearCGD' + % Convergence Criteria: - Epsilon = 1E-5; + epsilon = 1E-13; % Iteration Counter: i = 1; @@ -98,65 +97,80 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) % Initialize the PrematureExitFlag to false PrematureExitFlag = false; + % Live plotter if this.PlotLive Plotter.plotLive(psi,Params,Transf,Observ) drawnow end - + % Minimization Loop while true % Compute gradient - J = compute_gradient(psi, Params, Transf, VDk, V); - - % Check stopping criterion (Gradient norm) - if norm(J(:)) < Epsilon - disp('Tolerance reached: Gradient norm is below the specified epsilon.'); + J = compute_gradient(psi, Params, Transf, VDk, V); + % Calculate chemical potential + muchem = real(conj(psi(:))' * J(:)) / norm(psi(:))^2; + % Calculate residual + residual = J - (muchem * psi); + + % Compute energy difference between the last two saved energy values + if i == 1 + energydifference = NaN; + elseif mod(i,100) == 0 && length(Observ.EVec) > 1 + energydifference = abs(Observ.EVec(end) - Observ.EVec(end-1)); + end + % Convergence check - if energy difference is below set tolerance, then exit + if energydifference <= epsilon + disp('Tolerance reached: Energy difference is below the specified epsilon.'); PrematureExitFlag = true; % Set flag to indicate premature exit break; - elseif i >= this.MaxIterationsForCGD + elseif i >= this.MaxIterationsForGD % If set maximum number of iterations reached, then exit disp('Maximum number of iterations for CGD reached.'); PrematureExitFlag = true; % Set flag to indicate premature exit break; end - + % Initialize search direction if first iteration if i == 1 - S = -J; - else - % Update search direction - S = update_search_direction(S, J, J_old); + d = -residual; + else % Compute beta via Polak-Ribiere and create new direction + residual_new = residual; + beta = compute_beta(residual_new, residual_old); + d = -residual_new + beta * p_old; end - - % Step Size Optimization (Line Search) - alpha = optimize_step_size(f, psi, S, Params, Transf, VDk, V); - + + residual_old = residual; + p = d - (real(conj(d(:))' * psi(:)) .* psi); + p_old = p; + + % Compute optimal theta to generate psi in the direction of minimum in the energy landscape + theta = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V); + % Update solution - psi = psi + alpha * S; - + psi = (cos(theta).*psi) + (sin(theta).*(p / norm(p(:)))); % Normalize psi - Norm = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz; - psi = sqrt(Params.N) * psi / sqrt(Norm); - - % Store old gradient - J_old = J; - i = i + 1; - muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); - - if mod(i,500) == 0 + Norm = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz; + psi = sqrt(Params.N) * psi / sqrt(Norm); + i = i + 1; + + % Calculate chemical potential with new psi + muchem = real(conj(psi(:))' * J(:)) / norm(psi(:))^2; + + if mod(i,100) == 0 - % Change in Energy + % Collect Energy value E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); E = E/Norm; Observ.EVec = [Observ.EVec E]; - % Chemical potential + % Collect Chemical potential value Observ.mucVec = [Observ.mucVec muchem]; - % Normalized residuals + % Collect Normalized residuals res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem); Observ.residual = [Observ.residual res]; Observ.res_idx = Observ.res_idx + 1; + % Live plotter if this.PlotLive Plotter.plotLive(psi,Params,Transf,Observ) drawnow @@ -185,7 +199,7 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) end %% Modules -% Numerical Gradient Calculation using the finite differences method + function J = compute_gradient(psi, Params, Transf, VDk, V) % Operators @@ -210,44 +224,48 @@ function J = compute_gradient(psi, Params, Transf, VDk, V) H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w); - J = H(psi); + J = H(psi); end -% Backtracking Line Search (Step Size Optimization) -function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V) - alpha = 1; % Initial step size - rho = 0.5; % Step size reduction factor - c = 1E-4; % Armijo condition constant - max_iter = 100; % Max iterations for backtracking - tol = 1E-4; % Tolerance for stopping +function g = compute_g(psi, p, Params, VDk) + + rho = real(psi.*conj(p)); - grad = compute_gradient(X, Params, Transf, VDk, V); % Compute gradient once - f_X = f(X); % Evaluate f(X) once + % Contact interactions + C = Params.gs*Params.N; + gint = @(w)(C.*rho).*w; + + % DDIs + D = Params.gdd*Params.N; + rhotilde = fftn(rho); + Phi = real(ifftn(rhotilde.*VDk)); + gaddi = @(w)(D.*Phi).*w; - for k = 1:max_iter - % Evaluate Armijo condition with precomputed f(X) and grad - if f(X + alpha * S) <= f_X + c * alpha * (S(:)' * grad(:)) - break; - else - alpha = rho * alpha; % Reduce the step size - end - - % Early stopping if step size becomes too small - if alpha < tol - break; - end - end + % Quantum fluctuations + eps_dd = Params.add/Params.as; + Q = (4/(3*pi^2)) * (C^(5/2)/Params.N) * (1 + ((3*eps_dd^2)/2)); + gqf = @(w)(((3/2)*Q).*abs(psi).*rho).*w; + + gop = @(w) gint(w) + gaddi(w) + gqf(w); + + g = gop(psi); end -% Update Search Direction -function S_new = update_search_direction(S, J_new, J_old) - % (Fletcher-Reeves method) - % beta = (norm(J_new(:))^2) / (norm(J_old(:))^2); - % S_new = -J_new + beta * S; +% Optimal direction via line search +function theta = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V) + + Hpsi = compute_gradient(psi, Params, Transf, VDk, V); + Hp = compute_gradient(p, Params, Transf, VDk, V); + g = compute_g(psi, p, Params, VDk); + numerator = real(conj(p(:))' * Hpsi(:))/norm(p(:)); + denominator = muchem - ((conj(p(:))' * Hp(:)) + real(conj(g(:))' * p(:)))/norm(p(:))^2; - % (Polak-Ribiere method) - beta = max(0, (J_new(:)' * (J_new(:) - J_old(:))) / (norm(J_old(:))^2)); - S_new = -J_new + beta * S; + theta = numerator / denominator; +end + +% Optimal step size via Polak-Ribiere +function beta = compute_beta(residual_new, residual_old) + beta = max(0, (residual_new(:)' * (residual_new(:) - residual_old(:))) / (norm(residual_old(:))^2)); end \ No newline at end of file