From e67f82b564745fc55bdd232bcb352e24c464c89f Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 17 Jun 2024 20:24:00 +0200 Subject: [PATCH] Expanded the Potential class. --- Dipolar Gas Simulator/+Scripts/test.m | 19 +- .../+Simulator/@Calculator/Calculator.m | 4 +- .../+Simulator/@DipolarGas/DipolarGas.m | 57 +++--- .../+Simulator/@DipolarGas/initialize.m | 9 +- .../@DipolarGas/propagateWavefunction.m | 170 ++++++++++++++++++ .../@DipolarGas/{runSimulation.m => run.m} | 12 +- .../+Simulator/@DipolarGas/setupParameters.m | 14 +- .../+Simulator/@DipolarGas/setupSpace.m | 2 +- .../+Simulator/@DipolarGas/solver.m | 89 --------- .../+Simulator/@Potentials/Potentials.m | 136 ++++++++++++-- .../+Simulator/@Potentials/setupParameters.m | 36 ++++ 11 files changed, 385 insertions(+), 163 deletions(-) create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m rename Dipolar Gas Simulator/+Simulator/@DipolarGas/{runSimulation.m => run.m} (71%) delete mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Potentials/setupParameters.m diff --git a/Dipolar Gas Simulator/+Scripts/test.m b/Dipolar Gas Simulator/+Scripts/test.m index 5fa5416..835e7a9 100644 --- a/Dipolar Gas Simulator/+Scripts/test.m +++ b/Dipolar Gas Simulator/+Scripts/test.m @@ -10,28 +10,31 @@ OptionsStruct.NumberOfAtoms = 1E6; OptionsStruct.DipolarPolarAngle = pi/2; OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.ScatteringLength = 86; + OptionsStruct.TrapFrequencies = [125, 125, 250]; +OptionsStruct.TrapDepth = 5; +OptionsStruct.BoxSize = 15; +OptionsStruct.TrapPotentialType = 'Harmonic'; + OptionsStruct.NumberOfGridPoints = [64, 64, 48]; OptionsStruct.Dimensions = [40, 40, 20]; OptionsStruct.CutoffType = 'Cylindrical'; -OptionsStruct.TrapPotentialType = 'Harmonic'; - OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStep = 50E-6; % in s -OptionsStruct.SimulationTime = 2E6; % in s +OptionsStruct.TimeStepSize = 50E-6; % in s +OptionsStruct.NumberOfTimeSteps = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Data'; - options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct -sim = Simulator.DipolarGas(options{:}); -pot = Simulator.Potentials(options{:}); +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); % + pot.repulsive_stirrer(); %-% Run Simulation %-% -[Params, Transf, psi, V, VDk] = sim.runSimulation(); +[Params, Transf, psi, V, VDk] = sim.run(); %% - Plot numerical grid Plotter.visualizeSpace(Transf) diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m index 93469f0..aae055e 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m @@ -42,8 +42,6 @@ classdef Calculator < handle & matlab.mixin.Copyable this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; this.CutoffType = p.Results.CutoffType; - this.CalculatorDefaults.CutoffType = this.CutoffType; - end function restoreDefaults(this) @@ -56,7 +54,7 @@ classdef Calculator < handle & matlab.mixin.Copyable this.CutoffType = this.CalculatorDefaults.CutoffType; end - end % - lifecycle + end % % methods diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m index 60b0855..e9a791e 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -8,14 +8,17 @@ classdef DipolarGas < handle & matlab.mixin.Copyable TrapFrequencies; NumberOfGridPoints; Dimensions; - + Potential; + SimulationMode; - TimeStep; - SimulationTime; + TimeStepSize; + NumberOfTimeSteps; EnergyTolerance; - MinimumTimeStep; + MinimumTimeStepSize; Calculator; + + SimulationParameters; %Flags @@ -47,16 +50,16 @@ classdef DipolarGas < handle & matlab.mixin.Copyable addParameter(p, 'Dimensions', 10 * ones(1,3),... @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',... - @(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))); + @(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})))); addParameter(p, 'CutoffType', 'Cylindrical',... - @(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))); - addParameter(p, 'TimeStep', 5E-4,... + @(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'})))); + addParameter(p, 'TimeStepSize', 5E-4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SimulationTime', 2e6,... + addParameter(p, 'NumberOfTimeSteps', 2e6,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'EnergyTolerance', 1e-10,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'MinimumTimeStep', 1e-6,... + addParameter(p, 'MinimumTimeStepSize', 1e-6,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'DebugMode', false,... @islogical); @@ -74,12 +77,12 @@ 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.SimulationMode = p.Results.SimulationMode; - this.TimeStep = p.Results.TimeStep; - this.SimulationTime = p.Results.SimulationTime; + this.TimeStepSize = p.Results.TimeStepSize; + this.NumberOfTimeSteps = p.Results.NumberOfTimeSteps; this.EnergyTolerance = p.Results.EnergyTolerance; - this.MinimumTimeStep = p.Results.MinimumTimeStep; + this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize; this.DebugMode = p.Results.DebugMode; this.DoSave = p.Results.SaveData; @@ -87,29 +90,25 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.Calculator = Simulator.Calculator('CutoffType', p.Results.CutoffType); - switch this.SimulationMode - case "ImaginaryTimeEvolution" - % Development In progress - case "RealTimeEvolution" - % Development In progress - end + this.SimulationParameters = this.setupParameters(); + end end methods - function set.TimeStep(this, val) + function set.TimeStepSize(this, val) assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!'); - this.TimeStep = val; + this.TimeStepSize = val; end - function ret = get.TimeStep(this) - ret = this.TimeStep; + function ret = get.TimeStepSize(this) + ret = this.TimeStepSize; end - function set.SimulationTime(this, val) -% assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!'); - this.SimulationTime = val; + function set.NumberOfTimeSteps(this, val) + assert(val <= 2E6, 'Not time efficient to compute for time spans longer than 2E6 seconds!'); + this.NumberOfTimeSteps = val; end - function ret = get.SimulationTime(this) - ret = this.SimulationTime; + function ret = get.NumberOfTimeSteps(this) + ret = this.NumberOfTimeSteps; end function set.NumberOfAtoms(this, val) assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!'); @@ -118,7 +117,6 @@ classdef DipolarGas < handle & matlab.mixin.Copyable function ret = get.NumberOfAtoms(this) ret = this.NumberOfAtoms; end - function set.DebugMode(this, val) this.DebugMode = val; end @@ -137,7 +135,6 @@ classdef DipolarGas < handle & matlab.mixin.Copyable function ret = get.SaveDirectory(this) ret = this.SaveDirectory; end - end % - setters and getters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m index ef860d5..1cc6a9d 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m @@ -1,9 +1,8 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad) - format long - X = Transf.X; Y = Transf.Y; Z = Transf.Z; - % == Trap potential == % - V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^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')) @@ -13,7 +12,7 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad) VDk = this.Calculator.calculateVDCutoff(Params,Transf,TransfRad); save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); end - fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', this.Calculator.CutoffType) + fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.\n', this.Calculator.CutoffType) % == Setting up the initial wavefunction == % psi = this.setupWavefunction(Params,Transf); diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m new file mode 100644 index 0000000..1cc2d12 --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m @@ -0,0 +1,170 @@ +function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ) + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + + switch this.SimulationMode + case 'ImaginaryTimeEvolution' + dt=-1j*abs(this.TimeStepSize); + + KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + Observ.residual = 1; Observ.res = 1; + + muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V); + AdaptIdx = 0; + + while t_idx < Params.sim_time_cut_off + %kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); + + %DDI + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + + %Real-space + psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem)); + + %kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); + + %Renorm + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; + psi = sqrt(Params.N)*psi/sqrt(Norm); + + muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V); + + if mod(t_idx,1000) == 0 + + %Change in Energy + E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; + + %Chemical potential + Observ.mucVec = [Observ.mucVec muchem]; + + %Normalized residuals + res = this.Calculator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem); + Observ.residual = [Observ.residual res]; + + Observ.res_idx = Observ.res_idx + 1; + + save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + + %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 AdaptIdx > 4 && abs(dt) > Params.mindt + dt = dt / 2; + fprintf('Time step changed to '); disp(dt); + AdaptIdx = 0; + elseif AdaptIdx > 4 && abs(dt) < Params.mindt + break + else + AdaptIdx = AdaptIdx + 1; + end + else + AdaptIdx = 0; + end + end + if any(isnan(psi(:))) + disp('NaNs encountered!') + break + end + t_idx=t_idx+1; + end + + %Change in Energy + E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; + + % Phase coherence + [PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params); + Observ.PCVec = [Observ.PCVec PhaseC]; + + Observ.res_idx = Observ.res_idx + 1; + save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + + case 'RealTimeEvolution' + + dt = abs(this.TimeStepSize); + + KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + + muchem = chemicalpotential(psi,Params,Transf,VDk,V); + + while t_idx < Params.sim_time_cut_off + % Parameters at time t + + %kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); + + %DDI + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + + %Real-space + psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem)); + + %kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); + + muchem = chemicalpotential(psi,Params,Transf,VDk,V); + + if mod(t_idx,1000)==0 + %Change in Normalization + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation + Observ.NormVec = [Observ.NormVec Norm]; + + %Change in Energy + E = energytotal(psi,Params,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; + + % Phase coherence + [PhaseC] = PhaseCoherence(psi,Transf); + Observ.PCVec = [Observ.PCVec PhaseC]; + + Observ.tVecPlot = [Observ.tVecPlot tVal]; + Observ.res_idx = Observ.res_idx + 1; + + save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',njob,Observ.res_idx),'psi','muchem','Observ','t_idx'); + end + if any(isnan(psi(:))) + disp('NaNs encountered!') + break + end + t_idx=t_idx+1; + end + + %Change in Normalization + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation + Observ.NormVec = [Observ.NormVec Norm]; + + %Change in Energy + E = energytotal(psi,Params,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; + + % Phase coherence + [PhaseC] = PhaseCoherence(psi,Transf); + Observ.PCVec = [Observ.PCVec PhaseC]; + + Observ.tVecPlot = [Observ.tVecPlot tVal]; + + Observ.res_idx = Observ.res_idx + 1; + save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',njob,Observ.res_idx),'psi','muchem','Observ','t_idx'); + + otherwise + disp('Choose a valid DDI cutoff type!') + return + end +end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m similarity index 71% rename from Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m rename to Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m index aba48d2..4db3dda 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m @@ -1,7 +1,8 @@ -function [Params, Transf, psi,V,VDk] = runSimulation(this) +function [Params, Transf, psi,V,VDk] = run(this) % --- Obtain simulation parameters --- [Params] = this.setupParameters(); - + this.SimulationParameters = Params; + % --- Set up spatial grids and transforms --- [Transf] = this.setupSpace(Params); [TransfRad] = this.setupSpaceRadial(Params); @@ -15,10 +16,7 @@ function [Params, Transf, psi,V,VDk] = runSimulation(this) t_idx = 1; %Start at t = 0; Observ.res_idx = 1; - % --- Job Settings --- - njob = 6; - mkdir(sprintf('./Data/Run_%03i',njob)) - % --- Run Simulation --- - % [psi] = this.solver(psi,Params,Transf,VDk,V,njob,t_idx,Observ); + % mkdir(sprintf('./Data/Run_%03i',Params.njob)) + % [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,njob,t_idx,Observ); end diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m index e0f516d..6d780ff 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -8,7 +8,7 @@ function [Params] = setupParameters(this) m0 = CONSTANTS.AtomicMassUnit; % [kg] w0 = 2*pi*100; % 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) + % mu0=mu0factor *hbar^2*a0/(m0*muB^2) % Number of points in each direction Params.Nx = this.NumberOfGridPoints(1); Params.Ny = this.NumberOfGridPoints(2); @@ -19,7 +19,7 @@ function [Params] = setupParameters(this) Params.Ly = this.Dimensions(2); Params.Lz = this.Dimensions(3); - % Masses + % Mass, length scale Params.m = CONSTANTS.Dy164Mass; l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length @@ -45,10 +45,12 @@ function [Params] = setupParameters(this) Params.gamma_S = 7.5*10^(-3); % gamma for the stochastic GPE Params.muchem = 12.64*Params.wz/w0; % fixing the chemical potential for the stochastic GPE - Params.Etol = this.EnergyTolerance; % Tolerances - Params.cut_off = this.SimulationTime; % sometimes the imaginary time gets a little stuck - % even though the solution is good, this just stops it going on forever - Params.mindt = this.MinimumTimeStep; % Minimum size for a time step using adaptive dt + Params.Etol = this.EnergyTolerance; % Tolerances + Params.sim_time_cut_off = this.NumberOfTimeSteps; % sometimes the imaginary time gets a little stuck + % 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.njob = 6; % ================ Parameters defined by those above ================ % diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m index 56996bd..bf58f64 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m @@ -1,4 +1,4 @@ -function [Transf] = setupSpace(this,Params) +function [Transf] = setupSpace(~,Params) Transf.Xmax = 0.5*Params.Lx; Transf.Ymax = 0.5*Params.Ly; Transf.Zmax = 0.5*Params.Lz; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m deleted file mode 100644 index 088f6b4..0000000 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m +++ /dev/null @@ -1,89 +0,0 @@ -function [psi] = propagate(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ) - set(0,'defaulttextInterpreter','latex') - set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); - - dt=-1j*abs(this.TimeStep); - - KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); - Observ.residual = 1; Observ.res = 1; - - muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V); - AdaptIdx = 0; - - while t_idx < Params.cut_off - %kin - psi = fftn(psi); - psi = psi.*exp(-0.5*1i*dt*KEop); - psi = ifftn(psi); - - %DDI - frho = fftn(abs(psi).^2); - Phi = real(ifftn(frho.*VDk)); - - %Real-space - psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem)); - - %kin - psi = fftn(psi); - psi = psi.*exp(-0.5*1i*dt*KEop); - psi = ifftn(psi); - - %Renorm - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; - psi = sqrt(Params.N)*psi/sqrt(Norm); - - muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V); - - if mod(t_idx,1000) == 0 - - %Change in Energy - E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; - Observ.EVec = [Observ.EVec E]; - - %Chemical potential - Observ.mucVec = [Observ.mucVec muchem]; - - %Normalized residuals - res = this.Calculator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem); - Observ.residual = [Observ.residual res]; - - Observ.res_idx = Observ.res_idx + 1; - - save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); - - %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 AdaptIdx > 4 && abs(dt) > Params.mindt - dt = dt / 2; - fprintf('Time step changed to '); disp(dt); - AdaptIdx = 0; - elseif AdaptIdx > 4 && abs(dt) < Params.mindt - break - else - AdaptIdx = AdaptIdx + 1; - end - else - AdaptIdx = 0; - end - end - if any(isnan(psi(:))) - disp('NaNs encountered!') - break - end - t_idx=t_idx+1; - end - - %Change in Energy - E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; - Observ.EVec = [Observ.EVec E]; - - % Phase coherence - [PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params); - Observ.PCVec = [Observ.PCVec PhaseC]; - - Observ.res_idx = Observ.res_idx + 1; - save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); -end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m b/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m index 4718c85..4920134 100644 --- a/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m +++ b/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m @@ -3,12 +3,19 @@ classdef Potentials < handle & matlab.mixin.Copyable properties (Access = private) PotentialDefaults = struct('TrapPotentialType', 'Harmonic', ... - 'TrapFrequency', 100e3); + 'TrapFrequencies', 100 * ones(1,3), ... + 'TrapDepth', 10, ... + 'BoxSize', 5); end properties (Access = public) TrapPotentialType; - TrapFrequency; + TrapFrequencies; + TrapDepth; + BoxSize; + NumberOfGridPoints; + Dimensions; + SimulationParameters; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -16,18 +23,107 @@ classdef Potentials < handle & matlab.mixin.Copyable methods function this = Potentials(varargin) - this.TrapFrequency = this.PotentialDefaults.TrapFrequency; - end - - function restoreDefaults(this) - this.TrapFrequency = this.PotentialDefaults.TrapFrequency; - end - - end % - lifecycle + + 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{:}); - % methods + 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 - % end % - setters and getters + 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); + z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz); + [X,Y,Z] = ndgrid(x,y,z); + ret = 0.0 * (X+Y+Z); + 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); + z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz); + [X,Y,Z] = ndgrid(x,y,z); + ret = 0.5*(this.SimulationParameters.gx.*X.^2+this.SimulationParameters.gy.*Y.^2+this.SimulationParameters.gz*Z.^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); + z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz); + [X,Y,Z] = ndgrid(x,y,z); + + heaviside_X = this.custom_heaviside(-abs(X) + this.SimulationParameters.boxsize/2); + heaviside_Y = this.custom_heaviside(-abs(Y) + this.SimulationParameters.boxsize/2); + heaviside_Z = this.custom_heaviside(-abs(Z) + this.SimulationParameters.boxsize/2); + + ret = this.SimulationParameters.boxdepth * (1 - heaviside_X .* heaviside_Y .* heaviside_Z); + 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); + z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz); + [X,Y,Z] = ndgrid(x,y,z); + + ret = this.SimulationParameters.boxdepth * (1 - this.custom_heaviside((this.SimulationParameters.boxsize/2)^2 - X.^2 - Y.^2 - Z.^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 @@ -50,7 +146,19 @@ classdef Potentials < handle & matlab.mixin.Copyable 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) @@ -63,4 +171,4 @@ classdef Potentials < handle & matlab.mixin.Copyable end end -end +end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Potentials/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@Potentials/setupParameters.m new file mode 100644 index 0000000..0c4797a --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@Potentials/setupParameters.m @@ -0,0 +1,36 @@ +function [Params] = setupParameters(this) + CONSTANTS = Helper.PhysicsConstants; + hbar = CONSTANTS.PlanckConstantReduced; % [J.s] + w0 = 2*pi*100; % 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); + Params.Nz = this.NumberOfGridPoints(3); + + % Dimensions (in units of l0) + Params.Lx = this.Dimensions(1); + Params.Ly = this.Dimensions(2); + Params.Lz = this.Dimensions(3); + + % 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/(Params.wz/w0); % 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