From 9bc1a1735a7176d1c2be63f07097acfa917fb857 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Thu, 1 Sep 2022 14:24:26 +0300
Subject: [PATCH] core: fix J
---
src/core.cpp | 15 +++++++++++++--
usadelndsoc/solver.py | 2 +-
2 files changed, 14 insertions(+), 3 deletions(-)
diff --git a/src/core.cpp b/src/core.cpp
index b2a7640..dc79a75 100644
--- a/src/core.cpp
+++ b/src/core.cpp
@@ -367,11 +367,22 @@ public:
Array<ADComplex, Shape<Dynamic,Dynamic,4,4,4> >
U({nx_, ny_, 4, 4, 4});
+ const size_t invdir[4] = {2, 3, 0, 1};
+ const int neigh[4][2] = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};
+
for (size_t i = 0; i < nx_; ++i) {
for (size_t j = 0; j < ny_; ++j) {
for (size_t k = 0; k < 4; ++k) {
- U.part(i,j,k).matrix() =
- U_.part(i,j,k).matrix() + J.part(i,j,k).matrix() * U_.part(i,j,k).matrix();
+ U.part(i,j,k).matrix() = U_.part(i,j,k).matrix();
+
+ const int i2 = i + neigh[k][0];
+ const int j2 = j + neigh[k][1];
+ if (i2 >= (int)nx_ || j2 >= (int)ny_)
+ continue;
+
+ U.part(i,j,k).matrix() = U_.part(i,j,k).matrix()
+ + J.part(i,j,k).matrix() * U_.part(i,j,k).matrix()
+ - U_.part(i,j,k).matrix() * J.part(i2,j2,invdir[k]).matrix();
}
}
}
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index 7d27edd..3d5cc75 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -596,7 +596,7 @@ class Result:
@property
def J(self):
- return self._core.grad_J(self.Phi).transpose(0, 1, 2, 4, 3) * (-1 / 8)
+ return self._core.grad_J(self.Phi).transpose(0, 1, 2, 4, 3) * (-1 / 16)
@property
def S(self):
--
GitLab