From 53894fb77ba94f95a0287d8c9d896e4bcfecfe23 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 14 Nov 2024 11:13:07 +0100 Subject: [PATCH] Fixed bugs to have a working version of the Variational2D solver. --- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 39 +++++++++++++++++++ .../@Calculator/calculateVDkWithCutoff.m | 13 ++++--- .../@Calculator/calculateVariationalEnergy.m | 5 ++- .../+Variational2D/@Calculator/find_nk_fit.m | 6 ++- .../+Variational2D/@DipolarGas/DipolarGas.m | 2 +- .../+Variational2D/@DipolarGas/initialize.m | 2 +- .../@DipolarGas/propagateWavefunction.m | 20 +++++----- .../+Variational2D/@DipolarGas/run.m | 11 +++--- 8 files changed, 71 insertions(+), 27 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index cd55b4b..9b5252e 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -39,6 +39,45 @@ sim.Potential = pot.trap(); % + pot.repulsive_chopstick(); %-% Run Simulation %-% [Params, Transf, psi, V, VDk] = sim.run(); +%% - Plot numerical grid +Plotter.visualizeSpace(Transf) +%% - Plot trap potential +Plotter.visualizeTrapPotential(sim.Potential,Params,Transf) +%% - Plot initial wavefunction +Plotter.visualizeWavefunction(psi,Params,Transf) +%% - Plot GS wavefunction +Plotter.visualizeGSWavefunction(Params.njob) + +%% - Create Variational2D and Calculator object with specified options + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 1.24E5; +OptionsStruct.DipolarPolarAngle = 0; +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 90; + +OptionsStruct.TrapFrequencies = [44.97, 10.4, 126.3]; + +OptionsStruct.NumberOfGridPoints = [128, 256, 128]; +OptionsStruct.Dimensions = [50, 120, 150]; +OptionsStruct.TimeStepSize = 500E-6; % in s +OptionsStruct.NumberOfTimeSteps = 100; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; + +OptionsStruct.JobNumber = 1; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Data'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = Variational2D.DipolarGas(options{:}); + +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + %% - Plot numerical grid Plotter.visualizeSpace(Transf) %% - Plot trap potential diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m index e7573ad..01eb82a 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m @@ -7,17 +7,18 @@ function VDk = calculateVDkWithCutoff(~, rcut, Transf, Params, VParams) Qsq = QX.^2 + QY.^2; absQ = sqrt(Qsq); - QDsq = QX.^2*cos(Params.eta)^2 + QY.^2*sin(Params.eta)^2; + QDsq = QX.^2*cos(Params.phi)^2 + QY.^2*sin(Params.phi)^2; - %Bare interaction + % Bare interaction Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; Fperp = 2 - 3*sqrt(pi).*absQ.*erfcx(absQ); - Fpar(absQ == 0) = -1; %Fpar(isinf(Fpar)) = -1; + Fpar(absQ == 0) = -1; + % Fpar(isinf(Fpar)) = -1; - %Full DDI - VDk = (Fpar*sin(Params.alpha)^2 + Fperp*cos(Params.alpha)^2); + % Full DDI + VDk = (Fpar*sin(Params.theta)^2 + Fperp*cos(Params.theta)^2); - %Implementing a cutoff: + % Implementing a cutoff: VDr = ifftn(VDk); VDr = fftshift(VDr); diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m index fa7f170..ba52d75 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m @@ -1,4 +1,4 @@ -function E = calculateVariationalEnergy(~, psi, Params, VarArray, Transf, VDk, V) +function E = calculateVariationalEnergy(this, psi, Params, VarArray, Transf, VDk, V) VParams.ell = VarArray(1); VParams.nu = VarArray(2); @@ -12,7 +12,7 @@ function E = calculateVariationalEnergy(~, psi, Params, VarArray, Transf, VDk, V normfac = Params.Lx*Params.Ly/numel(psi); % DDIs - [VParams] = find_nk_fit(VParams); % Not totally sure this should be here + [VParams] = this.find_nk_fit(VParams); % Not totally sure this should be here frho = fftn(abs(psi).^2); Phi = real(ifftn(frho.*VDk)); Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell_eff);% @@ -34,5 +34,6 @@ function E = calculateVariationalEnergy(~, psi, Params, VarArray, Transf, VDk, V Eqf = 0.4*gamma_eff*abs(psi).^5; Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; + % Total E = EVar*Params.N + Ekin + Epot + Eint + Eddi + Eqf; end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m index 086c94d..01286ee 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m @@ -1,4 +1,4 @@ -function [VParams] = find_nk_fit(VParams) +function [ret] = find_nk_fit(~,VParams) Lz = 20; Nz = 4000; @@ -20,5 +20,7 @@ function [VParams] = find_nk_fit(VParams) f = fit(kz',nk_forfit','gauss1'); % f(x) = a1*exp(-((x-b1)/c1)^2) - VParams.ell_eff = 1/f.c1; + ret.ell = VParams.ell; + ret.nu = VParams.nu; + ret.ell_eff = 1/f.c1; end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m index 3224d42..5cd658d 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m @@ -91,7 +91,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.DoSave = p.Results.SaveData; this.SaveDirectory = p.Results.SaveDirectory; - this.Calculator = Simulator.Calculator(); + this.Calculator = Variational2D.Calculator(); this.SimulationParameters = this.setupParameters(); diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m index 5abe950..519f053 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m @@ -15,7 +15,7 @@ function [psi,V,VDk] = initialize(this,Params,VParams,Transf) VDk = this.Calculator.calculateVDkWithCutoff(rcut, Transf, Params, VParams); save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); end - fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.\n', this.Calculator.CutoffType) + fprintf('Computed and saved DDI potential in Fourier space with cutoff.\n') % == Setting up the initial wavefunction == % psi = this.setupWavefunction(Params,Transf); diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m index 01f68b4..47f1fd5 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m @@ -2,21 +2,21 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk, set(0,'defaulttextInterpreter','latex') set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); - dt =-1j*abs(this.TimeStepSize); % Imaginary Time + dt =-1j*abs(this.TimeStepSize); % Imaginary Time - KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); - Observ.res = 1; - AdaptIdx = 0; + Observ.res = 1; + AdaptIdx = 0; % Change in Energy - E = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V); - E = E/Params.N; - Observ.EVec = [Observ.EVec E]; + E = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V); + E = E/Params.N; + Observ.EVec = [Observ.EVec E]; % Chemical Potential - muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); - Observ.mucVec = [Observ.mucVec muchem]; + muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); + Observ.mucVec = [Observ.mucVec muchem]; % Normalized residuals res = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem); @@ -107,7 +107,7 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk, Observ.residual = [Observ.residual res]; Observ.res_idx = Observ.res_idx + 1; - %Chemical potential + % Chemical potential Observ.mucVec = [Observ.mucVec muchem]; pb.run(' - Job Completed!'); diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m index 66e9017..7876214 100644 --- a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m @@ -12,7 +12,7 @@ function [Params, Transf, psi, V, VDk] = run(this) VarArray = Params.VarArray; VParams.ell = VarArray(1); VParams.nu = VarArray(2); - [VParams] = find_nk_fit(VParams); + [VParams] = this.Calculator.find_nk_fit(VParams); % --- Initialize --- mkdir(sprintf(this.SaveDirectory)) @@ -20,23 +20,23 @@ function [Params, Transf, psi, V, VDk] = run(this) fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8); - [psi,V,VDk,~] = this.initialize(Params,VParams,Transf); + [psi,V,VDk] = this.initialize(Params,VParams,Transf); ells(1) = VarArray(1); nus(1) = VarArray(2); E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi,Params,x,Transf,VDk,V)/Params.N; E_vs_iter(1) = E_Var([ells(1) nus(1)]); t_idx = 1; % Start at t = 0; - + for nn = 1:Params.SelfConIter Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.res_idx = 1; VParams.ell = VarArray(1); VParams.nu = VarArray(2); - [VParams] = find_nk_fit(VParams); + [VParams] = this.Calculator.find_nk_fit(VParams); - [psi,V,VDk,~] = this.initialize(Params,VParams,Transf); + [psi,V,VDk] = this.initialize(Params,VParams,Transf); % --- Adding some noise --- Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation @@ -69,4 +69,5 @@ function [Params, Transf, psi, V, VDk] = run(this) end save(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','Observ','Transf','Params','VDk','V','VarArray'); delete(sprintf('./Data/Run_%03i/psi_gs_*.mat',runIdx)) + end