Included BdG solver in to the Simulator, added plotting routines for the 2D solver, fixed bugs so that the 2D solver runs to completion.

Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction2D.m

@@ -0,0 +1,81 @@
+function visualizeGSWavefunction2D(run_index)
+    
+    set(0,'defaulttextInterpreter','latex')
+    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
+    
+    folder_path = './Data';
+
+    % Load data
+    Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ');
+    
+    Params = Data.Params;
+    Transf = Data.Transf; 
+    Observ = Data.Observ; 
+    if isgpuarray(Data.psi)
+        psi = gather(Data.psi);
+    else
+        psi = Data.psi;
+    end
+    
+    if isgpuarray(Data.Observ.residual)
+        Observ.residual = gather(Data.Observ.residual);
+    else
+        Observ.residual = Data.Observ.residual;
+    end
+
+    % Set long format for output
+    format long
+    
+    % Coordinates in micrometers
+    x = Transf.x * Params.l0 * 1e6; 
+    y = Transf.y * Params.l0 * 1e6; 
+    
+    % Plotting
+    figure('Position', [100, 100, 1600, 900]);
+    clf
+    
+    % Compute probability density |psi|^2
+    n = abs(psi).^2;
+    
+    % Plot |psi(x,y)|^2 (density in x and y plane)
+    subplot('Position', [0.05, 0.55, 0.28, 0.4])
+    pcolor(x, y, n');
+    shading interp;
+    colorbar;
+    xlabel('$x$ [$\mu$m]', 'FontSize', 14); ylabel('$y$ [$\mu$m]', 'FontSize', 14);
+    title('$|\Psi_{xy}|^2$', 'FontSize', 14);
+    
+    % Plot real part of psi in the x-y plane
+    subplot('Position', [0.36, 0.55, 0.28, 0.4])
+    pcolor(x, y, real(psi)');
+    shading interp;
+    colorbar;
+    xlabel('$x$ [$\mu$m]', 'FontSize', 14); ylabel('$y$ [$\mu$m]', 'FontSize', 14);
+    title('Re$\{\Psi_{xy}\}$', 'FontSize', 14);
+
+    % Plot imaginary part of psi in the x-y plane
+    subplot('Position', [0.67, 0.55, 0.28, 0.4])
+    pcolor(x, y, imag(psi)');
+    shading interp;
+    colorbar;
+    xlabel('$x$ [$\mu$m]', 'FontSize', 14); ylabel('$y$ [$\mu$m]', 'FontSize', 14);
+    title('Im$\{\Psi_{xy}\}$', 'FontSize', 14);
+
+    % Plot residual (time steps vs -log10(residual))
+    subplot('Position', [0.05, 0.05, 0.26, 0.4]);
+    plot(-log10(Observ.residual), '-b')
+    ylabel('$-\mathrm{log}_{10}(r)$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
+    title('Residual', 'FontSize', 14);
+
+    % Plot total energy over time
+    subplot('Position', [0.36, 0.05, 0.26, 0.4]);
+    plot(Observ.EVec, '-b')
+    ylabel('$E$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
+    title('Total Energy', 'FontSize', 14);
+
+    % Plot chemical potential over time
+    subplot('Position', [0.67, 0.05, 0.26, 0.4]);
+    plot(Observ.mucVec, '-b')
+    ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
+    title('Chemical Potential', 'FontSize', 14);
+end

Dipolar-Gas-Simulator/+Plotter/visualizeSpace.m

@@ -23,19 +23,19 @@ function visualizeSpace(Transf)
     end
     
     % set the axes labels' properties
-    xlabel(gca, {'$x / l_o ~ (m)$'}, ...
+    xlabel(gca, {'$x / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
             'FontWeight', 'normal', ...
             'FontAngle', 'normal')
-    ylabel(gca, {'$y / l_o ~ (m)$'}, ...
+    ylabel(gca, {'$y / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
             'FontWeight', 'normal', ...
             'FontAngle', 'normal')
-    zlabel(gca, {'$z / l_o ~ (m)$'}, ...
+    zlabel(gca, {'$z / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
@@ -63,19 +63,19 @@ function visualizeSpace(Transf)
     end
     
     % set the axes labels' properties
-    xlabel(gca, {'$k_x / l_o ~ (m^{-1})$'}, ...
+    xlabel(gca, {'$k_x / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
             'FontWeight', 'normal', ...
             'FontAngle', 'normal')
-    ylabel(gca, {'$k_y / l_o ~ (m^{-1})$'}, ...
+    ylabel(gca, {'$k_y / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
             'FontWeight', 'normal', ...
             'FontAngle', 'normal')
-    zlabel(gca, {'$k_z / l_o ~ (m^{-1})$'}, ...
+    zlabel(gca, {'$k_z / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...

Dipolar-Gas-Simulator/+Plotter/visualizeSpace2D.m

@@ -26,13 +26,13 @@ function visualizeSpace2D(Transf)
     axis equal  % Ensures equal scaling for both axes
     
     % Set axes labels
-    xlabel(gca, {'$x / l_o ~ (m)$'}, ...
+    xlabel(gca, {'$x / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
             'FontWeight', 'normal', ...
             'FontAngle', 'normal')
-    ylabel(gca, {'$y / l_o ~ (m)$'}, ...
+    ylabel(gca, {'$y / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
@@ -62,13 +62,13 @@ function visualizeSpace2D(Transf)
     axis equal  % Ensure equal scaling for both axes
     
     % Set axes labels
-    xlabel(gca, {'$k_x / l_o ~ (m^{-1})$'}, ...
+    xlabel(gca, {'$k_x / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...
             'FontWeight', 'normal', ...
             'FontAngle', 'normal')
-    ylabel(gca, {'$k_y / l_o ~ (m^{-1})$'}, ...
+    ylabel(gca, {'$k_y / l_o$'}, ...
             'Interpreter', 'latex', ...
             'FontName', 'Times New Roman', ...
             'FontSize', 14, ...

Dipolar-Gas-Simulator/+Plotter/visualizeWavefunction2D.m

@@ -0,0 +1,41 @@
+function visualizeWavefunction2D(psi, Params, Transf)
+
+    set(0,'defaulttextInterpreter','latex')
+    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
+    
+    format long
+    
+    height = 10; 
+    width = 30;
+    figure(1)
+    clf
+    set(gcf, 'Units', 'centimeters')
+    set(gcf, 'Position', [2 8 width height])
+    set(gcf, 'PaperPositionMode', 'auto')
+    
+    % Axes scaling and coordinates in micrometers
+    x = Transf.x * Params.l0 * 1e6; 
+    y = Transf.y * Params.l0 * 1e6; 
+
+    % Compute probability density |psi|^2
+    n = abs(psi).^2;
+
+    % Plot |psi(x,y)|^2
+    subplot(1, 2, 1)
+    pcolor(x, y, n')
+    shading interp;  % Smooth shading
+    colorbar;        % Show color bar
+    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)
+
+    % Plot real part of psi
+    subplot(1, 2, 2)
+    pcolor(x, y, real(psi)')
+    shading interp;
+    colorbar;
+    xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
+    ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
+    title('Re$\{\Psi(x,y)\}$', 'Interpreter', 'latex', 'FontSize', 14)
+
+end

Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m

@@ -14,8 +14,8 @@ function [psi,V,VDk] = initialize(this,Params,VParams,Transf)
     else
         VDk = this.Calculator.calculateVDkWithCutoff(rcut, Transf, Params, VParams);
         save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
-    end
         fprintf('Computed and saved DDI potential in Fourier space with cutoff.\n')
+    end
     
     % == Setting up the initial wavefunction == %
     psi = this.setupWavefunction(Params,Transf);

Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m

@@ -27,7 +27,7 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
     EVar            = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu);
 
     pb = Helper.ProgressBar();
-    pb.run('Running simulation in imaginary time: ');
+    pb.run('Propagating in imaginary time: ');
 
     while t_idx < Params.sim_time_cut_off
 
@@ -69,15 +69,12 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
             Observ.residual = [Observ.residual res];
             Observ.res_idx  = Observ.res_idx + 1;
     
-            save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
-            
             %Adaptive time step -- Careful, this can quickly get out of control
             relres          = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
 
             if relres <1e-4
                 if AdaptIdx > 3 && abs(dt) > Params.mindt
                     dt = dt / 2;
-                    fprintf('Time step changed to '); disp(dt);
                     AdaptIdx = 0;
                 elseif AdaptIdx > 3 && abs(dt) < Params.mindt
                     break
@@ -95,6 +92,7 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
             break
         end
         t_idx=t_idx+1;
+        pb.run(100*t_idx/Params.sim_time_cut_off);
     end
     
     % Change in Energy
@@ -110,8 +108,5 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
     % Chemical potential
     Observ.mucVec   = [Observ.mucVec muchem];
 
-    pb.run(' - Job Completed!');
+    pb.run(' - Done!');
-    disp('Saving data...');
-    save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
-    disp('Save complete!');
 end

Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m

@@ -23,13 +23,13 @@ function [Params, Transf, psi, V, VDk] = run(this)
     [psi,V,VDk]               = this.initialize(Params,VParams,Transf);
     ells(1)                   = VarArray(1); 
     nus(1)                    = VarArray(2);
-    E_Var                     = @(x) this.Calculator.calculateVariationalEnergy(psi,Params,x,Transf,VDk,V)/Params.N;
+    E_Var                     = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
     E_vs_iter(1)              = E_Var([ells(1) nus(1)]);
 
     t_idx                     = 1; % Start at t = 0;
        
     for nn = 1:Params.SelfConIter
-        Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
+        Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = [];
     Observ.res_idx = 1;
 
         VParams.ell           = VarArray(1);
@@ -47,11 +47,11 @@ function [Params, Transf, psi, V, VDk] = run(this)
         psi                   = psi + 4*noise;
 
         % --- Run --- 
-        [psi,Observ]          = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ);
+        [psi]                 = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ);
         psi                   = gather(psi);
         
         % --- Constrained minimization --- 
-        E_Var                 = @(x) this.Calculator.calculateVariationalEnergy(psi,Params,x,Transf,V)/Params.N;
+        E_Var                 = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
         VarArray              = fmincon(E_Var,VarArray,[],[],[],[],[Params.ell_lower, Params.nu_lower],[Params.ell_upper, Params.nu_upper],[],fminconoptions);
         
         % --- Convergence check --- 
@@ -61,13 +61,16 @@ function [Params, Transf, psi, V, VDk] = run(this)
         relnudiff             = abs(nus(nn+1)-nus(nn))/nus(nn);
         E_vs_iter             = [E_vs_iter E_Var(VarArray)];
 
-        save(sprintf('./Data/Run_%03i/psi_gs_%i.mat',runIdx,nn),'psi','Observ','Transf','Params','VDk','V','VarArray');
+        save(sprintf('./Data/Run_%03i/psi_gs_%i.mat',Params.njob),'psi','Observ','Transf','Params','VDk','V','VarArray');
 
         if relelldiff < Params.ellcutoff && relnudiff < Params.nucutoff
             break
         end
     end
-    save(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','Observ','Transf','Params','VDk','V','VarArray');
+    
-    delete(sprintf('./Data/Run_%03i/psi_gs_*.mat',runIdx))
+    disp('Saving data...');
+    save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','Observ','Transf','Params','VDk','V','VarArray');
+    disp('Save complete!');
+    delete(sprintf('./Data/Run_%03i/psi_gs_*.mat',Params.njob))
 
 end