Added Potentials class, removed in-plane trap

Dipolar-Gas-Simulator/+Scripts/run_locally.m

@@ -58,6 +58,7 @@ OptionsStruct.DipolarAzimuthAngle     = 0;
 OptionsStruct.ScatteringLength        = 98.0676;     % 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
@@ -74,9 +75,11 @@ options                               = Helper.convertstruct2cell(OptionsStruct)
 clear OptionsStruct
 
 solver                                = VariationalSolver2D.DipolarGas(options{:});
+pot                                   = VariationalSolver2D.Potentials(options{:});
+solver.Potential                      = pot.trap();
 
 %-% Run Solver %-%
-[Params, Transf, psi, V, VDk] = solver.run();
+[Params, Transf, psi, V, VDk]         = solver.run();
 
 %% - Plot numerical grid
 % Plotter.visualizeSpace2D(Transf)

Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m

@@ -1,4 +1,4 @@
-function VDk = calculateVDkWithCutoff(~, rcut, Transf, Params, VParams)
+function VDk = calculateVDkWithCutoff(~, Transf, Params, VParams)
 % == Calculating the DDI potential in Fourier space with appropriate cutoff == %
     
     % Interaction in K space
@@ -21,7 +21,7 @@ function VDk = calculateVDkWithCutoff(~, rcut, Transf, Params, VParams)
     % Implementing a cutoff:
     VDr             = ifftn(VDk);
     VDr             = fftshift(VDr);
-    
+    rcut            = min(Transf.Xmax,Transf.Ymax);
     R               = sqrt(Transf.X.^2+Transf.Y.^2) < 0.9*rcut;
     VDr             = VDr.*double(R);
     VDr             = ifftshift(VDr);

Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m

@@ -8,6 +8,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
         TrapFrequencies;
         NumberOfGridPoints;
         Dimensions;
+        Potential;
 
         SimulationMode; 
         TimeStepSize;
@@ -79,6 +80,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
             this.TrapFrequencies      = p.Results.TrapFrequencies;
             this.NumberOfGridPoints   = p.Results.NumberOfGridPoints;
             this.Dimensions           = p.Results.Dimensions;
+            this.Potential            = NaN;
             this.TimeStepSize         = p.Results.TimeStepSize;
             this.MinimumTimeStepSize  = p.Results.MinimumTimeStepSize;
             this.TimeCutOff           = p.Results.TimeCutOff;

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

@@ -1,18 +1,15 @@
 function [psi,V,VDk] = initialize(this,Params,VParams,Transf)
     
-    format long
+    % == User-defined potential == %
-    X = Transf.X; Y = Transf.Y;
+    V = this.Potential;
-    rcut = min(Transf.Xmax,Transf.Ymax);
+    assert(~anynan(V), 'Potential not defined! Specify as <SimulatorObject>.Potential = <PotentialsObject>.trap() + <AdditionalTerms>.');
-
-    % == Potential == %
-    V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2);
     
     % == Calculating the DDIs == %
     if isfile(strcat(this.SaveDirectory, '/VDk_M.mat'))
         VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')));
         VDk = VDk.VDk;
     else
