From b47d81e68d5ac12b71fb173fdf7ea86680afdf3a Mon Sep 17 00:00:00 2001 From: Karthik Date: Mon, 28 Apr 2025 11:46:54 +0200 Subject: [PATCH] Major corrections to gradient descent code and associated plotting. --- Dipolar-Gas-Simulator/+Plotter/plotLive.m | 64 ++++++++++++++++- .../@DipolarGas/runGradientDescent.m | 71 ++++++++++--------- 2 files changed, 101 insertions(+), 34 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Plotter/plotLive.m b/Dipolar-Gas-Simulator/+Plotter/plotLive.m index cb6e194..8d98648 100644 --- a/Dipolar-Gas-Simulator/+Plotter/plotLive.m +++ b/Dipolar-Gas-Simulator/+Plotter/plotLive.m @@ -15,7 +15,7 @@ function plotLive(psi,Params,Transf,Observ,SimulationMode) n = abs(psi).^2; switch SimulationMode - case {'ImaginaryTimeEvolution', 'EnergyMinimization'} + case 'ImaginaryTimeEvolution' %Plotting fig = figure(1); fig.WindowState = 'maximized'; @@ -76,6 +76,68 @@ function plotLive(psi,Params,Transf,Observ,SimulationMode) ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); title('Chemical Potential', 'FontSize', 14); grid on + + case 'EnergyMinimization' + %Plotting + fig = figure(1); + fig.WindowState = 'maximized'; + clf + % set(gcf,'Position', [100, 100, 1600, 900]) + t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid + + nexttile; + nxz = squeeze(trapz(n*dy,2)); + nyz = squeeze(trapz(n*dx,1)); + nxy = squeeze(trapz(n*dz,3)); + plotxz = pcolor(x,z,nxz'); + set(plotxz, 'EdgeColor', 'none'); + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + % cbar1.Ticks = []; % Disable the ticks + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(x,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plotyz = pcolor(y,z,nyz'); + set(plotyz, 'EdgeColor', 'none'); + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + % cbar1.Ticks = []; % Disable the ticks + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plotxy = pcolor(x,y,nxy'); + set(plotxy, 'EdgeColor', 'none'); + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + % cbar1.Ticks = []; % Disable the ticks + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plot(Observ.theta,'-b') + ylabel('$\theta_{opt}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Optimal angle', 'FontSize', 14); + grid on + + nexttile; + plot(Observ.EVec,'-b') + ylabel('$E_{tot}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Total Energy', 'FontSize', 14); + grid on + + nexttile; + plot(Observ.mucVec,'-b') + ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Chemical Potential', 'FontSize', 14); + grid on case 'RealTimeEvolution' diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m index c0b7110..e65457c 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m @@ -94,6 +94,7 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) i = 1; Observ.residual = 1; Observ.res = 1; + Observ.theta = 0; % Initialize the PrematureExitFlag to false PrematureExitFlag = false; @@ -109,24 +110,23 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) % Compute gradient J = compute_gradient(psi, Params, Transf, VDk, V); % Calculate chemical potential - muchem = real(psi(:)' * J(:)) / sum(abs(psi(:)).^2); + muchem = real(inner_product(psi, J, Transf)) / inner_product(psi, psi, Transf); % 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 + % Energy convergence check every 100 steps + if mod(i,100) == 0 && length(Observ.EVec) > 1 energydifference = abs(Observ.EVec(end) - Observ.EVec(end-1)); + if energydifference <= epsilon + disp('Tolerance reached: Energy difference is below the specified epsilon.'); + PrematureExitFlag = true; + break; + end 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.MaxIterationsForGD % If set maximum number of iterations reached, then exit + + if i >= this.MaxIterationsForGD disp('Maximum number of iterations for CGD reached.'); - PrematureExitFlag = true; % Set flag to indicate premature exit + PrematureExitFlag = true; break; end @@ -135,33 +135,34 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) d = -residual; else % Compute beta via Polak-Ribiere and create new direction residual_new = residual; - beta = compute_beta(residual_new, residual_old); + beta = compute_beta(residual_new, residual_old, Transf); d = -residual_new + beta * p_old; end residual_old = residual; - p = d - (real(d(:)' * psi(:) * Transf.dx * Transf.dy * Transf.dz) .* psi); + proj = inner_product(d, psi, Transf); + p = d - (proj * 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 - gamma = 1 / sqrt(sum(abs(p(:)).^2) * Transf.dx * Transf.dy * Transf.dz); - psi = (cos(theta).*psi) + (gamma*sin(theta).*p); + gamma = 1 / sqrt(inner_product(p, p, Transf)); + psi = (cos(theta).*psi) + (sin(theta).*(p*gamma)); + % Normalize psi - Norm = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz; - psi = sqrt(Params.N) * psi / sqrt(Norm); + psi = sqrt(Params.N) * psi / sqrt(abs(inner_product(psi, psi, Transf))); i = i + 1; % Calculate chemical potential with new psi - muchem = real(psi(:)' * J(:)) / sum(abs(psi(:)).^2); + J = compute_gradient(psi, Params, Transf, VDk, V); + muchem = real(inner_product(psi, J, Transf)) / inner_product(psi, psi, Transf); if mod(i,100) == 0 % Collect Energy value - E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; + E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / inner_product(psi, psi, Transf); Observ.EVec = [Observ.EVec E]; % Collect Chemical potential value @@ -171,13 +172,15 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem); Observ.residual = [Observ.residual res]; Observ.res_idx = Observ.res_idx + 1; + + % Collect calculated thetas + Observ.theta = [Observ.theta theta]; % Live plotter if this.PlotLive Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode) drawnow end - save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V'); end end @@ -190,13 +193,11 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ) end % Change in Energy - E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; + E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / inner_product(psi, psi, Transf); Observ.EVec = [Observ.EVec E]; - 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!'); + disp('Completed and saved!'); end end @@ -261,15 +262,19 @@ 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); - gamma = 1 / sqrt(sum(abs(p(:)).^2) * Transf.dx * Transf.dy * Transf.dz); - numerator = gamma * real(p(:)' * Hpsi(:) * Transf.dx * Transf.dy * Transf.dz); - denominator = muchem - (gamma^2 * (real(p(:)' * Hp(:) * Transf.dx * Transf.dy * Transf.dz) + real(g(:)' * p(:) * Transf.dx * Transf.dy * Transf.dz))); + gamma = 1 / sqrt(inner_product(p, p, Transf)); + numerator = gamma * real(inner_product(p, Hpsi, Transf)); + denominator = muchem - (gamma^2 * (real(inner_product(p, Hp, Transf)) + real(inner_product(g, p, Transf)))); theta = numerator / denominator; end % Optimal step size via Polak-Ribiere -function beta = compute_beta(residual_new, residual_old) - beta = (residual_new(:)' * (residual_new(:) - residual_old(:))) / sum(abs(residual_old(:)).^2); - beta = max(0,beta); -end \ No newline at end of file +function beta = compute_beta(residual_new, residual_old, Transf) + beta = inner_product(residual_new, (residual_new - residual_old), Transf) / inner_product(residual_old, residual_old, Transf); + beta = max(0, real(beta)); +end + +function s = inner_product(u, v, Transf) + s = sum(conj(u(:)) .* v(:)) * Transf.dx * Transf.dy * Transf.dz; +end