From 4b72731bd65dd3bc5342f6e4a24684376558ad33 Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Wed, 26 Oct 2022 15:37:22 +0200
Subject: [PATCH] bug fixes in gen_data; doc string in final_plots

gen_data: bug fixes in automatic selection of nmax (preset was ignored)
final_plots: added doc string for export_omega5_interp (npz export)
---
 experiments.py |  2 ++
 final_plots.py | 36 ++++++++++++++++++++++++++++++++++++
 gen_data.py    | 22 ++++++++++++----------
 plot.py        | 34 ++++++++++++++++++++++++++++++++++
 4 files changed, 84 insertions(+), 10 deletions(-)

diff --git a/experiments.py b/experiments.py
index 922d00e..3bcc3b8 100644
--- a/experiments.py
+++ b/experiments.py
@@ -242,6 +242,8 @@ def kogan04_calibrate(dm, **kwargs):
             yscale = np.pi,
             calibrate_min = 1,
             calibrate_max = 1.6,
+            include_Ga = True,
+            integral_method = -15,
             **kwargs,
             )
 
diff --git a/final_plots.py b/final_plots.py
index d3af580..df0c95a 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -245,6 +245,42 @@ def export_omega5_interp(
         ac_res = 201,
         korder = 2,
         ):
+    """
+    Export interpolated data for given parameters.
+    Data are stored with units e=hbar=kB=Tkrg=1.
+    All results are exported as 2d arrays of the same shape.
+    Results have the naming convention
+
+        {observable}_{order}{padding}_{method}_{suffix}
+
+    observable:
+        gdc   DC differential conductance, may have suffix "_ispline"
+        gac   AC differential conductance, must have suffix "_ispline"
+        idc   average current, may have suffix "_spline"
+        iac   oscillating current, may have suffix "_spline"
+        phase AC phase
+        gamma Γ(E=0) (not an observable)
+
+    order:
+        o2    2nd order RG equations
+        o3    3rd order RG equations without Ga and with approximated integral
+        o3a   3rd order RG equations with approximated integral
+        o3p   3rd order RG equations with full integral
+
+    padding:
+        ""    no Floquet matrix extrapolation
+        "_p"  Floquet matrix extrapolation to mitigate truncation effects
+
+    method:
+        mu    no unitary transformation, oscillating chemical potentials
+        J     include oscillating voltage in coupling by unitary transformation
+
+    suffix
+        ""    interpolation is done using griddata
+        "_spline"   interpolation is done using spline
+        "_ispline"  data is derived from spline interpolation for the
+                    current (dc or ac)
+    """
     dm = DataManager()
     data = dm.list(
             omega = omega,
diff --git a/gen_data.py b/gen_data.py
index 86eb40b..3c2daef 100644
--- a/gen_data.py
+++ b/gen_data.py
@@ -18,13 +18,14 @@ from data_management import DataManager
 from kondo import Kondo
 
 
-def get_ref_nmax_omega5(dm, omega, vdc, vac, preset="normal"):
-    assert omega == 16.5372
+def get_ref_nmax(dm, omega, vdc, vac, preset="precise"):
     parameters = dict(
-            precise = dict(omega=16.5372, solver_tol_rel=1e-8, solver_tol_abs=1e-10, d=1e9, method="mu", padding=0, good_flags=0x1000, bad_flags=0x0ffc),
-            normal = dict(omega=16.5372, solver_tol_rel=1e-9, solver_tol_abs=1e-11, d=1e5, method="mu", padding=0, good_flags=0x1000, bad_flags=0x0ffc),
+            old = dict(omega=omega, solver_tol_rel=1e-8, solver_tol_abs=1e-10, d=1e9, method="mu", padding=0, good_flags=0x0000, bad_flags=0x1ffc, xL=0.5, voltage_branches=4),
+            precise = dict(omega=omega, solver_tol_rel=1e-8, solver_tol_abs=1e-10, d=1e9, method="mu", padding=0, good_flags=0x1000, bad_flags=0x0ffc, xL=0.5, voltage_branches=4),
+            normal = dict(omega=omega, solver_tol_rel=1e-9, solver_tol_abs=1e-11, d=1e5, method="mu", padding=0, good_flags=0x1000, bad_flags=0x0ffc, xL=0.5, voltage_branches=4),
             )
     data = dm.list(**parameters[preset])
+    assert data.size > 0
     vdc_arr, vac_arr = np.broadcast_arrays(vdc, vac)
     nmax = -np.ones(vdc_arr.shape, dtype=np.int64).reshape((-1,))
     for i, (vdc, vac) in enumerate(zip(vdc_arr.flat, vac_arr.flat)):
@@ -244,8 +245,8 @@ python -m frtrg_kondo.gen_data \\
     method_group.add_argument("--truncation_order", metavar="int", type=int,
             choices=(2,3), default=3,
             help = "Truncation order of RG equations.")
-    method_group.add_argument("--preset", metavar="str", type=str, nargs='+',
-            choices=("normal", "precise"), default="precise",
+    method_group.add_argument("--preset", metavar="str", type=str,
+            choices=("normal", "precise", "old"), default="precise",
             help = "Preset for choosing nmax if set explicitly")
 
     # Numerical parameters concerning Floquet matrices
@@ -338,6 +339,7 @@ python -m frtrg_kondo.gen_data \\
             atol = options.pop("atol"),
             method = options.pop("solver_method"),
             )
+    preset = options.pop("preset", "precise")
 
     # Detect number of CPUs (if necessary)
     if threads == 0:
@@ -353,7 +355,7 @@ python -m frtrg_kondo.gen_data \\
                 assert np.abs(kondo_options["vdc"] - vdc) < 1e-9
                 kondo_options["vdc"] = vdc
             if kondo_options["nmax"] < 0:
-                kondo_options["nmax"] = get_ref_nmax_omega5(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"])
+                kondo_options["nmax"] = get_ref_nmax(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"], preset)
             if kondo_options["nmax"] < 0:
                 settings.logger.error(f"Invalid value nmax={nmax}: must be ≥0")
                 continue
@@ -373,7 +375,7 @@ python -m frtrg_kondo.gen_data \\
         lock = mp.Lock()
         queue = mp.Queue()
         # create processes
-        processes = [mp.Process(target=save_data_mp, args=(queue, lock, solver_options, filename, include, find_pole)) for i in range(threads)]
+        processes = [mp.Process(target=save_data_mp, args=(queue, lock, solver_options, filename, include, find_pole, False, preset)) for i in range(threads)]
         # start processes
         for p in processes:
             p.start()
@@ -385,7 +387,7 @@ python -m frtrg_kondo.gen_data \\
             queue.put(None)
 
 
-def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole=False, overwrite=False):
+def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole=False, overwrite=False, preset="precise"):
     """
     Generate data points in own process and save them to HDF5 file.
     Each process owns one DataManager instance.
@@ -401,7 +403,7 @@ def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole
             assert np.abs(kondo_options["vdc"] - vdc) < 1e-9
             kondo_options["vdc"] = vdc
         if kondo_options["nmax"] < 0:
-            kondo_options["nmax"] = get_ref_nmax_omega5(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"])
+            kondo_options["nmax"] = get_ref_nmax(dm, kondo_options["omega"], kondo_options["vdc"], kondo_options["vac"], preset)
         if kondo_options["nmax"] < 0:
             settings.logger.error(f"Invalid value nmax={kondo_options['nmax']}: must be ≥0")
             continue
diff --git a/plot.py b/plot.py
index aec9255..8611536 100644
--- a/plot.py
+++ b/plot.py
@@ -1721,6 +1721,40 @@ def RGeq_comparison_contours(
     plt.show()
 
 
+def RGeq_comparison_contours2(
+        dm = None,
+        filename = "figdata/omega5_interp.npz",
+        observable = "gdc",
+        prefactor = np.pi,
+        **trashparams):
+    omega = 16.5372
+    kwargs = dict(linewidths=0.8)
+    colors = dict(o2="#ff0000", o3="#6000ff", o3a="#00f000", o3p="#000000", ref="#808080")
+    zorder = dict(o2=1, o3=2, o3a=3, o3p=4)
+    approximations = list(zorder.keys())
+    data = np.load(filename)
+    vdc = data["vdc"] / omega
+    vac = data["vac"] / omega
+    levels = (1/np.linspace(1, 12, 12)[::-1]) / prefactor
+    fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True, sharey=True)
+    for approx in approximations:
+        ax.contour(vdc, vac, data[f"{observable}_mu_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+        try:
+            ax.contour(vdc, vac, data[f"{observable}_J_{approx}"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+        except KeyError:
+            settings.logger.warning(f"No data available for {observable}_J_{approx}")
+        try:
+            ax.contour(vdc, vac, data[f"{observable}_J_{approx}_p"], levels, zorder=zorder[approx], colors=colors[approx], **kwargs)
+        except KeyError:
+            settings.logger.warning(f"No data available for {observable}_J_{approx}_p")
+
+    #ax.contour(-vdc, vac, data[f"{observable}_mu_o3p"], levels, zorder=-1, colors=colors["ref"], **kwargs)
+    #ax.vlines((0,), 0, vac[-1], color="#808080", linewidth=0.5)
+    ax.set_xlabel(r"$V_\mathrm{avg}~(\Omega)$")
+    ax.set_ylabel(r"$V_\mathrm{osc}~(\Omega)$")
+    plt.show()
+
+
 def plot_Ga_max(
         dm,
         omega = 16.5372,
-- 
GitLab