diff --git a/.gitignore b/.gitignore index 5c6461f..1cc13fe 100644 --- a/.gitignore +++ b/.gitignore @@ -506,3 +506,4 @@ results/**.npy /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/10_antizeno_engine/taurus/ +/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/taurus/ diff --git a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.envrc b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.envrc index c7a853c..d3b3289 100644 --- a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.envrc +++ b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.envrc @@ -1,2 +1,2 @@ -use_flake +use_flake --impure eval "$shellHook" diff --git a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/first_experiments.org b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/first_experiments.org index 065f2d0..04a639e 100644 --- a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/first_experiments.org +++ b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/first_experiments.org @@ -12,6 +12,7 @@ with the markovian result. import qutip as qt import utilities as ut import stocproc + import matplotlib.pyplot as plt #+end_src Init ray and silence stocproc. @@ -23,7 +24,7 @@ Init ray and silence stocproc. #+end_src #+RESULTS: -: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '192.168.100.170', 'raylet_ip_address': '192.168.100.170', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-15_13-45-28_971048_59572/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-15_13-45-28_971048_59572/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-15_13-45-28_971048_59572', 'metrics_export_port': 53281, 'gcs_address': '192.168.100.170:61047', 'address': '192.168.100.170:61047', 'node_id': '19344eac29e15d6404c14920662bd8259981954f0a9fde73ef5f8b4b'}) +: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '192.168.100.170', 'raylet_ip_address': '192.168.100.170', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-15_16-33-19_035369_76479/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-15_16-33-19_035369_76479/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-15_16-33-19_035369_76479', 'metrics_export_port': 56879, 'gcs_address': '192.168.100.170:64999', 'address': '192.168.100.170:64999', 'node_id': 'cd3c35a56063d205e53b29934473b89b3baea526868d0ff8e806b08e'}) #+begin_src jupyter-python :results none from hops.util.logging_setup import logging_setup @@ -55,7 +56,7 @@ Let's lay down some basic functionality. #+begin_src jupyter-python :results none t_compression = .1 t_expansion = .1 - comp_ratio = .25 + comp_ratio = .5 hot_switch_ratio = .8 cold_switch_ratio = .8 @@ -75,8 +76,8 @@ Let's lay down some basic functionality. #+RESULTS: :RESULTS: -| | -[[file:./.ob-jupyter/3cc22c46295a7b959768d60bfea0e4c0142e2b43.svg]] +| | +[[file:./.ob-jupyter/8326a30a0b67dac1fcbfa726d3677ad081f50681.svg]] :END: ** Hot Thermalization @@ -89,8 +90,8 @@ Let's lay down some basic functionality. #+RESULTS: :RESULTS: -| | -[[file:./.ob-jupyter/98aec1955458432d8e29d21e1675b698b47b50fb.svg]] +| | +[[file:./.ob-jupyter/beef0051823ad8d69a03749d66188767c06f00d6.svg]] :END: ** Compression @@ -101,8 +102,8 @@ Let's lay down some basic functionality. #+RESULTS: :RESULTS: -| | -[[file:./.ob-jupyter/8ebe76752707b723232c25684baa8b2b3a0af60e.svg]] +| | +[[file:./.ob-jupyter/2c1c45c42776db17721492efacd6f4856f8fe3b3.svg]] :END: ** Cold Thermalization @@ -115,8 +116,8 @@ Let's lay down some basic functionality. #+RESULTS: :RESULTS: -| | -[[file:./.ob-jupyter/ffdfdc73f371b6485c38c7714e3d6f2075350606.svg]] +| | +[[file:./.ob-jupyter/a425435a0d34373dd92cf2a415ceec18a1fd69af.svg]] :END: ** Full Cycle @@ -125,29 +126,73 @@ Now we turn the system around after each fill-cycle. H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_expansion, 1]), 1) L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)] period = 1 - scale = .1 + scale = .05 H_cyc, *L = [ScaleTime(op, scale) for op in [H_cyc, *L]] ω_mod = 2*np.pi / (period) * scale - plt.plot(t, H_cyc.operator_norm(t)) - plt.plot(t, L[0].operator_norm(t)) - plt.plot(t, L[1].operator_norm(t)) ω_mod t_cycle = period/scale - periods = 20 + periods = 10 t_max = t_cycle * periods - dt = .01 + dt = .1 + + t = np.linspace(0, t_cycle, 1000) + + plt.plot(t, H_cyc.operator_norm(t), label=r"$||H_\mathrm{S}||$") + plt.plot(t, L[0].operator_norm(t), label=r"$||L_\mathrm{c}||$") + plt.plot(t, L[1].operator_norm(t), label=r"$||L_\mathrm{h}||$") + plt.xlabel(r"$\tau$") + plt.legend() + fs.export_fig("modulation") #+end_src #+RESULTS: -[[file:./.ob-jupyter/c6fcbae6c5d076e45f48b0d4e003265146c48409.svg]] +[[file:./.ob-jupyter/2b79808a3127f0a93ce6cbde00ae4888c83b2696.svg]] + +** Shifts +#+begin_src jupyter-python + from scipy.optimize import minimize_scalar + + ref_model = QubitModelMutliBath( + ω_c=[1] * 2, + T=[1, 20], + ) + + shifts = [] + scales = [] + + for i, t_exp in enumerate([0, t_expansion / scale]): + ω_exp = H_cyc.operator_norm(i) + + def objective(ω_s): + ref_model.ω_s[i] = ω_s + + return -ref_model.full_thermal_spectral_density(i)(ω_exp) + + res = minimize_scalar( + objective, 1, method="bounded", bounds=(0.01, ω_exp), options=dict(maxiter=100) + ) + shifts.append(res.x) + scales.append(-res.fun) + ω = np.linspace(0.1, 10, 1000) + plt.plot(ω, ref_model.full_thermal_spectral_density(i)(ω) / -res.fun) + + print(shifts, scales) +#+end_src + +#+RESULTS: +:RESULTS: +: [0.010005778765750251, 0.5000003902186526] [0.3656483223378975, 3.198956858189882] +[[file:./.ob-jupyter/97d3960ba8288d117b1bf5774e57ac8bb1500b83.svg]] +:END: + * Model #+begin_src jupyter-python :results none model = QubitModelMutliBath( - δ=[.05*4, .01*4], - ω_c=[1] * 2, - ω_s=[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1], + δ=list(1/np.array(scales) * .1/2),#[.05*4, .01*4], + ω_c=ref_model.ω_c, + ω_s=shifts,#[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1], t=ut.linspace_with_strobe(0, t_max, int(t_max // dt), ω_mod), ψ_0=(qt.basis([2], [0])*0 + qt.basis([2], [1])*1), description=f"A simple otto cycle.", @@ -156,95 +201,60 @@ Now we turn the system around after each fill-cycle. truncation_scheme="simplex", driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2, thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2, - T = [1, 20], + T = ref_model.T, #= [1, 20], L = L, H = H_cyc, - therm_methods=["tanhsinh", "fft"], + therm_methods=["fft", "fft"], ) #+end_src #+begin_src jupyter-python - aux.integrate(model, 1000) + aux.integrate(model, 10000) #+end_src #+RESULTS: -: [INFO hops.core.integration 59572] Choosing the nonlinear integrator. -: [INFO hops.core.integration 59572] Using 8 integrators. -: [INFO hops.core.integration 59572] Some 0 trajectories have to be integrated. -: [INFO hops.core.integration 59572] Using 1820 hierarchy states. -: 0it [00:00, ?it/s] +: [INFO hops.core.integration 76479] Choosing the nonlinear integrator. +: [INFO hops.core.integration 76479] Using 8 integrators. +: [INFO hops.core.integration 76479] Some 30 trajectories have to be integrated. +: [INFO hops.core.integration 76479] Using 1820 hierarchy states. +: 100% 30/30 [00:22<00:00, 1.32it/s] #+begin_src jupyter-python - f, a = pu.plot_energy_overview(model) - a.plot(model.t, model.H.operator_norm(model.t)) - a.plot(model.t, model.L[0].operator_norm(model.t)) - a.plot(model.t, model.L[1].operator_norm(model.t)) + f, a = pu.plot_energy_overview(model, strobe_frequency=ω_mod) + # a.plot(model.t, model.H.operator_norm(model.t)) + # a.plot(model.t, model.L[0].operator_norm(model.t)) + # a.plot(model.t, model.L[1].operator_norm(model.t)) + ax.set_xlabel(r"$\tau$") a.legend() + #fs.export_fig("energy_strobe", y_scaling=.6) #+end_src #+RESULTS: -:RESULTS: -#+begin_example - Loading: 100% 8/8 [00:00<00:00, 192.76it/s] - [INFO root 59572] Writing cache to: results/5e8c8c2b39bd3b8ad75fb0da753ca3989cec2ca758e5b60d542cfb0c39e6009d_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_op_exp_task_33_None_7b4b9e9cfa38e0c95465b4087684ba9103d527ed24c440718c98f01d3693d041.npy - Loading: 100% 16/16 [00:00<00:00, 50.04it/s] - [INFO root 59572] Writing cache to: results/flow_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_flow_worker_33_None_c15d067b912279587be3d3939223173fe5368154834a5a60648fc0a46950f8cf.npy - Loading: 100% 16/16 [00:00<00:00, 86.42it/s] - [INFO root 59572] Writing cache to: results/interaction_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_interaction_energy_task_33_None_58455bf5ce6e35b96d04516e2a659a907e0c44b588167e00faa3fdf3e226bb11.npy - Loading: 100% 16/16 [00:00<00:00, 79.44it/s] - [INFO root 59572] Writing cache to: results/interaction_power_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_interaction_energy_task_33_None_58455bf5ce6e35b96d04516e2a659a907e0c44b588167e00faa3fdf3e226bb11.npy - Loading: 100% 8/8 [00:00<00:00, 245.72it/s] - [INFO root 59572] Writing cache to: results/143e514d782a0bf3369be86475c494b0df4baa9952cf06bbf576d262f2dd27a4_12db9adeba3a9d05eee93dba02bb24bbd2bd553fa5e391b4607f3101447ce076_op_exp_task_33_None_7b4b9e9cfa38e0c95465b4087684ba9103d527ed24c440718c98f01d3693d041.npy -#+end_example -: -[[file:./.ob-jupyter/229ace105ac0a848d7b968ebdfd79d9fd38dbe9e.svg]] -:END: +[[file:./.ob-jupyter/5d6c6959a02462528cdb3d541f8dd5ecf303f362.svg]] #+begin_src jupyter-python fig, ax = plt.subplots() with aux.get_data(model) as data: - power = model.interaction_power(data).for_bath(0) + power = model.interaction_power(data).sum_baths() power_sys = model.system_power(data).sum_baths() energy = model.total_energy_from_power(data) - pu.plot_with_σ(model.t, power, ax=ax) - pu.plot_with_σ(model.t, power_sys, ax=ax, label="sys") - pu.plot_with_σ(model.t, energy, ax=ax) + pu.plot_with_σ(model.t, power, ax=ax, label="Interaction") + pu.plot_with_σ(model.t, power_sys, ax=ax, label="System") ax.legend() + ax.set_title("Power") + ax.set_xlabel(r"$\tau$") + fs.export_fig("power", y_scaling=.4) #+end_src #+RESULTS: -:RESULTS: -# [goto error] -#+begin_example - --------------------------------------------------------------------------- - RuntimeError Traceback (most recent call last) - Input In [97], in () -  1 fig, ax = plt.subplots() - ----> 2 with aux.get_data(model) as data: -  3 power = model.interaction_power(data).for_bath(0) -  4 power_sys = model.system_power(data).sum_baths() - - File ~/src/two_qubit_model/hiro_models/model_auxiliary.py:146, in get_data(model, data_path, read_only, **kwargs) -  135 return HIData( -  136 path, -  137 hi_key=model.hops_config, -  (...) -  142 **kwargs, -  143 ) -  145 else: - --> 146 raise RuntimeError(f"No data found for model with hash '{hexhash}'.") - - RuntimeError: No data found for model with hash '2676c736eb2bf1d2a1f37e21cda1823a3dc6b3f2cedb82509cfa6231a5232b2f'. -#+end_example -[[file:./.ob-jupyter/5285f6b51867d45830747c7d1501193929fe8fd1.svg]] -:END: +[[file:./.ob-jupyter/4938909f0a3699e7c9912ab201d2f04199b355e7.svg]] #+begin_src jupyter-python - steady_index = 6 + steady_index = 3 cycle_times, cycle_indices = ut.strobe_times(model.t, ω_mod) @@ -257,7 +267,7 @@ Now we turn the system around after each fill-cycle. ) e_diff = total_energies.slice(0) - total_energies.slice(-1) - mean_power = e_diff * (1 / (model.t[-1] - model.t[steady_index])) + mean_power = e_diff * (1 / (model.t[-1] - cycle_times[steady_index])) efficiency = e_diff * ( 1 / (hot_bath_energies.slice(0) - hot_bath_energies.slice(-1)).value ) @@ -265,13 +275,13 @@ Now we turn the system around after each fill-cycle. #+end_src #+RESULTS: -| -0.011704395396932077 | 8.941753951852788e-05 | -0.6487259446355118 | +| 0.0019806021074448648 | 2.3085106663538818e-05 | 0.1975630167587971 | #+begin_src jupyter-python with aux.get_data(model) as data: ρ = data.rho_t_accum.mean[:] σ_ρ = data.rho_t_accum.ensemble_std[:] - plt.plot(model.t, ρ[:, 1, 1]) + plt.plot(model.t, ρ[:, 0, 0]) 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()) @@ -279,9 +289,9 @@ Now we turn the system around after each fill-cycle. #+RESULTS: :RESULTS: -: /nix/store/dbz51ji4p5pcppra5h7yqhwcrbqgq40v-python3-3.9.13-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part +: /nix/store/yf7vf1pic6xfwr035xywcczjm08psvpa-python3-3.9.13-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part : return np.asarray(x, float) -[[file:./.ob-jupyter/9b6dbfd5d863d78fc20ce574cf7e282bc8341b0c.svg]] +[[file:./.ob-jupyter/69c7da5a4fa9adc0fb49b0ff097ea6fbc28c58ca.svg]] :END: @@ -299,9 +309,9 @@ Now we turn the system around after each fill-cycle. #+RESULTS: :RESULTS: -: /nix/store/dbz51ji4p5pcppra5h7yqhwcrbqgq40v-python3-3.9.13-env/lib/python3.9/site-packages/qutip/bloch.py:639: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword +: /nix/store/yf7vf1pic6xfwr035xywcczjm08psvpa-python3-3.9.13-env/lib/python3.9/site-packages/qutip/bloch.py:639: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword : self.fig.canvas.draw() -: /nix/store/dbz51ji4p5pcppra5h7yqhwcrbqgq40v-python3-3.9.13-env/lib/python3.9/site-packages/IPython/core/pylabtools.py:151: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword +: /nix/store/yf7vf1pic6xfwr035xywcczjm08psvpa-python3-3.9.13-env/lib/python3.9/site-packages/IPython/core/pylabtools.py:151: UserWarning: There are no gridspecs with layoutgrids. Possibly did not call parent GridSpec with the "figure" keyword : fig.canvas.print_figure(bytes_io, **kw) -[[file:./.ob-jupyter/7806ea3a984dc21f51c7f09f8720ae8e361465a5.svg]] +[[file:./.ob-jupyter/633b45a276b59c015a575e42b611b52283311ebe.svg]] :END: diff --git a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/mount_taurus.sh b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/mount_taurus.sh new file mode 100755 index 0000000..baf19b2 --- /dev/null +++ b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/mount_taurus.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env bash +sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_11/project/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/.data .data +sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_11/project/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/results results + + + +sshfs -oIdentityFile=~/.ssh/id_ed25519_taurus s8896854@taurusexport.hrsk.tu-dresden.de:/lustre/ssd/ws/s8896854-m_11/project/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/ taurus diff --git a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor new file mode 100644 index 0000000..316f1a3 --- /dev/null +++ b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor @@ -0,0 +1,125 @@ +import figsaver as fs +import plot_utils as pu +from hiro_models.one_qubit_model import QubitModelMutliBath, StocProcTolerances +import hiro_models.model_auxiliary as aux +import numpy as np +import qutip as qt +import utilities as ut +import stocproc + +import ray +ray.shutdown() +#ray.init(address='auto') +ray.init() + +from hops.util.logging_setup import logging_setup +import logging +logging_setup(logging.INFO) + +from hops.util.dynamic_matrix import SmoothStep, Periodic, ConstantMatrix, ScaleTime, Shift, Piecewise + +H_op = 2*(1/2 * (qt.sigmaz() + qt.identity(2))).full() +L_op = (1/2 * qt.sigmax()).full() +L_op_2 = (1/2 * qt.sigmax()).full() +print(H_op, L_op, sep="\n") + +t_compression = .1 +t_expansion = .1 +comp_ratio = .5 +hot_switch_ratio = .8 +cold_switch_ratio = .8 + +t_hot = (1 - t_expansion - t_compression) / 2 +t_hot_switch = t_hot * hot_switch_ratio / 2 +t_cold = (1 - t_expansion - t_compression) / 2 +t_cold_switch = t_cold * cold_switch_ratio / 2 +t = np.linspace(0, 10, 1000) + +%matplotlib inline +H_exp = SmoothStep(comp_ratio * H_op, 0, t_expansion, 2) + H_op * (1 - comp_ratio) +tt = np.linspace(0, 1, 1000) +plt.plot(tt, H_exp.operator_norm(tt)) + +L_hot = SmoothStep(L_op, t_expansion, t_expansion + t_hot_switch, 2) - SmoothStep( + L_op, t_expansion + t_hot - t_hot_switch, t_expansion + t_hot, 2 + ) +plt.plot(tt, L_hot.operator_norm(tt)) + +H_comp = ConstantMatrix(H_op) - SmoothStep(comp_ratio * H_op, t_expansion + t_hot, t_expansion + t_hot + t_compression, 2) +plt.plot(tt, H_comp.operator_norm(tt)) + +L_cold = SmoothStep(L_op_2, t_expansion + t_hot + t_compression, t_expansion + t_hot + t_compression + t_cold_switch, 3) - SmoothStep( + L_op_2, t_expansion + t_hot + t_compression + t_cold - t_cold_switch, 1, 2 +) +plt.plot(tt, L_cold.operator_norm(tt)) + +H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_expansion, 1]), 1) +L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)] +period = 1 +scale = .1 +H_cyc, *L = [ScaleTime(op, scale) for op in [H_cyc, *L]] + +ω_mod = 2*np.pi / (period) * scale +plt.plot(t, H_cyc.operator_norm(t)) +plt.plot(t, L[0].operator_norm(t)) +plt.plot(t, L[1].operator_norm(t)) +ω_mod +t_cycle = period/scale +periods = 20 +t_max = t_cycle * periods +dt = .01 + +model = QubitModelMutliBath( + δ=[.5, .1], + ω_c=[1] * 2, + ω_s=[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1], + t=ut.linspace_with_strobe(0, t_max, int(t_max // dt), ω_mod), + ψ_0=(qt.basis([2], [0])*0 + qt.basis([2], [1])*1), + description=f"A simple otto cycle.", + k_max=4, + bcf_terms=[6]*2, + truncation_scheme="simplex", + driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2, + thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2, + T = [0, 20], + L = L, + H = H_cyc, + therm_methods=["tanhsinh", "fft"], +) + +aux.integrate(model, 200) + +%matplotlib inline +f, a = pu.plot_energy_overview(model) +a.plot(model.t, model.H.operator_norm(model.t)) +a.plot(model.t, model.L[0].operator_norm(model.t)) +a.plot(model.t, model.L[1].operator_norm(model.t)) +a.legend() + +fig, ax = plt.subplots() +with aux.get_data(model) as data: + power = model.interaction_power(data).sum_baths() + power_sys = model.system_power(data).sum_baths() + energy = model.total_energy_from_power(data).sum_baths() + pu.plot_with_σ(model.t, power, ax=ax) + pu.plot_with_σ(model.t, power_sys, ax=ax, label="sys") + pu.plot_with_σ(model.t, energy, ax=ax) + +ax.legend() + +with aux.get_data(model) as data: + ρ = data.rho_t_accum.mean[:] + σ_ρ = data.rho_t_accum.ensemble_std[:] + plt.plot(model.t, ρ[:, 1, 1]) + 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()) + +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 diff --git a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor.py b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor.py index a4e69f2..4e0bedc 100644 --- a/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor.py +++ b/python/energy_flow_proper/09_dynamic_two_bath_one_qubit/otto_motor.py @@ -6,6 +6,7 @@ import numpy as np import qutip as qt import utilities as ut import stocproc +import matplotlib.pyplot as plt import ray ray.shutdown() @@ -25,7 +26,7 @@ print(H_op, L_op, sep="\n") t_compression = .1 t_expansion = .1 -comp_ratio = .25 +comp_ratio = .5 hot_switch_ratio = .8 cold_switch_ratio = .8 @@ -55,23 +56,57 @@ plt.plot(tt, L_cold.operator_norm(tt)) H_cyc = Periodic(Piecewise([H_exp, H_comp], [0, t_expansion, 1]), 1) L = [Periodic(L_i, 1) for L_i in (L_cold, L_hot)] period = 1 -scale = .1 +scale = .05 H_cyc, *L = [ScaleTime(op, scale) for op in [H_cyc, *L]] ω_mod = 2*np.pi / (period) * scale -plt.plot(t, H_cyc.operator_norm(t)) -plt.plot(t, L[0].operator_norm(t)) -plt.plot(t, L[1].operator_norm(t)) ω_mod t_cycle = period/scale -periods = 20 +periods = 10 t_max = t_cycle * periods -dt = .01 +dt = .1 + +t = np.linspace(0, t_cycle, 1000) + +plt.plot(t, H_cyc.operator_norm(t), label=r"$||H_\mathrm{S}||$") +plt.plot(t, L[0].operator_norm(t), label=r"$||L_\mathrm{c}||$") +plt.plot(t, L[1].operator_norm(t), label=r"$||L_\mathrm{h}||$") +plt.xlabel(r"$\tau$") +plt.legend() +fs.export_fig("modulation") + +from scipy.optimize import minimize_scalar + +ref_model = QubitModelMutliBath( + ω_c=[1] * 2, + T=[1, 20], +) + +shifts = [] +scales = [] + +for i, t_exp in enumerate([0, t_expansion / scale]): + ω_exp = H_cyc.operator_norm(i) + + def objective(ω_s): + ref_model.ω_s[i] = ω_s + + return -ref_model.full_thermal_spectral_density(i)(ω_exp) + + res = minimize_scalar( + objective, 1, method="bounded", bounds=(0.01, ω_exp), options=dict(maxiter=100) + ) + shifts.append(res.x) + scales.append(-res.fun) + ω = np.linspace(0.1, 10, 1000) + plt.plot(ω, ref_model.full_thermal_spectral_density(i)(ω) / -res.fun) + +print(shifts, scales) model = QubitModelMutliBath( - δ=[.05*4, .01*4], - ω_c=[1] * 2, - ω_s=[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1], + δ=list(1/np.array(scales) * .1/2),#[.05*4, .01*4], + ω_c=ref_model.ω_c, + ω_s=shifts,#[H_cyc.operator_norm(0) - 1, H_cyc.operator_norm(t_expansion/scale) - 1], t=ut.linspace_with_strobe(0, t_max, int(t_max // dt), ω_mod), ψ_0=(qt.basis([2], [0])*0 + qt.basis([2], [1])*1), description=f"A simple otto cycle.", @@ -80,32 +115,36 @@ model = QubitModelMutliBath( truncation_scheme="simplex", driving_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2, thermal_process_tolerances=[StocProcTolerances(1e-3, 1e-3)] * 2, - T = [1, 20], + T = ref_model.T, #= [1, 20], L = L, H = H_cyc, - therm_methods=["tanhsinh", "fft"], + therm_methods=["fft", "fft"], ) -aux.integrate(model, 1000) +aux.integrate(model, 10000) -f, a = pu.plot_energy_overview(model) -a.plot(model.t, model.H.operator_norm(model.t)) -a.plot(model.t, model.L[0].operator_norm(model.t)) -a.plot(model.t, model.L[1].operator_norm(model.t)) +f, a = pu.plot_energy_overview(model, strobe_frequency=ω_mod) +# a.plot(model.t, model.H.operator_norm(model.t)) +# a.plot(model.t, model.L[0].operator_norm(model.t)) +# a.plot(model.t, model.L[1].operator_norm(model.t)) +ax.set_xlabel(r"$\tau$") a.legend() +#fs.export_fig("energy_strobe", y_scaling=.6) fig, ax = plt.subplots() with aux.get_data(model) as data: - power = model.interaction_power(data).for_bath(0) + power = model.interaction_power(data).sum_baths() power_sys = model.system_power(data).sum_baths() energy = model.total_energy_from_power(data) - pu.plot_with_σ(model.t, power, ax=ax) - pu.plot_with_σ(model.t, power_sys, ax=ax, label="sys") - pu.plot_with_σ(model.t, energy, ax=ax) + pu.plot_with_σ(model.t, power, ax=ax, label="Interaction") + pu.plot_with_σ(model.t, power_sys, ax=ax, label="System") ax.legend() +ax.set_title("Power") +ax.set_xlabel(r"$\tau$") +fs.export_fig("power", y_scaling=.4) -steady_index = 6 +steady_index = 3 cycle_times, cycle_indices = ut.strobe_times(model.t, ω_mod) @@ -118,7 +157,7 @@ with aux.get_data(model) as data: ) e_diff = total_energies.slice(0) - total_energies.slice(-1) -mean_power = e_diff * (1 / (model.t[-1] - model.t[steady_index])) +mean_power = e_diff * (1 / (model.t[-1] - cycle_times[steady_index])) efficiency = e_diff * ( 1 / (hot_bath_energies.slice(0) - hot_bath_energies.slice(-1)).value ) @@ -127,7 +166,7 @@ mean_power.value, mean_power.σ, efficiency.value with aux.get_data(model) as data: ρ = data.rho_t_accum.mean[:] σ_ρ = data.rho_t_accum.ensemble_std[:] - plt.plot(model.t, ρ[:, 1, 1]) + plt.plot(model.t, ρ[:, 0, 0]) 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())