diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 34e5aef..0929700 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/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) diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m index 8723d20..20170dd 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m @@ -1,6 +1,6 @@ -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 QX = Transf.KX*VParams.ell/sqrt(2); QY = Transf.KY*VParams.ell/sqrt(2); @@ -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); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m index 00a4da6..b2446c7 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m +++ b/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; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m index d958c62..75a536b 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m @@ -1,18 +1,15 @@ function [psi,V,VDk] = initialize(this,Params,VParams,Transf) - format long - X = Transf.X; Y = Transf.Y; - rcut = min(Transf.Xmax,Transf.Ymax); - - % == Potential == % - V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2); + % == User-defined potential == % + V = this.Potential; + assert(~anynan(V), 'Potential not defined! Specify as .Potential = .trap() + .'); % == 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 diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m index 06ad9ff..bca3b99 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m @@ -25,6 +25,9 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk, g_eff = Params.gs * (1/(sqrt(2*pi)*VParams.ell)); 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 @@ -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-4 - if AdaptIdx > 3 && abs(dt) > Params.mindt + if relres <1e-5 + if AdaptIdx > 4 && 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 diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m index 3e488a9..2c3759e 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m +++ b/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; - Ry = 6*sqrt(2)*elly; + Rx = 4*sqrt(2)*ellx; + Ry = 4*sqrt(2)*elly; X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m new file mode 100644 index 0000000..ae4b150 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/setupParameters.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/setupParameters.m new file mode 100644 index 0000000..29cbc26 --- /dev/null +++ b/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 \ No newline at end of file