90 lines
4.2 KiB
Matlab
90 lines
4.2 KiB
Matlab
function [Params, Transf, psi, V, VDk] = run(this)
|
|
format long;
|
|
|
|
% --- Obtain simulation parameters ---
|
|
[Params] = this.setupParameters();
|
|
this.SimulationParameters = Params;
|
|
|
|
% --- Set up spatial grids and transforms ---
|
|
[Transf] = this.setupSpace(Params);
|
|
|
|
% --- Set up for Variational method ---
|
|
VParams.SelfConIter = Params.SelfConIter; % Max number of iterations to perform self-consistent calculation
|
|
VParams.ell = Params.ell; % initial [ell], ell is the "width" - psi ~ e^(z^2/ell^2)
|
|
|
|
% Window of optimization
|
|
VParams.ell_lower = Params.ell_lower;
|
|
VParams.ell_upper = Params.ell_upper;
|
|
|
|
% Relative cutoffs
|
|
VParams.ellcutoff = Params.ellcutoff;
|
|
|
|
% --- Initialize ---
|
|
mkdir(sprintf(this.SaveDirectory))
|
|
mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob))
|
|
fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8);
|
|
|
|
[psi,V,VDk] = this.initialize(Params,VParams,Transf);
|
|
ells(1) = VParams.ell;
|
|
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
|
|
E_vs_iter(1) = E_Var(ells(1));
|
|
|
|
for nn = 1:Params.SelfConIter
|
|
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = []; Observ.res_idx = 1;
|
|
|
|
t_idx = 1; % Start at t = 0;
|
|
|
|
[~,V,VDk] = this.initialize(Params,VParams,Transf); % Do not overwrite psi, susbequent iterations should use psi generated at the end of the loop to converge faster
|
|
% This is still needed however to recalculate 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
|
|
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation
|
|
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
|
r = normrnd(0,1,size(psi));
|
|
theta = rand(size(psi));
|
|
noise = r.*exp(2*pi*1i*theta);
|
|
psi = psi + Params.nsf*noise;
|
|
|
|
% --- Run ---
|
|
[psi, Observ] = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ, nn);
|
|
psi = gather(psi);
|
|
|
|
% --- Constrained minimization ---
|
|
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
|
|
VParams.ell = fmincon(E_Var,VParams.ell,[],[],[],[],Params.ell_lower,Params.ell_upper,[],fminconoptions);
|
|
|
|
% --- Convergence check ---
|
|
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');
|
|
|
|
%Plotting
|
|
if this.PlotLive
|
|
figure(2);
|
|
plot(ells,E_vs_iter,'LineStyle','none','Marker','o','MarkerSize',3,'Color','b','MarkerFaceColor','b')
|
|
xlabel('$\ell$', 'FontSize', 14)
|
|
ylabel('$E_{var}$', 'FontSize', 14)
|
|
title('Variational Energy', 'FontSize', 14);
|
|
drawnow
|
|
end
|
|
|
|
if relelldiff < Params.ellcutoff
|
|
break
|
|
end
|
|
end
|
|
|
|
VParams.ells = ells;
|
|
VParams.E_vs_iter = E_vs_iter;
|
|
|
|
fprintf('\n')
|
|
disp('Saving data...');
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','Observ','Transf','Params','VDk','V','VParams');
|
|
disp('Save complete!');
|
|
delete(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs_*.mat'),Params.njob))
|
|
|
|
end
|