solver: fixup initial comparison residual with terminals

usadelndsoc/solver.py

@@ -159,6 +159,9 @@ class Solver:
         A_shape = (A_size, A_size)
         A_dtype = np.dtype(float)
 
+        Phi00 = Phi00.copy()
+        self._core.fix_terminals(Phi00)
+
         F0 = eval_rhs(Phi00)
         F0_norm = _norm(F0)
 
@@ -622,9 +625,15 @@ class Result:
             G = self.G
             return s * np.linalg.solve(I + gt @ g, I - gt @ g)
 
+    def _get_J(self, Phi):
+        return self._core.grad_A(Phi).transpose(0, 1, 2, 4, 3) * (1 / 16)
+
     @property
     def J(self):
-        return self._core.grad_A(self.Phi).transpose(0, 1, 2, 4, 3) * (1 / 16)
+        if np.asarray(self.omega).ndim == 0:
+            return self._get_J(self.Phi)
+        else:
+            return np.array([self._get_J(P) for P in self.Phi])
 
     @property
     def J_c(self):