Merge pull request #93 from FusionComputingLab/fix-normal-mode-greens

Correct normal-mode Green transformation

Author: kpentland

examples/example10 - growth_rates.ipynb

@@ -337,7 +337,17 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can visualise the poloidal flux produced by the unstable mode using the following code. We can also do this for other (stable) modes if required (just change the index `i` in the cell below)."
+    "We can visualise the poloidal flux produced by the unstable mode using the following code. We can also do this for other (stable) modes if required (just change the index `i` in the cell below).\n",
+    "\n",
+    "The physical metal currents $\\vec{I}_m$ and normal-mode currents $\\vec{I}_d$ are related by\n",
+    "\n",
+    "$$\\vec{I}_m = P\\vec{I}_d.$$\n",
+    "\n",
+    "Writing the physical-metal Green maps as $G_m$, the flux is therefore\n",
+    "\n",
+    "$$\\psi = \\vec{I}_m^T G_m = \\vec{I}_d^T P^T G_m.$$\n",
+    "\n",
+    "Thus the Green maps per unit normal-mode current are $G_d=P^T G_m$. Notice that this is a covector transformation: $P^{-1}$ is used to recover modal currents from physical currents, but it is not the transformation for the Green maps."
    ]
   },
   {
@@ -348,16 +358,16 @@
    "source": [
     "# mode number (choose which one you want to visualise)\n",
     "i = idx # default is unstable mode\n",
-    "mode_currents = np.real(modes[:,i])\n",
+    "mode_currents = np.real(modes[:-1, i])\n",
     "\n",
     "# the associated instability timescale and growth rate\n",
     "print(f\"Mode {i} ---> {'stable' if np.real(timescales[i]) < 0 else 'unstable'}\")\n",
     "print(f\"Growth rate = {np.real(growth_rates[i]):.2e} [1/s]\")\n",
     "print(f\"Timescale = {np.real(timescales[i]):.2e} [s]\")\n",
     "\n",
-    "# multiply each metal current (from the eigenvector) with its corresponding Greens matrix and sum\n",
-    "# (don't forget to omit the plasma current mode, i.e. the final element)\n",
-    "flux = np.sum(mode_currents[0:-1, np.newaxis, np.newaxis] * nonlinear_solver.vessel_modes_greens[i], axis=0)\n",
+    "# Contract the modal currents with the modal Green maps G_d = P.T @ G_m.\n",
+    "# The final eigenvector element, representing plasma current, was omitted above.\n",
+    "flux = np.einsum(\"i,ijk->jk\", mode_currents, nonlinear_solver.vessel_modes_greens)\n",
     "\n",
     "# plot\n",
     "fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n",
@@ -743,7 +753,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Again, we can plot how these look using similar code as before, except we need to remember that now we've done a mode decomposition, the eigenmode currents $\\vec{I}_d$ need to be transformed back into the original metal currents $\\vec{I}_m$. This is done in the plotting cell below."
+    "Again, we can plot how these look using similar code as before. Here we explicitly transform the retained normal-mode currents back to physical metal currents using $\\vec{I}_m=P\\vec{I}_d$, then contract those currents with the physical-metal Green maps $G_m$. This is equivalent to contracting $\\vec{I}_d$ with the corresponding retained rows of $G_d=P^T G_m$."
    ]
   },
   {
@@ -784,15 +794,17 @@
     "\n",
     "# this function transforms the (decomposed) mode currents back to regular metal currents\n",
     "# (don't forget to omit the plasma current mode, i.e. the final element)\n",
-    "mode_currents = nonlinear_solver_option_1.evol_metal_curr.IdtoIvessel(np.real(modes[0:-1,i]))\n",
+    "mode_currents = nonlinear_solver_option_1.evol_metal_curr.IdtoIvessel(\n",
+    "    np.real(modes[:-1, i])\n",
+    ")\n",
     "\n",
     "# the associated instability timescale and growth rate\n",
     "print(f\"Mode {i} ---> {'stable' if np.real(timescales[i]) < 0 else 'unstable'}\")\n",
     "print(f\"Growth rate = {np.real(growth_rates[i]):.2e} [1/s]\")\n",
     "print(f\"Timescale = {np.real(timescales[i]):.2e} [s]\")\n",
     "\n",
-    "# multiply each metal current (from the eigenvector) with its corresponding Greens matrix and sum\n",
-    "flux = np.sum(mode_currents[:, np.newaxis, np.newaxis] * nonlinear_solver_option_1.vessel_modes_greens[i], axis=0)\n",
+    "# Contract reconstructed physical currents with the physical-metal Green maps.\n",
+    "flux = np.einsum(\"i,ijk->jk\", mode_currents, eq._vgreen)\n",
     "\n",
     "# plot\n",
     "fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n",

freegsnke/normal_modes.py

@@ -114,8 +114,21 @@ def __init__(self, coil_resist, coil_self_ind, n_coils, n_active_coils):
 
     def normal_modes_greens(self, eq_vgreen):
         """
-        Calculates the green functions of the vessel normal modes,
-        i.e. the psi flux per unit current for each mode.
+        Calculate the Green functions of the vessel normal modes.
+
+        If ``I_m`` contains the physical metal currents and ``I_d`` contains
+        the normal-mode currents, the current transformation is
+
+        ``I_m = Pmatrix @ I_d``.
+
+        Since the physical-coil Green functions satisfy
+        ``psi = I_m.T @ G_m``, substitution gives
+
+        ``psi = I_d.T @ Pmatrix.T @ G_m``.
+
+        The modal Green functions are therefore ``Pmatrix.T @ G_m``.  The
+        inverse transformation is instead used to recover modal currents from
+        physical currents and must not be applied to the Green functions.
 
         Parameters
         ----------
@@ -124,10 +137,6 @@ def normal_modes_greens(self, eq_vgreen):
             Can be found at eq._vgreen. np.shape(eq_vgreen)=(n_coils, nx, ny)
         """
 
-        dgreen = np.sum(
-            eq_vgreen[np.newaxis, :, :, :]
-            * self.Pmatrix_inverse[:, :, np.newaxis, np.newaxis],
-            axis=1,
-        )
-
-        return dgreen
+        grid_shape = eq_vgreen.shape[1:]
+        physical_greens = eq_vgreen.reshape(self.n_coils, -1)
+        return (self.Pmatrix.T @ physical_greens).reshape(self.n_coils, *grid_shape)

freegsnke/tests/test_dynamics.py

@@ -96,7 +96,7 @@ def create_machine():
 
 def test_linearised_growth_rate(create_machine):
     tokamak, eq, profiles, stepping = create_machine
-    true_GR = 0.0586
+    true_GR = 0.05867
     # check that
     assert (
         abs((stepping.linearised_sol.instability_timescale[0] - true_GR) / true_GR)

freegsnke/tests/test_normal_modes.py

@@ -0,0 +1,67 @@
+"""Tests for passive-structure normal-mode transformations."""
+
+import numpy as np
+
+from freegsnke.normal_modes import mode_decomposition
+
+
+def _mode_decomposition():
+    """Return a small decomposition whose passive modes are not orthogonal."""
+
+    resistances = np.array([0.8, 1.1, 1.6, 2.0])
+    inductances = np.array(
+        [
+            [4.0, 0.2, 0.1, 0.0],
+            [0.2, 3.2, 0.5, 0.1],
+            [0.1, 0.5, 2.8, 0.4],
+            [0.0, 0.1, 0.4, 2.5],
+        ]
+    )
+    return mode_decomposition(
+        coil_resist=resistances,
+        coil_self_ind=inductances,
+        n_coils=4,
+        n_active_coils=1,
+    )
+
+
+def test_current_transform_uses_inverse_not_transpose():
+    """The non-orthogonal mode basis requires its inverse for currents."""
+
+    modes = _mode_decomposition()
+    modal_currents = np.array([1.0, -0.4, 0.7, 0.2])
+    physical_currents = modes.Pmatrix @ modal_currents
+
+    assert not np.allclose(modes.Pmatrix.T, modes.Pmatrix_inverse)
+    np.testing.assert_allclose(
+        modes.Pmatrix_inverse @ physical_currents,
+        modal_currents,
+        rtol=1e-13,
+        atol=1e-13,
+    )
+
+
+def test_mode_greens_reproduce_physical_current_flux():
+    """Modal and reconstructed physical currents must produce the same flux."""
+
+    modes = _mode_decomposition()
+    rng = np.random.default_rng(42)
+    physical_greens = rng.normal(size=(modes.n_coils, 5, 6))
+    modal_currents = np.array([1.0, -0.4, 0.7, 0.2])
+
+    physical_currents = modes.Pmatrix @ modal_currents
+    flux_from_physical_currents = np.einsum(
+        "i,ijk->jk", physical_currents, physical_greens
+    )
+    flux_from_modal_currents = np.einsum(
+        "i,ijk->jk",
+        modal_currents,
+        modes.normal_modes_greens(physical_greens),
+    )
+
+    np.testing.assert_allclose(
+        flux_from_modal_currents,
+        flux_from_physical_currents,
+        rtol=1e-13,
+        atol=1e-13,
+    )