From d5523fd0d90be7d81f1a999b562286387f1d598b Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 18 Nov 2024 18:06:14 +0100 Subject: [PATCH] Latest working version. --- Dipolar-Gas-Simulator/+Plotter/plotLive2D.m | 2 +- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 18 +++++++++------- .../+Scripts/run_on_cluster.m | 21 +++++++++++++------ .../@DipolarGas/DipolarGas.m | 8 +++++-- .../@DipolarGas/propagateWavefunction.m | 9 ++++---- .../+VariationalSolver2D/@DipolarGas/run.m | 2 +- .../@DipolarGas/setupParameters.m | 6 ++++-- .../@DipolarGas/setupWavefunction.m | 16 +++++++------- 8 files changed, 50 insertions(+), 32 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Plotter/plotLive2D.m b/Dipolar-Gas-Simulator/+Plotter/plotLive2D.m index d9cafa1..0791992 100644 --- a/Dipolar-Gas-Simulator/+Plotter/plotLive2D.m +++ b/Dipolar-Gas-Simulator/+Plotter/plotLive2D.m @@ -23,7 +23,7 @@ function plotLive2D(psi, Params, Transf, Observ, vrun) plotxy = pcolor(x,y,nxyScaled'); set(plotxy, 'EdgeColor', 'none'); % shading interp; % Smooth shading - clim(ax1); + clim(ax1,[0.00,0.3]); cbar1 = colorbar; cbar1.Label.Interpreter = 'latex'; ylabel(cbar1,'$na_{dd}^2$','FontSize',16,'Rotation',270) diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 2d18a99..66f584d 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -52,7 +52,7 @@ Plotter.visualizeGSWavefunction(Params.njob) OptionsStruct = struct; -OptionsStruct.NumberOfAtoms = 1E5; +OptionsStruct.NumberOfAtoms = 1.6E8; OptionsStruct.DipolarPolarAngle = 0; OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.ScatteringLength = 102.2518; % Critical point: 102.515; Triangular phase: 98.0676; Stripe phase: 102.2518; Honeycomb phase: 102.6441 @@ -60,18 +60,20 @@ OptionsStruct.ScatteringLength = 102.2518; % Critical point: 102.5 OptionsStruct.TrapFrequencies = [10, 10, 72.4]; OptionsStruct.TrapPotentialType = 'None'; -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [6.972, 6.972]; % Critical point: 6.996; Triangular phase: 7.5; Stripe phase: 6.972; Honeycomb phase: 6.239 for both for Atom Number fixed to 1E5 -OptionsStruct.TimeStepSize = 5E-6; % in s -OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [100, 100]; % Critical point: 6.996; Triangular phase: 7.5; Stripe phase: 6.972; Honeycomb phase: 6.239 for both for Atom Number fixed to 1E5 +OptionsStruct.TimeStepSize = 500E-6; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-04; +OptionsStruct.NoiseScaleFactor = 8; OptionsStruct.MaxIterations = 20; OptionsStruct.VariationalWidth = 4; OptionsStruct.WidthLowerBound = 0.2; OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-4; +OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.PlotLive = true; OptionsStruct.JobNumber = 1; @@ -90,6 +92,8 @@ solver.Potential = pot.trap(); %% - Plot numerical grid % Plotter.visualizeSpace2D(Transf) +%% - Plot trap potential +% Plotter.visualizeTrapPotential2D(solver.Potential,Params,Transf) %% - Plot initial wavefunction Plotter.visualizeWavefunction2D(psi,Params,Transf) %% - Plot GS wavefunction diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 820d96f..a74a7c3 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -2,21 +2,30 @@ OptionsStruct = struct; -OptionsStruct.NumberOfAtoms = 1E5; +OptionsStruct.NumberOfAtoms = 1.6E8; OptionsStruct.DipolarPolarAngle = 0; OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 98.0676; % Critical point: 102.515; Triangular phase: 98.0676; Stripe phase: 102.2518; Honeycomb phase: 102.6441 +OptionsStruct.ScatteringLength = 102.2518; % Critical point: 102.515; Triangular phase: 98.0676; Stripe phase: 102.2518; Honeycomb phase: 102.6441 OptionsStruct.TrapFrequencies = [10, 10, 72.4]; OptionsStruct.TrapPotentialType = 'None'; -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [7.5, 7.5]; % Critical point: 6.996; Triangular phase: 7.5; Stripe phase: 6.972; Honeycomb phase: 6.239 for both for Atom Number fixed to 1E5 -OptionsStruct.TimeStepSize = 500E-6; % in s -OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [100, 100]; % Critical point: 6.996; Triangular phase: 7.5; Stripe phase: 6.972; Honeycomb phase: 6.239 for both for Atom Number fixed to 1E5 +OptionsStruct.TimeStepSize = 500E-6; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-04; +OptionsStruct.NoiseScaleFactor = 8; +OptionsStruct.MaxIterations = 20; +OptionsStruct.VariationalWidth = 4; +OptionsStruct.WidthLowerBound = 0.2; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 1e-2; + +OptionsStruct.PlotLive = false; OptionsStruct.JobNumber = 1; OptionsStruct.RunOnGPU = true; OptionsStruct.SaveData = true; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m index 777e078..017a418 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m @@ -16,6 +16,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable TimeCutOff; EnergyTolerance; ResidualTolerance; + NoiseScaleFactor; MaxIterations; VariationalWidth; @@ -58,7 +59,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); addParameter(p, 'TimeStepSize', 5E-4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'MinimumTimeStepSize', 1e-6,... + addParameter(p, 'MinimumTimeStepSize', 1e-5,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'TimeCutOff', 2e6,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @@ -66,6 +67,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'ResidualTolerance', 1e-10,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'NoiseScaleFactor', 4,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'MaxIterations', 20,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'VariationalWidth', 4,... @@ -104,6 +107,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.TimeCutOff = p.Results.TimeCutOff; this.EnergyTolerance = p.Results.EnergyTolerance; this.ResidualTolerance = p.Results.ResidualTolerance; + this.NoiseScaleFactor = p.Results.NoiseScaleFactor; this.MaxIterations = p.Results.MaxIterations; this.VariationalWidth = p.Results.VariationalWidth; @@ -141,7 +145,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable ret = this.TimeCutOff; end function set.NumberOfAtoms(this, val) - assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!'); + assert(val <= 1E9, '!!Not time efficient to compute for atom numbers larger than 1,000,000,000!!'); this.NumberOfAtoms = val; end function ret = get.NumberOfAtoms(this) diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m index 2924d77..9bc3114 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m @@ -84,10 +84,10 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans % %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); - relE = abs(Observ.EVec(Observ.res_idx)-Observ.EVec(Observ.res_idx-1))/Observ.EVec(Observ.res_idx); - relmu = abs(Observ.mucVec(Observ.res_idx)-Observ.mucVec(Observ.res_idx-1))/Observ.mucVec(Observ.res_idx); - + relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx); + relE = abs(Observ.EVec(Observ.res_idx)-Observ.EVec(Observ.res_idx-1))/Observ.EVec(Observ.res_idx); + relmu = abs(Observ.mucVec(Observ.res_idx)-Observ.mucVec(Observ.res_idx-1))/Observ.mucVec(Observ.res_idx); + if relres <1e-4 if AdaptIdx > 3 && abs(dt) > Params.mindt dt = dt / 2; @@ -103,7 +103,6 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans AdaptIdx = 0; end end - if any(isnan(psi(:))) disp('NaNs encountered!') break diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m index be697ab..b044eb0 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m @@ -45,7 +45,7 @@ function [Params, Transf, psi, V, VDk] = run(this) r = normrnd(0,1,size(psi)); theta = rand(size(psi)); noise = r.*exp(2*pi*1i*theta); - psi = psi + 0.25*noise; + psi = psi + Params.nsf*noise; % --- Run --- [psi, Observ] = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ, nn); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m index 4d6226f..3b9d150 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m @@ -6,7 +6,7 @@ function [Params] = setupParameters(this) muB = CONSTANTS.BohrMagneton; % [J/T] a0 = CONSTANTS.BohrRadius; % [m] m0 = CONSTANTS.AtomicMassUnit; % [kg] - w0 = 2*pi*61.6582; % Angular frequency unit [s^-1] + w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] mu0factor = 0.3049584233607396; % =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const. % mu0=mu0factor *hbar^2*a0/(m0*muB^2) % Number of points in each direction @@ -46,6 +46,8 @@ function [Params] = setupParameters(this) % even though the solution is good, this just stops it going on forever Params.mindt = this.MinimumTimeStepSize; % Minimum size for a time step using adaptive dt + Params.nsf = this.NoiseScaleFactor; + Params.njob = this.JobNumber; @@ -70,7 +72,7 @@ function [Params] = setupParameters(this) Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2); % DDI strength - Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined + Params.gdd = 4*pi*Params.add/l0; % sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined % Trap gamma Params.gx = (Params.wx/w0)^2; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m index 2c3759e..3e38a4c 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m @@ -4,18 +4,18 @@ function [psi] = setupWavefunction(~,Params,Transf) Y = Transf.Y; ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0; - elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0; + elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0; - Rx = 4*sqrt(2)*ellx; - Ry = 4*sqrt(2)*elly; + Rx = 6*sqrt(2)*ellx; + Ry = 6*sqrt(2)*elly; X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; - psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2); + psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2); cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; - psi = psi/sqrt(cur_norm); + psi = psi/sqrt(cur_norm); + + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; + psi = sqrt(Params.N)*psi/sqrt(Norm); - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; - psi = sqrt(Params.N)*psi/sqrt(Norm); - end \ No newline at end of file