diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index f88b3c0..7eb96d0 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -62,16 +62,16 @@ 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 = 1E4; % in s +OptionsStruct.TimeStepSize = 5E-6; % in s +OptionsStruct.TimeCutOff = 1E5; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.MaxIterations = 20; OptionsStruct.VariationalWidth = 4; -OptionsStruct.WidthLowerBound = 5; -OptionsStruct.WidthUpperBound = 8; -OptionsStruct.WidthCutoff = 1e-3; +OptionsStruct.WidthLowerBound = 0.2; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 1e-4; OptionsStruct.PlotLive = true; OptionsStruct.JobNumber = 1; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m index 155c71a..2b08ccb 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m @@ -2,7 +2,7 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V) 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); + Ez = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2); % Parameters normfac = Params.Lx*Params.Ly/numel(psi); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m index e4eb9d2..67e9523 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m @@ -2,7 +2,7 @@ function res = calculateNormalizedResiduals(~,psi,Params,VParams,Transf,VDk,V,m 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); + Ez = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2); KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m index 43a91e7..0d9823a 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m @@ -2,7 +2,7 @@ function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V) 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); + Ez = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2); % Parameters KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m index 20170dd..1f98cbc 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m @@ -10,10 +10,9 @@ function VDk = calculateVDkWithCutoff(~, Transf, Params, VParams) QDsq = QX.^2*cos(Params.phi)^2 + QY.^2*sin(Params.phi)^2; % Bare interaction - Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; + Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x) Fperp = 2 - 3*sqrt(pi).*absQ.*erfcx(absQ); Fpar(absQ == 0) = -1; - % Fpar(isinf(Fpar)) = -1; % Full DDI VDk = (Fpar*sin(Params.theta)^2 + Fperp*cos(Params.theta)^2); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m index b599f9f..f3466ae 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m @@ -2,7 +2,7 @@ function E = calculateVariationalEnergy(~, psi, Params, ell, Transf, VDk, V) g_eff = Params.gs * (1/(sqrt(2*pi)*ell)); gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*ell^(3/2))); - Ez = (0.25*ell^2) + (0.25*Params.gz*ell^2); + Ez = (0.25/ell^2) + (0.25*Params.gz*ell^2); % Parameters KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m index d02fcb8..2924d77 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m @@ -29,7 +29,7 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans 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); + Ez = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2); pb = Helper.ProgressBar(); fprintf('\n') @@ -61,7 +61,7 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); %Plotting loop - if mod(t_idx,250) == 0 + if mod(t_idx,500) == 0 % Change in Energy E = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V); @@ -85,11 +85,14 @@ 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); - if relres <1e-5 - if AdaptIdx > 4 && abs(dt) > Params.mindt + 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; AdaptIdx = 0; - elseif AdaptIdx > 4 && abs(dt) < Params.mindt + elseif AdaptIdx > 3 && abs(dt) < Params.mindt break elseif Observ.residual(Observ.res_idx) < Params.rtol break