From 2ace65385a3fc63bbaaeff8f9332ce6930a5ff91 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Wed, 21 Sep 2022 15:35:33 +0300
Subject: [PATCH] examples: make cpr_sns produce a plot
---
examples/cpr_sns.py | 35 ++++++++++++++++++-----------------
1 file changed, 18 insertions(+), 17 deletions(-)
diff --git a/examples/cpr_sns.py b/examples/cpr_sns.py
index 371a2ce..5c1a543 100644
--- a/examples/cpr_sns.py
+++ b/examples/cpr_sns.py
@@ -15,6 +15,8 @@ from usadelndsoc.matsubara import get_matsubara_sum
from usadelndsoc.solver import *
from usadelndsoc.util import vectorize_parallel
+import collections
+
mem = joblib.Memory("cache")
@@ -49,9 +51,13 @@ def get_solver(soc_alpha, eta, phi, L=10, h=0.5, D=1.0, n=5):
return sol
-# @usadelndsoc.with_log_level(logging.WARNING)
+Res = collections.namedtuple("Res", ["x", "y", "J", "Jx", "Jy"])
+
+
@vectorize_parallel(returns_object=True, noarray=True)
+@usadelndsoc.with_log_level(logging.WARNING)
def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, perturbative=False):
+
sol = get_solver(soc_alpha=alpha_soc, eta=eta, L=L, n=n, phi=phi, h=h)
E_typical = 10 + 10 * T
@@ -67,29 +73,24 @@ def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, perturbative=False):
Jx = (J[:, :, 0] + J[:, :, 2]).real / 2
Jy = (J[:, :, 1] + J[:, :, 3]).real / 2
- return sol.x, sol.y, J, Jx, Jy
+ return Res(sol.x, sol.y, J, Jx, Jy)
def main():
- T = 0.02
- alpha_soc = 0.05
- Da2 = alpha_soc**2
- h = np.linspace(0, 1.75, 29)
+ T = 0.1
+ alpha_soc = 0.1
+ h = 0.5
+ phi = np.linspace(-pi, pi, 37)
- res = j_anom(T, h, alpha_soc=alpha_soc, perturbative=False, mem=mem)
- j = np.asarray([x[0] for x in res])
- j_an = np.asarray([x[1] for x in res])
+ res = j(T, h, phi, alpha_soc=alpha_soc, perturbative=False, mem=mem)
- t3 = np.diag([1, 1, -1, -1])
- j = tr(j @ t3)[:, 2, 2, UP]
+ Jx_mean = np.asarray([x.Jx[1:-1].sum(axis=1).mean(axis=0) for x in res])
- plt.plot(h, -j.real, "-", label=r"numerics")
- plt.plot(h, -j_an.real, "--", label=r"analytic")
- plt.title(rf"$T = {T} \Delta$, $D\alpha^2 = {Da2:.4} \Delta$")
- plt.xlabel(r"$h/\Delta$")
- plt.ylabel(r"$j_{\mathrm{an}}$")
+ plt.plot(phi / pi, Jx_mean)
+ plt.xlabel(r"$\varphi / \pi$")
+ plt.ylabel(r"$I$")
plt.legend()
- plt.savefig("bulk_j_an.pdf")
+ plt.savefig("cpr_sns.pdf")
if __name__ == "__main__":
--
GitLab