karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Latest scripts - tweaked loop, switched to Polak-Ribiere, included live plotting for the energy minimizer routine.

Browse Source
This commit is contained in:
Karthik 2025-03-31 19:54:58 +02:00
parent abe402313e
commit 7aa4c7801b
4 changed files with 53 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Dipolar-Gas-Simulator/+Plotter/visualizeWavefunction.m
View File

@ -9,6 +9,10 @@ function visualizeWavefunction(psi,Params,Transf)
z = Transf.z*Params.l0*1e6;
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
if isgpuarray(psi)
psi = gather(psi);
end
%Plotting
height = 10;

15
Dipolar-Gas-Simulator/+Scripts/run_locally.m
View File

@ -17,13 +17,13 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [128, 128, 64];
OptionsStruct.NumberOfGridPoints = [64, 64, 32];
OptionsStruct.Dimensions = [18, 18, 18];
OptionsStruct.UseApproximationForLHY = true;
OptionsStruct.IncludeDDICutOff = true;
OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
OptionsStruct.MaxIterationsForCGD = 10;
OptionsStruct.MaxIterationsForCGD = 1E3;
OptionsStruct.TimeStepSize = 1E-4; % in s
OptionsStruct.MinimumTimeStepSize = 2E-10; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
@ -52,10 +52,6 @@ sim.Potential = pot.trap(); % + pot.repulsive_chopstick(
Plotter.visualizeTrapPotential(sim.Potential,Params,Transf)
%% - Plot initial wavefunction
Plotter.visualizeWavefunction(psi,Params,Transf)
%%
SaveDirectory = './Results/Data_3D/AspectRatio2';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% - Plot GS wavefunction
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7';
@ -71,12 +67,13 @@ 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_8';
SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_85';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD';
SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth';
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)

50
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m
View File

@ -1,18 +1,25 @@
function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V)
function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V,Observ)
format long;
% Define the function handle
f = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N;
f = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N;
% Convergence Criteria:
Epsilon = 1E-5;
Epsilon = 1E-5;
% Iteration Counter:
i = 1;
i = 1;
Observ.residual = 1;
Observ.res = 1;
% Initialize the PrematureExitFlag to false
PrematureExitFlag = false;
if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ)
drawnow
end
% Minimization Loop
while true
@ -51,6 +58,30 @@ function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V)
% Store old gradient
J_old = J;
i = i + 1;
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
if mod(i,10) == 0
% Change in Energy
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
% Normalized residuals
res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
Observ.residual = [Observ.residual res];
Observ.res_idx = Observ.res_idx + 1;
if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ)
drawnow
end
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
end
end
% Check if loop ended prematurely
@ -61,10 +92,13 @@ function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V)
fprintf('Number of Iterations for Convergence: %d\n\n', i);
end
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
% Change in Energy
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
disp('Saving data...');
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Transf','Params','VDk','V');
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
disp('Save complete!');
end
@ -105,7 +139,7 @@ function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V)
rho = 0.5; % Step size reduction factor
c = 1E-4; % Armijo condition constant
max_iter = 100; % Max iterations for backtracking
tol = 1E-8; % Tolerance for stopping
tol = 1E-4; % Tolerance for stopping
grad = compute_gradient(X, Params, Transf, VDk, V); % Compute gradient once
f_X = f(X); % Evaluate f(X) once

2
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m
View File

@ -20,7 +20,7 @@ function [Params, Transf, psi,V,VDk] = run(this)
mkdir(sprintf(this.SaveDirectory))
mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob))
if strcmp(this.SimulationMode, 'EnergyMinimization')
[psi] = this.minimizeEnergyFunctional(psi,Params,Transf,VDk,V);
[psi] = this.minimizeEnergyFunctional(psi,Params,Transf,VDk,V,Observ);
else
[psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ);
end