solver: improve docs for interpolate_J

usadelndsoc/solver.py

@@ -1530,12 +1530,6 @@ def interpolate_J(x0, y0, J, x, y, method="bilinear"):
     The current values are defined on links between cells, and not at
     the cell centers, so special interpolation is necessary.
 
-    The interpolant inside each cell is linear :math:`J(x,y) = (J_0x + x J_1,
-    J_0y + y J_2)`. The component perpendicular to the cell faces is
-    piecewise linear, but the component parallel to cell faces is
-    piecewise constant.  The interpolant is divergenceless inside
-    cells, if the matrix current values do not indicate divergence.
-
     Parameters
     ----------
     x0 : array, shape (nx,)
@@ -1551,12 +1545,34 @@ def interpolate_J(x0, y0, J, x, y, method="bilinear"):
     method : {'linear', 'bilinear'}
         Interpolation method.
 
+        - ``linear``: :math:`J_x` is a continuous linear function
+          along :math:`x` but piecewise constant along `y', and vice versa
+          for :math:`J_y`.
+
+        - ``bilinear``: :math:`J_x` and :math:`J_y` are continuous
+          piecewise linear functions both along :math:`x` and :math:`y`.
+
     Returns
     -------
     J_interp : array, shape broadcast_shape(x, y) + (2, ...)
         Array containing x and y components of the current,
         interpolated at the coordinates (x, y).
 
+    Notes
+    -----
+    For all of the interpolants, the integral of the interpolated
+    current across each cell edge equals the corresponding value in
+    the original *J* array.
+
+    Consequently, the average of the divergence of the interpolant
+    over a cell equals the sum of currents flowing out of the cell in
+    the original *J* array. Hence, the interpolated currents satisfy
+    similar conservation laws as the original ones.
+
+    Covariant conservation is however not exactly preserved.  Current
+    across an edge between cells is taken as the average of the
+    current on both sides.
+
     """
 
     dx = x0[1] - x0[0]
@@ -1615,7 +1631,8 @@ def interpolate_J(x0, y0, J, x, y, method="bilinear"):
 
         # Solve corner values Iyc: (Iyc(i,j) + Iyc(i+1,j))/2 == I3(i,j)
         r = np.empty((nx, ny + 1) + I0.shape[2:], dtype=I0.dtype)
-        r[:, :ny] = I3
+        r[:, 0] = I3[:, 0]
+        r[:, 1:ny] = (I3[:, 1:] + I1[:, :-1]) / 2
         r[:, ny] = I1[:, -1]
 
         rp = r.reshape(nx, -1)
@@ -1624,7 +1641,8 @@ def interpolate_J(x0, y0, J, x, y, method="bilinear"):
 
         # Solve corner values Ixc: (Ixc(i,j) + Ixc(i,j+1))/2 == I2(i,j)
         r = np.empty((nx + 1, ny) + I0.shape[2:], dtype=I0.dtype)
-        r[:nx, :] = I2
+        r[0, :] = I2[0, :]
+        r[1:nx, :] = (I2[1:] + I0[:-1]) / 2
         r[nx, :] = I0[-1, :]
 
         rp = np.moveaxis(r, 1, 0).reshape(ny, -1)