add tangled

python/energy_flow_proper/11_new_ho_comparison/zero_t.py

@@ -0,0 +1,182 @@
+import figsaver as fs
+import plot_utils as pu
+
+import hops.util.bcf
+import hops.util.bcf_fits
+from hops.core.hierarchy_parameters import *
+from stocproc import StocProc_FFT
+from hops.util.abstract_truncation_scheme import TruncationScheme_Simplex
+from hops.core.integration import HOPSSupervisor
+import qutip
+
+
+def ho_zero(max_HO_level=15, wc=2, s=1, bcf_terms=7, k_max=5, bcf_scale=0.2):
+    hops_bcf = hops.util.bcf.OhmicBCF_zeroTemp(
+        s,
+        1,
+        wc,
+        normed=True
+    )
+
+    g, w = hops_bcf.exponential_coefficients(bcf_terms)
+    integration = IntP(t_steps=(40, 100))
+
+    q = qutip.operators.create(max_HO_level) + qutip.operators.destroy(max_HO_level)
+    p = (
+        qutip.operators.destroy(max_HO_level) - qutip.operators.create(max_HO_level)
+    ) / 1j
+
+    system = SysP(
+        H_sys=0.25 * (p ** 2 + q ** 2).full(),
+        L=[0.5 * q.full()],
+        psi0=qutip.states.fock(max_HO_level, n=1).full().flatten(),
+        g=[g],
+        w=[w],
+        bcf_scale=[bcf_scale],
+        T=[0],
+        description="A single HO at zero Temperature",
+    )
+
+
+    params = HIParams(
+        SysP=system,
+        IntP=integration,
+        HiP=HiP(
+            seed=0,
+            nonlinear=True,
+            result_type=ResultType.ZEROTH_AND_FIRST_ORDER,
+            truncation_scheme=TruncationScheme_Simplex(k_max),
+            save_therm_rng_seed=True,
+            auto_normalize=True,
+            accum_only=False,
+            rand_skip=None,
+        ),
+        Eta=[
+            StocProc_FFT(
+                spectral_density=hops.util.bcf.OhmicSD_zeroTemp(
+                    s,
+                    1,
+                    wc,
+                    normed=True
+                ),
+                alpha=hops_bcf,
+                t_max=integration.t_max,
+                intgr_tol=1e-5,
+                intpl_tol=1e-5,
+                negative_frequencies=False,
+            )
+        ],
+        EtaTherm=[None],
+    )
+
+    return params
+
+model_keys = [dict(bcf_scale=1, wc=ω_c) for ω_c in [1, 2, 3]] + [
+    dict(bcf_scale=bcf_scale, wc=2) for bcf_scale in [1, 2]
+]
+
+multi_params = [ho_zero(**keys) for keys in model_keys]
+multi_params = multi_params
+
+import ray
+ray.shutdown()
+ray.init(address="auto")
+
+supervisors = []
+for params in multi_params:
+      supervisor = HOPSSupervisor(
+            params,
+            5000,
+            data_path="ho_data",
+            data_name="zero_t",
+        )
+      supervisor.integrate()
+      supervisors.append(supervisor)
+
+from hopsflow import hopsflow
+
+flow_hops = []
+for supervisor in supervisors:
+    hf_sys = hopsflow.SystemParams.from_hi_params(supervisor.params)
+    flow_hops.append(
+        hopsflow.heat_flow_from_data(
+            supervisor.get_data(True),
+            hf_sys,
+            every=1000,
+            save=f"flow_zero",
+        )
+    )
+
+fig, ax = plt.subplots()
+for params, flow, keys in zip(multi_params, flow_hops, model_keys):
+    pu.plot_with_σ(
+        params.IntP.t,
+        -1 * flow.for_bath(0),
+        ax=ax,
+        label=f"{len(params.SysP.g[0])}, {params.HiP.truncation_scheme.kmax}",
+    )
+ax.legend()
+
+from hopsflow import gaussflow as gf
+
+exact_flows = []
+for params in multi_params:
+    α_0 = gf.BCF(
+        params.IntP.t_max,
+        factors=params.SysP.g[0],
+        exponents=params.SysP.w[0],
+        # hops_bcf,
+        # num_terms=6,
+        # resolution=0.001,
+    )
+
+    α_0_dot = gf.BCF(
+        params.IntP.t_max,
+        factors=-params.SysP.w[0] * params.SysP.g[0],
+        exponents=params.SysP.w[0],
+        # lambda t: 2 / (1j * np.pi) * (wc / (1 + 1j * wc * t)) ** 3,
+        # num_terms=6,
+        # resolution=0.001,
+    )
+
+    sys = gf.SystemParams(Ω=1, η=1, α_0=α_0)
+    flow = gf.Flow(sys, α_0, α_0_dot, n=1)
+    flow_τ = flow(params.IntP.t)
+    exact_flows.append(flow_τ)
+
+fig, ax = plt.subplots()
+for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
+    consistency = (-1 * flow).consistency(ex_flow)
+    pu.plot_with_σ(
+        params.IntP.t,
+        -1 * flow,
+        bath=0,
+        ax=ax,
+        label=rf"$α(0)={params.SysP.g[0].sum().real:.2f}$ $ω_c={keys['wc']}$ ${consistency}\%$",
+    )
+    ax.plot(params.IntP.t, ex_flow, linestyle="dotted", color="black")
+
+    plt.xlabel("$τ$")
+    plt.ylabel("$-J$")
+
+# ax.plot(multi_params[-1].IntP.t, flow_τ, label="Analytic", linestyle="dotted", color="black")
+ax.legend()
+# fs.export_fig("convergence_zero")
+
+pu.plot_convergence(
+    multi_params[-1].IntP.t,
+    -1 * flow_hops[-1].for_bath(0)
+    -flow_τ,
+    reference=np.zeros_like(flow_τ)
+)
+plt.axhline(0)
+plt.legend()
+# plt.plot(multi_params[-1].IntP.t, -1 * flow_hops[-1].for_bath(0).σ)
+
+f, ax = pu.plot_σ_development(flow_hops[-1].for_bath(0), label="Numerics")
+smooth_n = np.linspace(flow_hops[-1].Ns[0], flow_hops[-1].Ns[-1], 1000)
+ax.plot(smooth_n, flow_hops[-1][0].σ.mean() * np.sqrt(flow_hops[-1][0].N)/np.sqrt(smooth_n), label=r"$1/\sqrt{N}$")
+ax.set_yscale("log")
+ax.set_xscale("log")
+ax.legend()
+fs.export_fig("sqrt_convergence", f)