Major bugfixes - fmincon still fails.

This commit is contained in:
Karthik 2024-11-22 00:04:27 +01:00
parent 9577bafff4
commit a9beb317c9
10 changed files with 49 additions and 43 deletions

View File

@ -74,21 +74,21 @@ OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [100, 100];
OptionsStruct.TimeStepSize = 100E-6; % in s
OptionsStruct.TimeStepSize = 500E-6; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.VariationalWidth = 5.5;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1;
OptionsStruct.JobNumber = 2;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data_StripePhase';

View File

@ -9,7 +9,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;
@ -32,9 +32,9 @@ OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.VariationalWidth = 4;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.PlotLive = false;
@ -51,7 +51,7 @@ solver.Potential = pot.trap();
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
%}
%% - Create Variational2D and Calculator object with specified options
OptionsStruct = struct;
@ -64,7 +64,7 @@ OptionsStruct.ScatteringLength = 100.0289;
OptionsStruct.TrapFrequencies = [10, 10, 72.4];
OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [256, 256];
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [100, 100];
OptionsStruct.TimeStepSize = 100E-6; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
@ -74,10 +74,10 @@ OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.VariationalWidth = 5;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 2;
@ -93,7 +93,7 @@ solver.Potential = pot.trap();
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
%{
%% - Create Variational2D and Calculator object with specified options
OptionsStruct = struct;
@ -116,9 +116,9 @@ OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.VariationalWidth = 6.5;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthUpperBound = 30;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.PlotLive = false;
@ -134,4 +134,5 @@ pot = VariationalSolver2D.Potentials(options{:
solver.Potential = pot.trap();
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
[Params, Transf, psi, V, VDk] = solver.run();
%}

View File

@ -13,20 +13,20 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V)
Phi = real(ifftn(frho.*VDk));
Eddi = (Params.gdd*Phi.*abs(psi).^2)/(sqrt(2*pi)*VParams.ell);
%Kinetic energy
% Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2;
%Potential energy
% Potential energy
Epot = V.*abs(psi).^2;
%Contact interactions
% Contact interactions
Eint = g_eff*abs(psi).^4;
%Quantum fluctuations
% Quantum fluctuations
Eqf = gamma_eff*abs(psi).^5;
%Total energy
% Total energy
muchem = Ekin + Ez*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; %
muchem = muchem / Params.N; %Only use if psi is normalized to N
end

View File

@ -11,7 +11,7 @@ function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V)
% DDIs
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);%
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);
Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy;
% Kinetic energy

View File

@ -22,7 +22,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
VariationalWidth;
WidthLowerBound;
WidthUpperBound;
WidthCutoff;
VariationalWidthTolerance;
VariationalEnergyTolerance;
Calculator;
@ -77,7 +78,9 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'WidthUpperBound', 12,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'WidthCutoff', 1e-2,...
addParameter(p, 'VariationalWidthTolerance', 1e-2,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'VariationalEnergyTolerance', 1e-2,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'JobNumber', 1,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
@ -113,7 +116,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
this.VariationalWidth = p.Results.VariationalWidth;
this.WidthLowerBound = p.Results.WidthUpperBound;
this.WidthUpperBound = p.Results.WidthUpperBound;
this.WidthCutoff = p.Results.WidthCutoff;
this.VariationalWidthTolerance = p.Results.VariationalWidthTolerance;
this.VariationalEnergyTolerance = p.Results.VariationalEnergyTolerance;
this.PlotLive = p.Results.PlotLive;
this.JobNumber = p.Results.JobNumber;

View File

@ -31,10 +31,6 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
Ez = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
pb = Helper.ProgressBar();
fprintf('\n')
pb.run('Propagating wavefunction in imaginary time: ');
while t_idx < Params.sim_time_cut_off
% kin
@ -60,7 +56,7 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans
muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
%Plotting loop
% Intermittent saving and plotting
if mod(t_idx,500) == 0
% Change in Energy
@ -108,7 +104,6 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans
break
end
t_idx=t_idx+1;
pb.run(100*t_idx/Params.sim_time_cut_off);
end
% Change in Energy
@ -124,6 +119,6 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans
% Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
pb.run(' - Run Completed!');
clear pb.run
disp('Run Completed!');
end

View File

@ -18,6 +18,7 @@ function [Params, Transf, psi, V, VDk] = run(this)
% Relative cutoffs
VParams.ellcutoff = Params.ellcutoff;
VParams.evarcutoff = Params.evarcutoff;
% --- Initialize ---
mkdir(sprintf(this.SaveDirectory))
@ -34,7 +35,7 @@ function [Params, Transf, psi, V, VDk] = run(this)
t_idx = 1; % Start at t = 0;
[~,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate VDk with new value for the variational parameter, keep new psi at the end of loop to increase likelihood of faster convergence
[~,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter
% --- Adding some noise ---
% Noise added in every iteration to ensure it is not stuck in some local minimum
@ -57,8 +58,9 @@ function [Params, Transf, psi, V, VDk] = run(this)
ells = [ells VParams.ell];
relelldiff = abs(ells(nn+1)-ells(nn))/ells(nn);
E_vs_iter = [E_vs_iter E_Var(VParams.ell)];
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs_%i.mat'),Params.njob),'psi','Observ','Transf','Params','VDk','V','VParams');
relevardiff = abs(E_vs_iter(nn+1)-E_vs_iter(nn))/E_vs_iter(nn);
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs_%i.mat'),Params.njob,nn),'psi','Observ','Transf','Params','VDk','V','VParams');
%Plotting
if this.PlotLive
@ -70,9 +72,10 @@ function [Params, Transf, psi, V, VDk] = run(this)
drawnow
end
if relelldiff < Params.ellcutoff
if relelldiff < Params.ellcutoff && relevardiff < Params.evarcutoff
break
end
end
VParams.ells = ells;

View File

@ -61,7 +61,8 @@ function [Params] = setupParameters(this)
Params.ell_upper = this.WidthUpperBound;
% Relative cutoffs
Params.ellcutoff = this.WidthCutoff;
Params.ellcutoff = this.VariationalWidthTolerance;
Params.evarcutoff = this.VariationalEnergyTolerance;
% ================ Parameters defined by those above ================ %

View File

@ -1,5 +1,7 @@
function [psi] = setupWavefunction(~,Params,Transf)
format long
X = Transf.X;
Y = Transf.Y;

View File

@ -5,10 +5,10 @@
# Request number of nodes and GPU for job
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --gres=gpu:A100:1
#SBATCH --gres=gpu:A40:1
#SBATCH --mem=8G
# Estimated wallclock time for job
#SBATCH --time=16:00:00
#SBATCH --time=03:00:00
#SBATCH --job-name=simulation
#SBATCH --error=simulation.err
#SBATCH --output=simulation.out