diff --git a/TPFA_ResSim/__init__.py b/TPFA_ResSim/__init__.py
index f246ec7..eb50d2c 100644
--- a/TPFA_ResSim/__init__.py
+++ b/TPFA_ResSim/__init__.py
@@ -127,7 +127,7 @@ def _set_Q(self, k):
             for xy, q in zip(xys, rates):
                 # Use += in case of superimposed wells (e.g. by optimzt)
                 Q[self.xy2ind(*xy)] += sign * q
-        assert np.isclose(Q.sum(), 0), "Inj - Prod does not sum to 0"
+        assert np.isclose(Q.sum(), 0), "(Inj - Prod) does not sum to 0"
         self._Q = Q
 
     # Pres() -- listing 5
diff --git a/TPFA_ResSim/grid.py b/TPFA_ResSim/grid.py
index 76b05a0..36eab20 100644
--- a/TPFA_ResSim/grid.py
+++ b/TPFA_ResSim/grid.py
@@ -7,13 +7,13 @@
 The index ordering is "C-style" (numpy default).
 This choice means that `x` is the 1st coord., `y` is 2nd,
 and is hardcoded in the reservoir simulator model code
-in what takes place **between** `np.ravel` and `np.reshape`
-(both of which are configured to use row-major index ordering).
+(in what takes place **between** `np.ravel` and `np.reshape`,
+both of which are configured to use row-major index ordering.
+"F-style" (column-major) indexing implementation is perfectly possible,
+but would imply an undue amount hassle).
 Conveniently, it also means that `x` and `y` tend to occur in alphabetic order.
 Thus, in printing a matrix of a field, the `x` coordinate corresponds to the row index.
 By contrast, the plotting module depicts `x` from left to right, `y` from bottom to top.
-Implementing support for "F-style" (column-major) indexing is possible,
-but would imply an undue amount hassle.
 """
 
 from dataclasses import dataclass
@@ -68,7 +68,7 @@ def size(self):
 
     @property
     def domain(self):
-        """`(0, 0, Lx, Ly)`"""
+        """`((0, 0), (Lx, Ly))`"""
         return ((0, 0), (self.Lx, self.Ly))
 
     @property
diff --git a/TPFA_ResSim/plotting.py b/TPFA_ResSim/plotting.py
index 1603d16..9aa5238 100644
--- a/TPFA_ResSim/plotting.py
+++ b/TPFA_ResSim/plotting.py
@@ -56,16 +56,14 @@ def plt_field(self, ax, Z, style="default", wells=True,
 
         `kwargs` falls back to `styles[style]`, which falls back to `styles['defaults']`.
         """
-        # Set defaults
+        # Populate kwargs with fallback style
         kwargs = {**styles["default"], **styles[style], **kwargs}
         # Pop from kwargs. Remainder goes to countourf
         ax.set(**axprops(kwargs))
         cticks = kwargs.pop("cticks")
 
-        # Plotting with extent=(0, Lx, 0, Ly), rather than merely changing ticks
-        # has the advantage that set_aspect("equal") yields correct axes size,
-        # and that mouse hovering (with interactive backends) reports correct pos.
-        # Disadvantage: well_scatter must also account for coord_type.
+        # Why extent=(0, Lx, 0, Ly), rather than merely changing ticks?
+        # set_aspect("equal") and mouse hovering (reporting x,y).
         if "rel" in coord_type:
             Lx, Ly = 1, 1
         elif "abs" in coord_type:
@@ -75,19 +73,27 @@ def plt_field(self, ax, Z, style="default", wells=True,
         else:
             raise ValueError(f"Unsupported coord_type: {coord_type}")
 
-        # Need to transpose coz self assumes shape (Nx, Ny),
-        # and contour() uses the same orientation as array printing.
+        # Apply transform
         Z = kwargs.pop("transf")(Z)
+
+        # Need to transpose coz orientation is model.shape==(Nx, Ny),
+        # while contour() displays the same orientation as array printing.
         Z = Z.reshape(self.shape).T
 
         # Did we bother to specify set_over/set_under/set_bad ?
         has_out_of_range = getattr(kwargs["cmap"], "_rgba_over", None) is not None
 
-        # ax.imshow(Z[::-1])
+        # Unlike `ax.imshow(Z[::-1])`, `contourf` does not simply fill pixels/cells (but
+        # it does provide nice interpolation!) so there will be whitespace on the margins.
+        # No fix is needed, and anyway it would not be trivial/fast,
+        # ref https://github.com/matplotlib/basemap/issues/406 .
         collections = ax.contourf(
-            Z, **kwargs, extend="both" if has_out_of_range else "neither",
-            origin="lower", extent=(0, Lx, 0, Ly),
-            )
+            Z, **kwargs,
+            # origin=None,  # ⇒ NB: falsely stretches the field!!!
+            origin="lower",
+            extent=(0, Lx, 0, Ly),
+            extend="both" if has_out_of_range else "neither",
+        )
 
         # Contourf does not plot (at all) the bad regions. "Fake it" by facecolor
         if has_out_of_range: