karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Latest versions of multiple scripts - modified plotting routines, execution scripts used to generate the latest phase diagrams.

Browse Source
This commit is contained in:
Karthik 2025-05-26 16:54:03 +02:00
parent 211724a779
commit dc5ccb67ba
10 changed files with 142 additions and 156 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
View File

@ -35,7 +35,6 @@ function visualizeGSWavefunction(folder_path, run_index)
%Plotting %Plotting
fig = figure(1); fig = figure(1);
fig.WindowState = 'maximized';
clf clf
set(gcf,'Position', [100, 100, 1600, 900]) set(gcf,'Position', [100, 100, 1600, 900])
t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid

88
Dipolar-Gas-Simulator/+Scripts/plotPhaseDiagramWithBoundaries.m
View File

@ -1,54 +1,58 @@
function plotPhaseDiagramWithBoundaries(PhaseDiagramMatrix, NumberOfAtoms, ScatteringLengths, PhaseBoundary_Untilted, TitleString) function plotPhaseDiagramWithBoundaries(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, PhaseBoundary, TitleString)
% plotPhaseDiagramWithBoundaries - Plots a phase diagram with overlaid interpolated boundaries
% % plotPhaseDiagramWithBoundaries - Plots a phase diagram with overlaid interpolated boundaries
% Inputs: %
% PhaseDiagramMatrix - 2D matrix of phase values % Inputs:
% NumberOfAtoms - Vector of atom numbers (x-axis) % PhaseDiagramMatrix - 2D matrix of phase values
% ScatteringLengths - Vector of scattering lengths (y-axis) % NumberOfAtoms - Vector of atom numbers (x-axis)
% PhaseBoundary_Untilted - Object containing selected boundary points % ScatteringLengths - Vector of scattering lengths (y-axis)
% TitleString - Title string for the plot % PhaseBoundary_Untilted - Object containing selected boundary points
% % TitleString - Title string for the plot
% This function displays the interpolated phase diagram and overlays %
% boundary curves extracted from user-selected points. % This function displays the interpolated phase diagram and overlays
% boundary curves extracted from user-selected points.
% Extract boundary point sets % Extract boundary point sets
rawPoints = PhaseBoundary_Untilted.SelectedPoints.values; rawPoints = PhaseBoundary.SelectedPoints.values;
if iscell(rawPoints) if iscell(rawPoints)
pointSets = rawPoints; pointSets = rawPoints;
else else
pointSets = {rawPoints}; pointSets = {rawPoints};
end end
% Interpolate phase diagram regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"];
[X, Y] = meshgrid(NumberOfAtoms, ScatteringLengths); regionColors = [
[Xq, Yq] = meshgrid(linspace(min(NumberOfAtoms), max(NumberOfAtoms), 500), ... 0.8 0.8 0.8;
linspace(min(ScatteringLengths), max(ScatteringLengths), 500)); 0.2 0.6 1.0;
Mq = interp2(X, Y, PhaseDiagramMatrix, Xq, Yq, 'nearest'); 0.2 0.8 0.2;
1.0 0.6 0.2;
0.8 0.2 0.8
];
% Plot [X, Y] = meshgrid(NUM_ATOMS_LIST, SCATTERING_LENGTH_RANGE);
figure('Color','w'); [Xq, Yq] = meshgrid(linspace(min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST), 500), ...
clf linspace(min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE), 500));
set(gcf, 'Position', [100, 100, 950, 800]); Mq = interp2(X, Y, M, Xq, Yq, 'nearest');
set(gcf, 'Renderer', 'opengl');
imagesc([min(NumberOfAtoms), max(NumberOfAtoms)], ... fig = figure('Color', 'w'); clf;
[min(ScatteringLengths), max(ScatteringLengths)], Mq); set(fig, 'Position', [100, 100, 1050, 800], 'Renderer', 'opengl');
set(gca, 'YDir', 'normal'); ax = axes('Parent', fig, 'Position', [0.1 0.15 0.75 0.8]);
colormap([
0.8 0.8 0.8; % Unmodulated imagesc(ax, [min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST)], ...
0.2 0.6 1.0; % SSD [min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE)], Mq);
0.2 0.8 0.2; % Stripe set(ax, 'YDir', 'normal');
1.0 0.6 0.2; % Labyrinth colormap(ax, regionColors);
0.8 0.2 0.8 % Honeycomb cb = colorbar('Ticks', 0:4, 'TickLabels', regionNames, 'FontSize', 12);
]); cb.Color = 'k'; % set colorbar text color
cb = colorbar('Ticks', 0:4, ... clim(ax, [0 4]);
'TickLabels', {'Unmodulated','SSD','Stripe','Labyrinth','Honeycomb'}, ... hold(ax, 'on');
'FontSize', 12); xlabel(ax, 'Number of Atoms', 'FontSize', 16);
cb.Color = 'k'; ylabel(ax, 'Scattering Length (\times a_o)', 'FontSize', 16,'Interpreter', 'tex');
clim([0 4]); t = title(ax, TitleString, 'FontSize', 18, 'Interpreter', 'tex');
xlabel('Number of Atoms'); t.Color = 'k'; % force title color to black
ylabel('Scattering Length a_s (a_0)'); set(ax, 'FontSize', 16, 'Color', 'none');
t = title(TitleString); axis(ax, 'tight');
t.Color = 'k'; grid(ax, 'on');
% Overlay curves % Overlay curves
hold on; hold on;

6
Dipolar-Gas-Simulator/+Scripts/plotSmoothedContourOnPhaseDiagram.m
View File

@ -1,5 +1,4 @@
function plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ... function plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString, varargin)
titleString, varargin)
interpMethod = 'makima'; interpMethod = 'makima';
showPoints = true; showPoints = true;
matFilePaths = {}; matFilePaths = {};
@ -52,7 +51,8 @@ function plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS
xlabel(ax, 'Number of Atoms', 'FontSize', 16); xlabel(ax, 'Number of Atoms', 'FontSize', 16);
ylabel(ax, 'Scattering Length (\times a_o)', 'FontSize', 16); ylabel(ax, 'Scattering Length (\times a_o)', 'FontSize', 16);
title(ax, titleString, 'FontSize', 18, 'Interpreter', 'tex'); t = title(ax, TitleString, 'FontSize', 18, 'Interpreter', 'tex');
t.Color = 'k'; % force title color to black
set(ax, 'FontSize', 16, 'Color', 'none'); set(ax, 'FontSize', 16, 'Color', 'none');
axis(ax, 'tight'); axis(ax, 'tight');
grid(ax, 'on'); grid(ax, 'on');

4
Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m
View File

@ -1,4 +1,4 @@
function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, titleString) function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString)
% Region labels and corresponding colors % Region labels and corresponding colors
regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"]; regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"];
regionColors = [ regionColors = [
@ -34,7 +34,7 @@ function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ti
% Axis formatting % Axis formatting
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16); xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16); ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16);
t = title(titleString, 'Interpreter', 'tex', 'FontSize', 18); t = title(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
t.Color = 'k'; % force title color to black t.Color = 'k'; % force title color to black
set(gca, 'FontSize', 16, ... set(gca, 'FontSize', 16, ...

11
Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker.m
View File

@ -23,7 +23,7 @@ function run_hybrid_worker(batchParams, batchIdx)
% Create unique save directory % Create unique save directory
jobName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', a_s, theta_deg, phi_deg, N_atoms); jobName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', a_s, theta_deg, phi_deg, N_atoms);
saveDir = fullfile('./Results/Data_3D/GradientDescent', jobName); saveDir = fullfile('./Results/Data_3D/PhaseDiagram', jobName);
if ~exist(saveDir, 'dir') if ~exist(saveDir, 'dir')
mkdir(saveDir); mkdir(saveDir);
end end
@ -43,9 +43,14 @@ function run_hybrid_worker(batchParams, batchIdx)
OptionsStruct.UseApproximationForLHY = true; OptionsStruct.UseApproximationForLHY = true;
OptionsStruct.IncludeDDICutOff = true; OptionsStruct.IncludeDDICutOff = true;
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'EnergyMinimization'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; % 'HeavyBall' | 'NonLinearCGD'
OptionsStruct.MaxIterationsForGD = 15000; OptionsStruct.MaxIterationsForGD = 15000;
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 2E5; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.010; OptionsStruct.NoiseScaleFactor = 0.010;
OptionsStruct.PlotLive = false; OptionsStruct.PlotLive = false;

109
Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary.m
View File

@ -1,74 +1,60 @@
function run_hybrid_worker_for_phase_boundary(batchParams, batchIdx) function run_hybrid_worker_for_phase_boundary(batchParams, batchIdx)
% Set up local cluster for parallel pool % Set up local cluster for parallel pool
cluster = parcluster('local'); cluster = parcluster('local');
nprocs = str2double(getenv('SLURM_CPUS_PER_TASK')); nprocs = str2double(getenv('SLURM_CPUS_PER_TASK'));
if isnan(nprocs), nprocs = feature('numcores'); end if isnan(nprocs), nprocs = feature('numcores'); end
tmpdir = fullfile(getenv('TMPDIR'), sprintf('matlab_job_%d', batchIdx)); tmpdir = fullfile(getenv('TMPDIR'), sprintf('matlab_job_%d', batchIdx));
if ~exist(tmpdir, 'dir'); mkdir(tmpdir); end if ~exist(tmpdir, 'dir'), mkdir(tmpdir); end
cluster.JobStorageLocation = tmpdir; cluster.JobStorageLocation = tmpdir;
pool = parpool(cluster, nprocs); pool = parpool(cluster, nprocs);
nJobs = size(batchParams, 1); % Extract single job parameters
parfor k = 1:nJobs initial_a_s = batchParams(1);
% Unpack parameter tuple theta_deg = batchParams(2);
initial_a_s = batchParams(k, 1); phi_deg = batchParams(3);
theta_deg = batchParams(k, 2); N_atoms = batchParams(4);
phi_deg = batchParams(k, 3);
N_atoms = batchParams(k, 4);
theta_rad = deg2rad(theta_deg); theta_rad = deg2rad(theta_deg);
phi_rad = deg2rad(phi_deg); phi_rad = deg2rad(phi_deg);
adjusted_a_s = initial_a_s; % Setup save directory
jobName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', initial_a_s, theta_deg, phi_deg, N_atoms);
saveDir = fullfile('./Results/Data_3D/PhaseDiagram', jobName);
initfile = fullfile(saveDir, 'psi_init.mat')
if ~exist(saveDir, 'dir'), mkdir(saveDir); end
if ~exist(initfile, 'file'), warning('Running without user-provided initial state'); end
% Create unique save directory % Parallelize attempts
jobName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', initial_a_s, theta_deg, phi_deg, N_atoms); MaxAttempts = 10;
saveDir = fullfile('./Results/Data_3D/GradientDescent', jobName); parfor attempt = 1:MaxAttempts
if ~exist(saveDir, 'dir') new_a_s = initial_a_s - 0.5 * (attempt - 1);
mkdir(saveDir);
end
srcFile = './Results/Data_3D/GradientDescent/psi_init.mat'; try
if exist(srcFile, 'file')
copyfile(srcFile, fullfile(saveDir, 'psi_init.mat'));
end
MaxAttempts = 10;
SuccessFlag = false;
AttemptCount = 0;
while ~SuccessFlag && AttemptCount < MaxAttempts
AttemptCount = AttemptCount + 1;
% Options for this run
OptionsStruct = struct; OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = N_atoms; OptionsStruct.NumberOfAtoms = N_atoms;
OptionsStruct.DipolarPolarAngle = theta_rad; OptionsStruct.DipolarPolarAngle = theta_rad;
OptionsStruct.DipolarAzimuthAngle = phi_rad; OptionsStruct.DipolarAzimuthAngle = phi_rad;
OptionsStruct.ScatteringLength = adjusted_a_s; OptionsStruct.ScatteringLength = new_a_s;
OptionsStruct.TrapFrequencies = [50, 20, 150]; OptionsStruct.TrapFrequencies = [50, 20, 150];
OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [128, 256, 128]; OptionsStruct.NumberOfGridPoints = [128, 256, 128];
OptionsStruct.Dimensions = [40, 80, 40]; OptionsStruct.Dimensions = [20, 40, 20];
OptionsStruct.UseApproximationForLHY = true; OptionsStruct.UseApproximationForLHY = true;
OptionsStruct.IncludeDDICutOff = true; OptionsStruct.IncludeDDICutOff = true;
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution';
OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; % 'HeavyBall' | 'NonLinearCGD' OptionsStruct.GradientDescentMethod = 'NonLinearCGD';
OptionsStruct.MaxIterationsForGD = 15000; OptionsStruct.MaxIterationsForGD = 15000;
OptionsStruct.TimeStepSize = 1E-3; % in s OptionsStruct.TimeStepSize = 1E-3;
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s OptionsStruct.MinimumTimeStepSize = 1E-6;
OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.TimeCutOff = 2E5;
OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 0.010; OptionsStruct.NoiseScaleFactor = 0.010;
OptionsStruct.PlotLive = false; OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 0; OptionsStruct.JobNumber = attempt - 1;
OptionsStruct.RunOnGPU = true; OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true; OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = saveDir; OptionsStruct.SaveDirectory = saveDir;
@ -79,35 +65,12 @@ function run_hybrid_worker_for_phase_boundary(batchParams, batchIdx)
pot = Simulator.Potentials(options{:}); pot = Simulator.Potentials(options{:});
sim.Potential = pot.trap(); sim.Potential = pot.trap();
NumberOfOutputs = 5; [~, ~, ~, ~, ~] = sim.run();
try
[Params, Transf, psi, ~, ~, stats] = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, saveDir);
if Scripts.determineDensityModulation(psi, Params, Transf) fprintf('SUCCESS | Attempt %d: a_s = %.2f\n', attempt, new_a_s);
SuccessFlag = true;
runDir = fullfile(saveDir, sprintf('Run_%03d', OptionsStruct.JobNumber));
psiFile = fullfile(runDir, 'psi_gs.mat');
if exist(psiFile, 'file') catch ME
Scripts.saveUpdatedScatteringLength(psiFile, adjusted_a_s, AttemptCount); fprintf('FAILURE | Attempt %d: %s\n', attempt, getReport(ME, 'basic'));
else
warning('Expected file %s not found. Cannot save final a_s.', psiFile);
end
else
adjusted_a_s = adjusted_a_s - 0.2; % Tweak as per your needs
end
catch ME
fprintf('ERROR in job %d (attempt %d):\n%s\n', k, AttemptCount, getReport(ME, 'extended'));
adjusted_a_s = adjusted_a_s - 0.2;
end
end
if SuccessFlag
fprintf('Batch %d | Job %d: a_s = %.2f, theta = %d°, phi = %d°, N = %d | Time = %.2f s\n', ...
batchIdx, k, adjusted_a_s, theta_deg, phi_deg, N_atoms, stats.runtime);
else
fprintf('Batch %d | Job %d FAILED after %d tries.\n', batchIdx, k, MaxAttempts);
end end
end end

11
Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary_wrapper.m
View File

@ -30,3 +30,14 @@ function run_hybrid_worker_for_phase_boundary_wrapper(batchIdx)
% Call worker with this batch of parameters % Call worker with this batch of parameters
Scripts.run_hybrid_worker_for_phase_boundary(batchParams, batchIdx); Scripts.run_hybrid_worker_for_phase_boundary(batchParams, batchIdx);
end end
function vals = parse_environmental_variable(varName, default)
str = getenv(varName);
if isempty(str)
vals = default;
elseif startsWith(str, '[')
vals = str2num(str); %#ok<ST2NM>
else
vals = eval(str); % Trust only controlled environments
end
end

32
Dipolar-Gas-Simulator/+Scripts/run_locally.m
View File

@ -557,7 +557,7 @@ SaveDirectory = './Results/Data_Full3D/PhaseDiagram/ImagTimePropagation/aS_8.312
JobNumber = 0; JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% %%
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta20/HighN/aS_9.562000e+01_theta_020_phi_000_N_1325000'; SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_9.250000e+01_theta_040_phi_000_N_916667';
JobNumber = 0; JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% Identify and count droplets %% Identify and count droplets
@ -661,7 +661,7 @@ SCATTERING_LENGTH_RANGE ="[75.00 76.09 77.18 78.27 79.36 80.00 81.04 82.08 83.12
NUM_ATOMS_LIST ="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]"; NUM_ATOMS_LIST ="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]";
%% Explore the states %% Explore the states
SCATTERING_LENGTH_RANGE = [80.00]; SCATTERING_LENGTH_RANGE = [97.71];
NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]; NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];
% NUM_ATOMS_LIST = [4183333 4387500 4591667 4795833 5000000]; % NUM_ATOMS_LIST = [4183333 4387500 4591667 4795833 5000000];
@ -672,7 +672,7 @@ for i = 1:length(SCATTERING_LENGTH_RANGE)
for j = 1:length(NUM_ATOMS_LIST) for j = 1:length(NUM_ATOMS_LIST)
N = NUM_ATOMS_LIST(j); N = NUM_ATOMS_LIST(j);
SaveDirectory = sprintf('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta20/HighN/aS_%.6e_theta_020_phi_000_N_%d', aS, N); SaveDirectory = sprintf('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_%.6e_theta_040_phi_000_N_%d', aS, N);
fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory); fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory);
% Call the plotting function % Call the plotting function
@ -756,42 +756,46 @@ Scripts.editPhaseDiagram(PhaseDiagramMatrix, ScatteringLengths, NumberOfAtoms)
%% Smoothen phase diagram %% Smoothen phase diagram
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat'); % Load M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat'); % Load M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ"; TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ";
Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString); Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString);
%% Select boundary points from original phase diagram %% Select boundary points from original phase diagram
load('./Results/Data_Full3D/PhaseDiagramUntilted.mat') load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
M = M; M = M;
N_atoms = round(NUM_ATOMS_LIST * 1E-5, 2); N_atoms = round(NUM_ATOMS_LIST * 1E-5, 2);
scatt_lengths = SCATTERING_LENGTH_RANGE; scatt_lengths = SCATTERING_LENGTH_RANGE;
savePath = './Results/Data_Full3D//BoundaryPoints/SelectedPoints_Untilted_1.mat'; savePath = './Results/Data_Full3D//BoundaryPoints/SelectedPoints_Theta20_1.mat';
Scripts.selectPhaseDiagramPoints(M, N_atoms, scatt_lengths, savePath); Scripts.selectPhaseDiagramPoints(M, N_atoms, scatt_lengths, savePath);
%% Interactively modify selected boundary points %% Interactively modify selected boundary points
load('./Results/Data_Full3D/PhaseDiagramUntilted.mat') load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ"; TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ... Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ...
TitleString, ... TitleString, ...
'./Results/Data_Full3D/BoundaryPoints/SelectedPoints_Untilted_1.mat', './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Untilted_2.mat', './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Untilted_3.mat', ... './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Theta20_1', './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Theta20_2', ...
'interpMethod', 'pchip', 'showPoints', true); 'interpMethod', 'pchip', 'showPoints', true);
%% Edit modified points %% Edit modified points
load('./Results/Data_Full3D/PhaseDiagramUntilted.mat') load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ"; TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ... Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ...
TitleString, ... TitleString, ...
'./Results/Data_Full3D//BoundaryPoints/ModifiedPoints_Untilted_1.mat', './Results/Data_Full3D//BoundaryPoints/ModifiedPoints_Untilted_2.mat', './Results/Data_Full3D//BoundaryPoints/ModifiedPoints_Untilted_3.mat', ... './Results/Data_Full3D/BoundaryPoints/ModifiedPoints_Theta20_1', './Results/Data_Full3D/BoundaryPoints/ModifiedPoints_Theta20_2', ...
'interpMethod', 'pchip', 'showPoints', true); 'interpMethod', 'pchip', 'showPoints', true);
%% Plot with interpolated phase boundary %% Plot with interpolated phase boundary
load('./Results/Data_Full3D/PhaseDiagramUntilted.mat') load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
PhaseDiagramMatrix = M; PhaseDiagramMatrix = M;
ScatteringLengths = SCATTERING_LENGTH_RANGE; ScatteringLengths = SCATTERING_LENGTH_RANGE;
NumberOfAtoms = NUM_ATOMS_LIST; NumberOfAtoms = NUM_ATOMS_LIST;
PhaseBoundary_Untilted = load("./Results/Data_Full3D/BoundaryPoints/PhaseBoundary_Untilted.mat"); PhaseBoundary = load("./Results/Data_Full3D/BoundaryPoints/PhaseBoundary_Tilted_Theta20.mat");
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ"; TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
Scripts.plotPhaseDiagramWithBoundaries(PhaseDiagramMatrix, NumberOfAtoms, ScatteringLengths, PhaseBoundary_Untilted, TitleString); Scripts.plotPhaseDiagramWithBoundaries(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, PhaseBoundary, TitleString);
%% Density modulation determination %% Density modulation determination