-        VDk = this.Calculator.calculateVDkWithCutoff(rcut, Transf, Params, VParams);
+        VDk = this.Calculator.calculateVDkWithCutoff(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

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

@@ -26,6 +26,9 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
     gamma_eff       = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
     Ez              = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
     
+    pb = Helper.ProgressBar();
+    pb.run('Running simulation in imaginary time: ');
+
     while t_idx < Params.sim_time_cut_off
 
         % kin
@@ -68,12 +71,11 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
     
             %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-5
-            if relres <1e-4
+                if AdaptIdx > 4 && abs(dt) > Params.mindt
-                if AdaptIdx > 3 && abs(dt) > Params.mindt
                     dt = dt / 2;
                     AdaptIdx = 0;
-                elseif AdaptIdx > 3 && abs(dt) < Params.mindt
+                elseif AdaptIdx > 4 && abs(dt) < Params.mindt
                     break
                 elseif Observ.residual(Observ.res_idx) < Params.rtol
                     break
@@ -84,11 +86,14 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
                 AdaptIdx = 0;
             end
         end
+
         if any(isnan(psi(:)))
             disp('NaNs encountered!')
             break
         end
+        
         t_idx=t_idx+1;
+        pb.run(100*t_idx/Params.sim_time_cut_off);
     end
     
     % Change in Energy

Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m

@@ -6,8 +6,8 @@ function [psi] = setupWavefunction(~,Params,Transf)
     ellx     = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0; 
     elly     = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
     
-    Rx       = 6*sqrt(2)*ellx; 
+    Rx       = 4*sqrt(2)*ellx; 
-    Ry       = 6*sqrt(2)*elly;
+    Ry       = 4*sqrt(2)*elly;
     X0       = 0.0*Transf.Xmax; 
     Y0       = 0.0*Transf.Ymax;
     

Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m

@@ -0,0 +1,171 @@
+classdef Potentials < handle & matlab.mixin.Copyable
+
+    properties (Access = private)
+        
+        PotentialDefaults     = struct('TrapPotentialType',           'Harmonic', ...
+                                       'TrapFrequencies',        100 * ones(1,3), ...
+                                       'TrapDepth',                           10, ...
+                                       'BoxSize',                              5);
+    end
+
+    properties (Access = public)
+        TrapPotentialType;
+        TrapFrequencies;
+        TrapDepth;
+        BoxSize;
+        NumberOfGridPoints;
+        Dimensions;
+        SimulationParameters;
+    end
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %- Methods
+    
+    methods
+        function this  = Potentials(varargin)
+            
+            p = inputParser;
+            p.KeepUnmatched = true;
+            
+            addParameter(p, 'TrapPotentialType', this.PotentialDefaults.TrapPotentialType,... 
+                @(x) assert(any(strcmpi(x,{'None','Harmonic','SquareBox','RoundBox'}))));
+            addParameter(p, 'TrapFrequencies',   this.PotentialDefaults.TrapFrequencies,...
+                @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
+            addParameter(p, 'TrapDepth',     this.PotentialDefaults.TrapDepth,...
+                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
+            addParameter(p, 'BoxSize', this.PotentialDefaults.BoxSize,...
+                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
+            addParameter(p, 'NumberOfGridPoints',    128 * ones(1,3),...
+                @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
+            addParameter(p, 'Dimensions',        10 * ones(1,3),...
+                @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
+            
+            p.parse(varargin{:});
+
+            this.TrapPotentialType    = p.Results.TrapPotentialType;
+            this.TrapFrequencies      = p.Results.TrapFrequencies;
+            this.TrapDepth            = p.Results.TrapDepth;
+            this.BoxSize              = p.Results.BoxSize;
+            this.NumberOfGridPoints   = p.Results.NumberOfGridPoints;
+            this.Dimensions           = p.Results.Dimensions;
+
+            this.SimulationParameters = this.setupParameters();
+        end
+        
+        function [ret] = trap(this)
+            format long
+            switch this.TrapPotentialType
+                case 'None'
+                    x        = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
+                    y        = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
+                    [X,Y]    = ndgrid(x,y);
+
+                    ret      = 0.0 * (X+Y);
+                case 'Harmonic'
+                    x        = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
+                    y        = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
+                    [X,Y]    = ndgrid(x,y);
+
+                    ret      = 0.5*(this.SimulationParameters.gx.*X.^2+this.SimulationParameters.gy.*Y.^2); 
+                case 'SquareBox'
+                    x        = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
+                    y        = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
+                    [X,Y]    = ndgrid(x,y);
+                    
+                    heaviside_X = this.custom_heaviside(-abs(X) + this.SimulationParameters.boxsize/2);
+                    heaviside_Y = this.custom_heaviside(-abs(Y) + this.SimulationParameters.boxsize/2);
+                    
+                    ret = this.SimulationParameters.boxdepth * (1 - heaviside_X .* heaviside_Y);
+                case 'RoundBox'
+                    x        = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
+                    y        = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
+                    [X,Y]    = ndgrid(x,y);
+                    
+                    ret = this.SimulationParameters.boxdepth * (1 - this.custom_heaviside((this.SimulationParameters.boxsize/2)^2 - X.^2 - Y.^2));
+            end
+        end
+    
+        function restoreDefaults(this)
+            this.TrapPotentialType          = this.PotentialDefaults.TrapPotentialType;
+            this.TrapFrequencies            = this.PotentialDefaults.TrapFrequencies;
+        end
+    end %
+
+    methods
+        function set.TrapFrequencies(this, val)
+            assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined trap frequencies!');
+            this.TrapFrequencies = val;
+        end
+        function ret = get.TrapFrequencies(this)
+            ret = this.TrapFrequencies;
+        end
+        function set.TrapPotentialType(this, val)
+            assert(any(strcmpi(val,{'None','Harmonic','SquareBox','RoundBox'})), 'Trap potential of the specified type cannot be generated!');
+            this.TrapPotentialType = val;
+        end
+        function ret = get.TrapPotentialType(this)
+            ret = this.TrapPotentialType;
+        end
+        function set.NumberOfGridPoints(this, val)
+            assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined grid!');
+            this.NumberOfGridPoints = val;
+        end
+        function ret = get.NumberOfGridPoints(this)
+            ret = this.NumberOfGridPoints;
+        end
+        function set.Dimensions(this, val)
+            assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined dimensions!');
+            this.Dimensions = val;
+        end
+        function ret = get.Dimensions(this)
+            ret = this.Dimensions;
+        end
+    end % - setters and getters
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %- Methods
+    
+    methods(Access = protected)
+        function cp = copyElement(this)
+            % Shallow copy object
+            cp = copyElement@matlab.mixin.Copyable(this);
+            
+            % Forces the setter to redefine the function handles to the new copied object
+            
+            pl = properties(this);
+            for k = 1:length(pl)
+                sc = superclasses(this.(pl{k}));
+                if any(contains(sc,{'matlab.mixin.Copyable'}))
+                    cp.(pl{k}) = this.(pl{k}).copy();
+                end
+            end
+        end
+    end
+    
+    methods (Static)
+
+        function y = custom_heaviside(x)
+            % This function computes the Heaviside step function with a custom value for Heaviside(0).
+            % x: input array
+            % H0: value for Heaviside(0)
+            
+            % Use MATLAB's built-in heaviside function
+            y = heaviside(x);
+            
+            % Replace the default value for Heaviside(0) with the custom value H0
+            y(x == 0) = 1;
+        end
+
+        % Creates an Instance of Class, ensures singleton behaviour (that there
+        % can only be one Instance of this class
+        function singleObj = getInstance(varargin)
+            % Creates an Instance of Class, ensures singleton behaviour
+            persistent localObj;
+            if isempty(localObj) || ~isvalid(localObj)
+                localObj =  Simulator.Potentials(varargin{:});
+            end
+            singleObj = localObj;
+        end
+    end
+    
+end

Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/setupParameters.m

@@ -0,0 +1,34 @@
+function [Params] = setupParameters(this)
+    CONSTANTS      = Helper.PhysicsConstants;
+    hbar           = CONSTANTS.PlanckConstantReduced; % [J.s]
+    w0             = 2*pi*61.6582;                    % Angular frequency unit [s^-1]
+    
+    % Mass, length scale
+    Params.m       = CONSTANTS.Dy164Mass;
+    l0             = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
+
+    % Number of points in each direction
+    Params.Nx      = this.NumberOfGridPoints(1);
+    Params.Ny      = this.NumberOfGridPoints(2);
+    
+    % Dimensions (in units of l0)
+    Params.Lx      = this.Dimensions(1);
+    Params.Ly      = this.Dimensions(2);
+    
+    % Trapping frequencies
+    Params.wx      = 2*pi*this.TrapFrequencies(1);
+    Params.wy      = 2*pi*this.TrapFrequencies(2);
+    Params.wz      = 2*pi*this.TrapFrequencies(3);
+    
+    % Trap depth and box size for box potentials
+    Params.boxdepth = this.TrapDepth;  % The depth of the box
+    Params.boxsize  = this.BoxSize;         % The size of the box
+
+    % ================ Parameters defined by those above ================ %
+    
+    % Trap gamma
+    Params.gx      = (Params.wx/w0)^2;
+    Params.gy      = (Params.wy/w0)^2;
+    Params.gz      = (Params.wz/w0)^2;
+
+end