From abe402313e290a7c509e93676da498681d6a7833 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 31 Mar 2025 14:09:32 +0200 Subject: [PATCH] Latest scripts for testing of energy minimization via CGD. --- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 1 + .../+Scripts/run_on_cluster_beyondSSD.m | 8 +- .../+Simulator/@DipolarGas/DipolarGas.m | 4 + .../@DipolarGas/minimizeEnergyFunctional.m | 109 ++++++++++++------ 4 files changed, 80 insertions(+), 42 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index b499fcd..7c8a8ed 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -23,6 +23,7 @@ OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.MaxIterationsForCGD = 10; OptionsStruct.TimeStepSize = 1E-4; % in s OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m index c5438c7..862181b 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m @@ -45,10 +45,10 @@ sim.Potential = pot.trap(); OptionsStruct = struct; -OptionsStruct.NumberOfAtoms = 4E5; +OptionsStruct.NumberOfAtoms = 8E5; OptionsStruct.DipolarPolarAngle = deg2rad(0); OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 87; +OptionsStruct.ScatteringLength = 85; AspectRatio = 2.0; HorizontalTrapFrequency = 125; @@ -63,7 +63,7 @@ OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.TimeStepSize = 1E-4; % in s -OptionsStruct.MinimumTimeStepSize = 2E-10; % in s +OptionsStruct.MinimumTimeStepSize = 1E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-08; @@ -73,7 +73,7 @@ OptionsStruct.PlotLive = false; OptionsStruct.JobNumber = 0; OptionsStruct.RunOnGPU = true; OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_3D/BeyondSSD_Labyrinth'; +OptionsStruct.SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth'; options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m index c30b382..1214fbe 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -17,6 +17,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable EnergyTolerance; ResidualTolerance; NoiseScaleFactor; + MaxIterationsForCGD; Calculator; @@ -69,6 +70,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'NoiseScaleFactor', 4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'MaxIterationsForCGD', 100,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'IncludeDDICutOff', true,... @islogical); @@ -104,6 +107,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.EnergyTolerance = p.Results.EnergyTolerance; this.ResidualTolerance = p.Results.ResidualTolerance; this.NoiseScaleFactor = p.Results.NoiseScaleFactor; + this.MaxIterationsForCGD = p.Results.MaxIterationsForCGD; this.IncludeDDICutOff = p.Results.IncludeDDICutOff; this.UseApproximationForLHY = p.Results.UseApproximationForLHY; diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m index d4cfd63..069389e 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m @@ -1,4 +1,4 @@ -function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V) +function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V) format long; @@ -8,35 +8,58 @@ function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V) % Convergence Criteria: Epsilon = 1E-5; - % Gradient Calculation - J = compute_gradient(psi, Params, Transf, VDk, V); % Gradient at initial point - S = -J; % Initial search direction - % Iteration Counter: i = 1; - % Minimization - while norm(S(:)) > Epsilon % Halting Criterion + % Initialize the PrematureExitFlag to false + PrematureExitFlag = false; + + % 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.'); + PrematureExitFlag = true; % Set flag to indicate premature exit + break; + elseif i >= this.MaxIterationsForCGD + 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); + end + % Step Size Optimization (Line Search) - alpha = optimize_step_size(f, psi, S, Params, Transf, VDk, V); - - % Update the solution - x_o_new = psi + alpha * S; - - % Update the gradient and search direction + alpha = optimize_step_size(f, psi, S, Params, Transf, VDk, V); + + % Update solution + psi = psi + alpha * S; + + % 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; - J = compute_gradient(x_o_new, Params, Transf, VDk, V); - S = update_search_direction(S, J, J_old); - - % Update for next iteration - psi = x_o_new; - Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; - psi = sqrt(Params.N)*psi/sqrt(Norm); i = i + 1; end - disp('Minimum found!'); - fprintf('Number of Iterations for Convergence: %d\n\n', i); + % Check if loop ended prematurely + if PrematureExitFlag + disp('Optimizer ended prematurely without convergence to a minimum.'); + else + disp('Minimum found!'); + fprintf('Number of Iterations for Convergence: %d\n\n', i); + end muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); @@ -48,26 +71,31 @@ end %% Helper functions % Numerical Gradient Calculation using the finite differences method -function J = compute_gradient(psi, Params, Transf, VDk, V) +function J = compute_gradient(psi, Params, Transf, VDk, V) - % Kinetic energy - KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); - H_o = KEop + V; + % Operators - %Contact interactions - Hint = Params.gs*abs(psi).^2; + % Kinetic energy + KEop = 0.5 * (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + HKin = @(w) real(ifft(KEop.*fft(w))); + + % Trap Potential + HV = @(w) V.*w; + + % Contact interactions + Hint = @(w) (Params.gs*abs(psi).^2).*w; % DDIs frho = fftn(abs(psi).^2); Phi = real(ifftn(frho.*VDk)); - Hddi = Params.gdd*Phi; - - %Quantum fluctuations - Hqf = Params.gammaQF*abs(psi).^3; - - H = H_o + Hint + Hddi + Hqf; + Hddi = @(w) (Params.gdd*Phi).*w; - J = H.*psi; + % Quantum fluctuations + Hqf = @(w) (Params.gammaQF*abs(psi).^3).*w; + + H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w); + + J = H(psi); end @@ -81,7 +109,7 @@ function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V) grad = compute_gradient(X, Params, Transf, VDk, V); % Compute gradient once f_X = f(X); % Evaluate f(X) once - + 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(:)) @@ -98,8 +126,13 @@ function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V) end -% Update Search Direction (Polak-Ribiere method) +% Update Search Direction function S_new = update_search_direction(S, J_new, J_old) - beta = (norm(J_new(:))^2) / (norm(J_old(:))^2); + % (Fletcher-Reeves method) + % beta = (norm(J_new(:))^2) / (norm(J_old(:))^2); + % S_new = -J_new + beta * S; + + % (Polak-Ribiere method) + beta = max(0, (J_new(:)' * (J_new(:) - J_old(:))) / (norm(J_old(:))^2)); S_new = -J_new + beta * S; end \ No newline at end of file