12
Dipolar-Gas-Simulator/submit_jobs.sh
View File

@ -1,11 +1,11 @@
# ----------- Define scan ranges ----------- # ----------- Define scan ranges -----------
# Use space-separated floating-point/integer values # Use space-separated floating-point/integer values
SCATTERING_LENGTH_RANGE="[79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0]" SCATTERING_LENGTH_RANGE="[75.00 76.09 77.18 78.27 79.36 80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29]"
POLAR_ANGLE_RANGE="[0.0 20.0 30.0]" POLAR_ANGLE_RANGE="[40.0]"
AZIMUTHAL_ANGLE_RANGE="[0.0]" AZIMUTHAL_ANGLE_RANGE="[0.0]"
NUM_ATOMS_LIST="[50000 55000 60000 65000 70000 75000 80000 85000 90000 95000 100000 105000 110000 115000 120000 125000 130000 135000 140000 145000 150000 155000 160000 165000]" NUM_ATOMS_LIST="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]"
CHUNK_SIZE=4 CHUNK_SIZE=1
# ----------- Count total combinations for SLURM array ----------- # ----------- Count total combinations for SLURM array -----------
@ -31,8 +31,8 @@ sbatch --export=SCATTERING_LENGTH_RANGE="$SCATTERING_LENGTH_RANGE",POLAR_ANGLE_R
#SBATCH --ntasks-per-node=1 #SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=8 #SBATCH --cpus-per-task=8
#SBATCH --gres=gpu:1 #SBATCH --gres=gpu:1
#SBATCH --mem=16G #SBATCH --mem=32G
#SBATCH --time=04:00:00 #SBATCH --time=06:00:00
#SBATCH --job-name=simulation #SBATCH --job-name=simulation
#SBATCH --error=simulation_%A_%a.err #SBATCH --error=simulation_%A_%a.err
#SBATCH --output=simulation_%A_%a.out #SBATCH --output=simulation_%A_%a.out

