Debugging attempt, minor modifications to plotting routines.

This commit is contained in:
Karthik 2024-11-20 09:49:39 +01:00
parent b8953f79a8
commit fde359df9e
6 changed files with 76 additions and 45 deletions

View File

@ -33,15 +33,21 @@ function visualizeGSWavefunction2D(folder_path, run_index)
clf
% Compute probability density |psi|^2
n = abs(psi).^2;
nxy = abs(psi).^2;
nxyScaled = nxy*(Params.add*10^6)^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);
ax1 = subplot('Position', [0.05, 0.55, 0.28, 0.4]);
plotxy = pcolor(x,y,nxyScaled');
set(plotxy, 'EdgeColor', 'none');
% shading interp; % Smooth shading
clim(ax1,[0.00,0.3]);
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
ylabel(cbar1,'$na_{dd}^2$','FontSize',16,'Rotation',270)
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 in the x-y plane
subplot('Position', [0.36, 0.55, 0.28, 0.4])

View File

@ -3,12 +3,13 @@ function analyzeRun2D(folder_path, run_index)
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
% Load data
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Observ','Transf','Params','VParams');
Params = Data.Params;
Transf = Data.Transf;
Observ = Data.Observ;
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Transf','Observ','Params','VParams');
Transf = Data.Transf;
Observ = Data.Observ;
Params = Data.Params;
VParams = Data.VParams;
if isgpuarray(Data.psi)
psi = gather(Data.psi);
else
@ -33,38 +34,53 @@ function analyzeRun2D(folder_path, run_index)
clf
% Compute probability density |psi|^2
n = abs(psi).^2;
nxy = abs(psi).^2;
nxyScaled = nxy*(Params.add*10^6)^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');
ax1 = subplot('Position', [0.05, 0.55, 0.28, 0.4]);
plotxy = pcolor(x,y,nxyScaled');
set(plotxy, 'EdgeColor', 'none');
% shading interp; % Smooth shading
clim(ax1,[0.00,0.3]);
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
ylabel(cbar1,'$na_{dd}^2$','FontSize',16,'Rotation',270)
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 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('$|\Psi_{xy}|^2$', 'FontSize', 14);
% Plot residual (time steps vs -log10(residual))
subplot('Position', [0.36, 0.55, 0.28, 0.4])
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])
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.67, 0.55, 0.28, 0.4])
% Plot residual (time steps vs -log10(residual))
subplot('Position', [0.05, 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.05, 0.05, 0.26, 0.4]);
subplot('Position', [0.67, 0.05, 0.26, 0.4]);
% Plot total energy over time
subplot('Position', [0.36, 0.05, 0.26, 0.4]);
plot(Observ.mucVec, '-b')
ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
title('Chemical Potential', 'FontSize', 14);
subplot('Position', [0.36, 0.05, 0.26, 0.4]);
% Plot chemical potential over time
subplot('Position', [0.67, 0.05, 0.26, 0.4]);
plot(VParams.ells,VParams.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);
end

View File

@ -84,8 +84,8 @@ OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1;
@ -111,10 +111,19 @@ Plotter.visualizeWavefunction2D(psi,Params,Transf)
%% - Plot GS wavefunction
Plotter.visualizeGSWavefunction2D(solver.SaveDirectory, solver.JobNumber)
%%
%% - Cluster Runs - Analysis
SaveDirectory = './Data_TriangularPhase';
JobNumber = 1;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
% analyzeRun2D(SaveDirectory, JobNumber)
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
Scripts.analyzeRun2D(SaveDirectory, JobNumber)
%%
%% - Cluster Runs - Analysis
SaveDirectory = './Data_StripePhase';
JobNumber = 2;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
Scripts.analyzeRun2D(SaveDirectory, JobNumber)
%% - Cluster Runs - Analysis
SaveDirectory = './Data_HoneycombPhase';
JobNumber = 3;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
% Scripts.analyzeRun2D(SaveDirectory, JobNumber)

View File

@ -34,8 +34,8 @@ OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 1;
@ -52,6 +52,7 @@ solver.Potential = pot.trap();
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
%{
%% - Create Variational2D and Calculator object with specified options
OptionsStruct = struct;
@ -76,8 +77,8 @@ OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 2;
@ -118,8 +119,8 @@ OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 20;
OptionsStruct.VariationalWidth = 5.7;
OptionsStruct.WidthLowerBound = 0.2;
OptionsStruct.WidthUpperBound = 20;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-3;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 3;
@ -134,4 +135,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

@ -34,9 +34,7 @@ function [Params, Transf, psi, V, VDk] = run(this)
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
[psi,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

View File

@ -8,7 +8,7 @@
#SBATCH --gres=gpu:A40:1
#SBATCH --mem=8G
# Estimated wallclock time for job
#SBATCH --time=12:00:00
#SBATCH --time=4:00:00
#SBATCH --job-name=simulation
#SBATCH --error=simulation.err
#SBATCH --output=simulation.out