diff --git a/Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m b/Dipolar-Gas-Simulator/+BdGSolver/solveBogoliubovdeGennesIn2D.m similarity index 100% rename from Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m rename to Dipolar-Gas-Simulator/+BdGSolver/solveBogoliubovdeGennesIn2D.m diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction2D.m b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction2D.m new file mode 100644 index 0000000..27cd369 --- /dev/null +++ b/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 diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeSpace.m b/Dipolar-Gas-Simulator/+Plotter/visualizeSpace.m index dfe5c4b..bfe0c5f 100644 --- a/Dipolar-Gas-Simulator/+Plotter/visualizeSpace.m +++ b/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, ... diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeSpace2D.m b/Dipolar-Gas-Simulator/+Plotter/visualizeSpace2D.m index 57734dc..58d8513 100644 --- a/Dipolar-Gas-Simulator/+Plotter/visualizeSpace2D.m +++ b/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, ... diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeWavefunction2D.m b/Dipolar-Gas-Simulator/+Plotter/visualizeWavefunction2D.m new file mode 100644 index 0000000..7a78515 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m index 519f053..d958c62 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m +++ b/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'); + fprintf('Computed and saved DDI potential in Fourier space with cutoff.\n') end - fprintf('Computed and saved DDI potential in Fourier space with cutoff.\n') % == Setting up the initial wavefunction == % psi = this.setupWavefunction(Params,Transf); diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m index 47f1fd5..8e668f4 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m +++ b/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!'); - 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!'); + pb.run(' - Done!'); end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m index 6de64b8..3809cc3 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m @@ -23,14 +23,14 @@ 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.res_idx = 1; + Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = []; + Observ.res_idx = 1; VParams.ell = VarArray(1); VParams.nu = VarArray(2); @@ -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