10
Dipolar-Gas-Simulator/submit_jobs_for_phase_boundary.sh
View File

@ -1,11 +1,11 @@
# ----------- Define scan ranges ----------- # ----------- Define scan ranges -----------
# Use space-separated floating-point/integer values # Use space-separated floating-point/integer values
SCATTERING_LENGTH_RANGE="[82.0 84.0 86.0 88.0 90.0 92.0 94.0 96.0 98.0 100.0]" SCATTERING_LENGTH_RANGE="[82.08 83.12 84.17 85.21 86.25 87.29 88.33]"
POLAR_ANGLE_RANGE="[0.0]" POLAR_ANGLE_RANGE="[0.0]"
AZIMUTHAL_ANGLE_RANGE="[0.0]" AZIMUTHAL_ANGLE_RANGE="[0.0]"
NUM_ATOMS_LIST="[50000 500000 950000 1400000 1850000 2300000 2750000 3200000 3650000 4100000 4550000 5000000]" NUM_ATOMS_LIST="[50000 54545 59091 63636 68182 72727 77273]"
CHUNK_SIZE=4 CHUNK_SIZE=1
# ----------- Count total combinations for SLURM array ----------- # ----------- Count total combinations for SLURM array -----------
@ -39,7 +39,7 @@ sbatch --export=SCATTERING_LENGTH_RANGE="$SCATTERING_LENGTH_RANGE",POLAR_ANGLE_R
#SBATCH --cpus-per-task=8 #SBATCH --cpus-per-task=8
#SBATCH --gres=gpu:1 #SBATCH --gres=gpu:1
#SBATCH --mem=16G #SBATCH --mem=16G
#SBATCH --time=04:00:00 #SBATCH --time=12:00:00
#SBATCH --job-name=simulation #SBATCH --job-name=simulation
#SBATCH --error=simulation_%A_%a.err #SBATCH --error=simulation_%A_%a.err
#SBATCH --output=simulation_%A_%a.out #SBATCH --output=simulation_%A_%a.out
@ -51,7 +51,7 @@ module load math/matlab/R2023a
echo "Initiating Job..." echo "Initiating Job..."
# Run MATLAB wrapper with this batch index # Run MATLAB wrapper with this batch index
matlab -nodisplay -nosplash -r "Scripts.run_hybrid_worker_wrapper(\$SLURM_ARRAY_TASK_ID)" matlab -nodisplay -nosplash -r "Scripts.run_hybrid_worker_for_phase_boundary_wrapper(\$SLURM_ARRAY_TASK_ID)"
echo "Job terminated successfully" echo "Job terminated successfully"