From ac47adcce33ba1eca8f508c84afda2c8f281112e Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sat, 30 Nov 2024 08:26:12 +0100 Subject: [PATCH] Adapted elements from a different working version of the solver from Wyatt - minor changes such as replacing trapz with sum and a few others to make the versions similar. --- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 12 ++++---- .../@Calculator/calculateChemicalPotential.m | 7 +++-- .../@Calculator/calculateEnergyComponents.m | 29 ------------------- .../calculateNormalizedResiduals.m | 2 +- .../@Calculator/calculateTotalEnergy.m | 11 +++---- .../@Calculator/calculateVariationalEnergy.m | 14 ++++----- .../@DipolarGas/DipolarGas.m | 2 +- .../@DipolarGas/propagateWavefunction.m | 10 +++---- .../+VariationalSolver2D/@DipolarGas/run.m | 23 ++++++++++----- .../@DipolarGas/setupParameters.m | 9 ++++-- .../@DipolarGas/setupWavefunction.m | 9 ++---- .../@Potentials/Potentials.m | 4 +-- 12 files changed, 57 insertions(+), 75 deletions(-) delete mode 100644 Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateEnergyComponents.m diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 2bc0f81..6e1d6ec 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -69,7 +69,7 @@ OptionsStruct.DipolarPolarAngle = 0; OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.ScatteringLength = 98.0676; -OptionsStruct.TrapFrequencies = [10, 10, 72.4]; +OptionsStruct.TrapFrequencies = [0, 0, 72.4]; OptionsStruct.TrapPotentialType = 'None'; OptionsStruct.NumberOfGridPoints = [128, 128]; @@ -78,14 +78,14 @@ OptionsStruct.TimeStepSize = 500E-6; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-04; +OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 4; -OptionsStruct.MaxIterations = 1; +OptionsStruct.MaxIterations = 10; OptionsStruct.VariationalWidth = 5.0; -OptionsStruct.WidthLowerBound = 0.2; +OptionsStruct.WidthLowerBound = 1; OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-2; +OptionsStruct.WidthCutoff = 5e-3; OptionsStruct.PlotLive = true; OptionsStruct.JobNumber = 1; @@ -101,7 +101,7 @@ solver.Potential = pot.trap(); %-% Run Solver %-% [Params, Transf, psi, V, VDk] = solver.run(); - +%% % Solve BdG equations % Load data Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V'); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m index 1bcd078..f8b14c9 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateChemicalPotential.m @@ -15,7 +15,7 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V) % Kinetic energy Ekin = KEop.*abs(fftn(psi)*normfac).^2; - Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; + Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; % Potential energy Epot = V.*abs(psi).^2; @@ -27,6 +27,7 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V) Eqf = gamma_eff*abs(psi).^5; % Total energy - muchem = Ekin + Ez*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; % - muchem = muchem / Params.N; %Only use if psi is normalized to N + muchem = Ekin + Ez*Params.N + sum(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; + muchem = muchem / Params.N; % Only use if psi is normalized to N + end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateEnergyComponents.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateEnergyComponents.m deleted file mode 100644 index 1011665..0000000 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateEnergyComponents.m +++ /dev/null @@ -1,29 +0,0 @@ -function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V,VParams) - - % Parameters - KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); - normfac = Params.Lx*Params.Ly/numel(psi); - - % DDIs - frho = fftn(abs(psi).^2); - Phi = (ifftn(frho.*VDk)); - Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);% - E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy; - - % Kinetic energy - Ekin = KEop.*abs(fftn(psi)*normfac).^2; - E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2 + 0.25*(1/VParams.ell^2)*Params.N; - - % Potential energy - Epot = V.*abs(psi).^2; - E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy + 0.25*Params.gz*VParams.ell^2*Params.N; - - % Contact interactions - Eint = 0.5*Params.gs*abs(psi).^4/(sqrt(2*pi)*VParams.ell); - E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy; - - % Quantum fluctuations - Eqf = 0.4*Params.gammaQF*abs(psi).^5*sqrt(2/(5*VParams.ell^3*pi^(3/2))); - E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; - -end diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m index 67e9523..2420f2e 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateNormalizedResiduals.m @@ -24,6 +24,6 @@ function res = calculateNormalizedResiduals(~,psi,Params,VParams,Transf,VDk,V,m Eqf = gamma_eff*abs(psi).^3.*psi; % Total energy - res = trapz(abs(Ekin(:) + Ez*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy); + res = sum(abs(Ekin(:) + Ez*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/sum(abs(muchem*psi(:))*Transf.dx*Transf.dy); end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m index 1fed471..3c1aa42 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateTotalEnergy.m @@ -12,23 +12,24 @@ function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V) frho = fftn(abs(psi).^2); Phi = real(ifftn(frho.*VDk)); Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell); - Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy; + Eddi = sum(Eddi(:))*Transf.dx*Transf.dy; % Kinetic energy Ekin = KEop.*abs(fftn(psi)*normfac).^2; - Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; + Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; % Potential energy Epot = V.*abs(psi).^2; - Epot = trapz(Epot(:))*Transf.dx*Transf.dy; + Epot = sum(Epot(:))*Transf.dx*Transf.dy; % Contact interactions Eint = 0.5*g_eff*abs(psi).^4; - Eint = trapz(Eint(:))*Transf.dx*Transf.dy; + Eint = sum(Eint(:))*Transf.dx*Transf.dy; % Quantum fluctuations Eqf = 0.4*gamma_eff*abs(psi).^5; - Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; + Eqf = sum(Eqf(:))*Transf.dx*Transf.dy; E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf; + end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m index 55f1698..4a780f3 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m @@ -7,29 +7,29 @@ function E = calculateVariationalEnergy(this, psi, Params, ell, Transf, V) % Parameters KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); normfac = Params.Lx*Params.Ly/numel(psi); - + % DDIs VDk = this.calculateVDkWithCutoff(Transf, Params, ell); % VDk depends on the variational parameters. THIS HAS TO BE RECALCULATED HERE FOR THE CONSTRAINED OPTIMIZATION TO WORK! frho = fftn(abs(psi).^2); Phi = real(ifftn(frho.*VDk)); - Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*ell);% + Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*ell); Eddi = sum(Eddi(:))*Transf.dx*Transf.dy; % Kinetic energy Ekin = KEop.*abs(fftn(psi)*normfac).^2; - Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; + Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; % Potential energy Epot = V.*abs(psi).^2; - Epot = trapz(Epot(:))*Transf.dx*Transf.dy; - + Epot = sum(Epot(:))*Transf.dx*Transf.dy; + % Contact interactions Eint = 0.5*g_eff*abs(psi).^4; - Eint = trapz(Eint(:))*Transf.dx*Transf.dy; + Eint = sum(Eint(:))*Transf.dx*Transf.dy; % Quantum fluctuations Eqf = 0.4*gamma_eff*abs(psi).^5; - Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; + Eqf = sum(Eqf(:))*Transf.dx*Transf.dy; % Total E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m index 01b80e8..c39a6ce 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/DipolarGas.m @@ -53,7 +53,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable addParameter(p, 'ScatteringLength', 120,... @(x) assert(isnumeric(x) && isscalar(x) && (x > -150) && (x < 150))); addParameter(p, 'TrapFrequencies', 100 * ones(1,3),... - @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); + @(x) assert(isnumeric(x) && isvector(x) && all(x >= 0))); addParameter(p, 'NumberOfGridPoints', 128 * ones(1,2),... @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); addParameter(p, 'Dimensions', 10 * ones(1,2),... diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m index da639d6..d3218f9 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m @@ -51,7 +51,7 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans psi = ifftn(psi); % Renorm - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy; psi = sqrt(Params.N)*psi/sqrt(Norm); muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); @@ -80,11 +80,11 @@ 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 relres < 1e-4 if AdaptIdx > 3 && abs(dt) > Params.mindt dt = dt / 2; AdaptIdx = 0; diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m index cc7feaf..416bbe4 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m @@ -23,29 +23,38 @@ function [Params, Transf, psi, V, VDk] = run(this) % --- Initialize --- mkdir(sprintf(this.SaveDirectory)) mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob)) - fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8); + fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-10,'Algorithm','sqp','MaxIterations',2000,'MaxFunctionEvaluations',4000); [psi,V,VDk] = this.initialize(Params,VParams,Transf); ells(1) = VParams.ell; E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V)/Params.N; E_vs_iter(1) = E_Var(ells(1)); + + t_idx = 1; firstpoint = 1; for nn = 1:Params.SelfConIter Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = []; Observ.res_idx = 1; - t_idx = 1; % Start at t = 0; - - [~,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter + if firstpoint == 1 + [psi,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter + firstpoint = 0; + else + [~,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter + end % --- Adding some noise --- % Noise added in every iteration to ensure it is not stuck in some local minimum - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation - psi = sqrt(Params.N)*psi/sqrt(Norm); + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation + psi = psi/sqrt(Norm); r = normrnd(0,1,size(psi)); theta = rand(size(psi)); noise = r.*exp(2*pi*1i*theta); psi = psi + Params.nsf*noise; - + + % Renormalize wavefunction + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation + psi = sqrt(Params.N)*psi/sqrt(Norm); + % --- Run --- [psi, Observ] = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ, nn); psi = gather(psi); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m index 021a5e2..0f63560 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m @@ -17,9 +17,8 @@ function [Params] = setupParameters(this) Params.Lx = this.Dimensions(1); Params.Ly = this.Dimensions(2); - % Mass, length scale + % Mass Params.m = CONSTANTS.Dy164Mass; - l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length % Atom numbers Params.N = this.NumberOfAtoms; @@ -65,12 +64,16 @@ function [Params] = setupParameters(this) Params.evarcutoff = this.VariationalEnergyTolerance; % ================ Parameters defined by those above ================ % - + % Length scale + l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length + % Contact interaction strength (units of l0/m) Params.gs = 4*pi*Params.as/l0; % Dipole lengths Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2); + Params.ppum2 = Params.N/(Params.Lx*Params.Ly*(l0*1e6)^2);% Particles per squared micron + Params.ppadd2 = Params.ppum2*(Params.add*1e6)^2; % Particles per squared add % DDI strength Params.gdd = 4*pi*Params.add/l0; % sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m index df440f4..f01321f 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupWavefunction.m @@ -8,16 +8,13 @@ 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 = 8*sqrt(2)*ellx; + Ry = 8*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); - cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; - psi = psi/sqrt(cur_norm); - - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy; psi = sqrt(Params.N)*psi/sqrt(Norm); end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m index 3006d9b..99cd942 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m @@ -26,7 +26,7 @@ classdef Potentials < handle & matlab.mixin.Copyable 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))); + @(x) assert(isnumeric(x) && isvector(x) && all(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),... @@ -68,7 +68,7 @@ classdef Potentials < handle & matlab.mixin.Copyable methods function set.TrapFrequencies(this, val) - assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined trap frequencies!'); + assert(isnumeric(val) && isvector(val) && all(val >= 0), 'Incorrectly defined trap frequencies!'); this.TrapFrequencies = val; end function ret = get.TrapFrequencies(this)