add energy shovel motor

.gitignore

@@ -503,3 +503,5 @@ results/**.npy
 /hops_data/08/
 /python/energy_flow_proper/11_new_ho_comparison/ho_data/
 /python/energy_flow_proper/11_new_ho_comparison/taurus/
+/python/energy_flow_proper/11_new_ho_comparison/ho_data_local/
+/python/energy_flow_proper/11_new_ho_comparison/local_ho_data/

python/energy_flow_proper/09_dynamic_two_bath_one_qubit/energy_shovel_motor.org

@@ -0,0 +1,209 @@
+#+PROPERTY: header-args :session 09_energy_shovel :kernel python :pandoc no :async yes :tangle no
+Herein I'll try to apply the energy shovel from the previous project in a cycle.
+
+* Boilerplate
+#+begin_src jupyter-python :results none
+  import figsaver as fs
+  from hiro_models.one_qubit_model import QubitModelMutliBath, StocProcTolerances
+  import hiro_models.model_auxiliary as aux
+  import numpy as np
+  import qutip as qt
+#+end_src
+
+Init ray and silence stocproc.
+#+begin_src jupyter-python
+  import ray
+  ray.shutdown()
+  ray.init(address='auto')
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+: 2022-07-13 14:56:26,440       INFO worker.py:956 -- Connecting to existing Ray cluster at address: 141.30.17.16:6379
+: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '141.30.17.225', 'raylet_ip_address': '141.30.17.225', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-07-07_19-11-40_930040_398742/sockets/plasma_store.6', 'raylet_socket_name': '/tmp/ray/session_2022-07-07_19-11-40_930040_398742/sockets/raylet.2', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-07-07_19-11-40_930040_398742', 'metrics_export_port': 53497, 'gcs_address': '141.30.17.16:6379', 'address': '141.30.17.16:6379', 'node_id': 'a1b9acc4f221cb88f7c445762b2f1afb56a05388fd4a986377513e20'})
+:END:
+
+#+begin_src jupyter-python :results none
+  from hops.util.logging_setup import logging_setup
+  import logging
+  logging_setup(logging.INFO)
+#+end_src
+
+* Cycle
+#+begin_src jupyter-python :results none
+  from hops.util.dynamic_matrix import SmoothStep, Periodic, ConstantMatrix, ScaleTime, Shift, Piecewise
+#+end_src
+
+Now we build a basic otto cycle.
+Let's lay down some basic functionality.
+
+#+begin_src jupyter-python
+  H_op = (1/2 * (qt.sigmaz() + qt.identity(2))).full()
+  L_op = (1/2 * qt.sigmax()).full()
+  L_op_2 = (1/2 * qt.sigmay()).full()
+  print(H_op, L_op, sep="\n")
+#+end_src
+
+#+RESULTS:
+: [[1.+0.j 0.+0.j]
+:  [0.+0.j 0.+0.j]]
+: [[0. +0.j 0.5+0.j]
+:  [0.5+0.j 0. +0.j]]
+
+#+begin_src jupyter-python :results none
+  comp_ratio = .1
+  t_exp_end = .1
+  t_hot_modulation = .1
+  t_hot_end = .5
+  t_comp_end = .6
+  t_cold_modulation = .1
+  t_cold_end = 1
+  t = np.linspace(0, 10, 1000)
+#+end_src
+
+** Expansion
+#+begin_src jupyter-python
+  %matplotlib inline
+  H_exp = SmoothStep(comp_ratio * H_op, 0, t_exp_end, 3) + H_op * (1 - comp_ratio)
+  tt = np.linspace(0, 1, 1000)
+  plt.plot(tt, H_exp.derivative().operator_norm(tt))
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f92066c9e50> |
+[[file:./.ob-jupyter/586c97bae17d5d7684ecc3f9937e8183a3f3d3e1.svg]]
+:END:
+
+** Hot Thermalization
+#+begin_src jupyter-python
+  L_hot = SmoothStep(L_op, t_exp_end, t_exp_end + t_hot_modulation, 3) - SmoothStep(
+      L_op, t_hot_end - t_hot_modulation, t_hot_end, 3
+  )
+  plt.plot(tt, L_hot.operator_norm(tt))
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f9206382070> |
+[[file:./.ob-jupyter/7bed2d851ddc1cfc55317d2dff986662f609ca88.svg]]
+:END:
+
+** Compression
+#+begin_src jupyter-python
+  H_comp = ConstantMatrix(H_op) - SmoothStep(comp_ratio * H_op, t_hot_end, t_comp_end, 3)
+  plt.plot(tt, H_comp.operator_norm(tt))
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f9205fe38e0> |
+[[file:./.ob-jupyter/e6b4df328725cf5e6fca48e4d52a890dafefb35e.svg]]
+:END:
+
+** Cold Thermalization
+#+begin_src jupyter-python
+  L_cold = SmoothStep(L_op_2, t_comp_end, t_comp_end + t_cold_modulation, 3) - SmoothStep(
+      L_op_2, t_cold_end - t_cold_modulation, t_cold_end, 3
+  )
+  plt.plot(tt, L_cold.operator_norm(tt))
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f9205f7c3d0> |
+[[file:./.ob-jupyter/a4b3dcee48e76c03405b365f1513fae0fb337d47.svg]]
+:END:
+
+** Full Cycle
+Now we turn the system around after each fill-cycle.
+#+begin_src jupyter-python
+  H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_exp_end, 1]), 1)
+  L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)]
+  H_cyc, *L = [ScaleTime(op, 2) for op in [H_cyc, *L]]
+
+  plt.plot(t, H_cyc.operator_norm(t))
+  plt.plot(t, L[0].operator_norm(t))
+  plt.plot(t, L[1].operator_norm(t))
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7f9205edbca0> |
+[[file:./.ob-jupyter/bafb9c982efbee26edc0e3e20cd25059793ee4f3.svg]]
+:END:
+
+* Model
+#+begin_src jupyter-python :results none
+  model = QubitModelMutliBath(
+      δ=[1, 1],
+      ω_c=[1] * 2,
+      ω_s=[0, 10],
+      t=np.arange(0, 60, .01),
+      ψ_0=(qt.basis([2], [0])  + qt.basis([2], [1])*0),
+      description=f"A simple otto cycle.",
+      k_max=5,
+      bcf_terms=[5]*2,
+      truncation_scheme="simplex",
+      driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
+      thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2,
+      T = [0, 2],
+      L = L,
+      H = H_cyc * 10,
+      therm_methods=["tanhsinh", "fft"],
+  )
+#+end_src
+
+
+#+begin_src jupyter-python
+  aux.integrate(model, 100)
+#+end_src
+
+#+RESULTS:
+: fc74a430-fac5-467e-ab88-1bb8c9da7e74
+
+
+#+begin_src jupyter-python
+  %matplotlib inline
+  fig, ax = plt.subplots()
+  with aux.get_data(model) as data:
+      fs.plot_with_σ(model.t, model.total_energy_from_power(data), ax=ax, label=r"$\langle H(t)\rangle$")
+      fs.plot_with_σ(model.t, model.bath_energy(data).for_bath(0), ax=ax, label=r"$\langle H_{B}^{(1)}\rangle$")
+      fs.plot_with_σ(model.t, model.bath_energy(data).for_bath(1), ax=ax, label=r"$\langle H_{B}^{(2)}\rangle$")
+      fs.plot_with_σ(model.t, model.system_energy(data), ax=ax, label=r"$\langle H_{S}\rangle$")
+      # fs.plot_with_σ(model.t, model.total_power(data), ax=ax, label=r"Power")
+  # ax.plot(model.t, model.H(model.t)[:,0,0].real, label="$H_{\mathrm{sys}}(t)$")
+  # ax.plot(model.t, model.L[0](model.t)[:,0,1].real, label="$L_C(t)$")
+  # ax.plot(model.t, model.L[1](model.t)[:,0,1].real, label="$L_H(t)$")
+
+  ax.legend()
+#+end_src
+
+#+RESULTS:
+: 5179bed6-d693-43de-8969-bf0dedacc8f3
+
+#+begin_src jupyter-python :results none
+  with aux.get_data(model) as data:
+      ρ = data.rho_t_accum.mean[:]
+      σ_ρ = data.rho_t_accum.ensemble_std[:]
+      xs = np.einsum("tij,ji->t", ρ, qt.sigmax().full())
+      ys = np.einsum("tij,ji->t", ρ, qt.sigmay().full())
+      zs = np.einsum("tij,ji->t", ρ, qt.sigmaz().full())
+#+end_src
+
+
+#+begin_src jupyter-python
+  %matplotlib tk
+  b = qt.Bloch()
+  b.add_points([xs, ys, zs])
+  b.view = [20, 20]
+  b.point_size = [0.1]
+  b.sphere_alpha = 0.1
+  b.size = [20, 20]
+  b.render()
+  b.fig
+#+end_src
+
+
+#+RESULTS:
+[[file:./.ob-jupyter/ceddb27aeff1174b98ec49f0d1e712f58591b14c.svg]]