Merge branch 'main' into flux_averaging_example

Author: georgeholt1

examples/example05 - evolutive_forward_solve.ipynb

@@ -244,6 +244,44 @@
     ")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Set active coil voltages\n",
+    "\n",
+    "In this example, we will evolve a plasma in absence of any control policy or current drive.\n",
+    "\n",
+    "Just as an example, the following calculates active voltages to be applied to the poloidal field coils (and Solenoid) using $V = RI$, with current values as defined by the initial equilibrium (i.e. we have constant voltages).\n",
+    "\n",
+    "In most FreeGSNKE use cases, these active voltages will be determined by a control policy."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.coil_resist)[:stepping.evol_metal_curr.n_active_coils] "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To start, the solver is prepared by setting the initial conditions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stepping.initialize_from_ICs(eq, profiles)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -293,7 +331,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.R)[:stepping.evol_metal_curr.n_active_coils] "
+    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.coil_resist)[:stepping.evol_metal_curr.n_active_coils] "
    ]
   },
   {
@@ -443,21 +481,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [
-    {
-     "ename": "NameError",
-     "evalue": "name 'plt' is not defined",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn[1], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Plot evolution of tracked values and compare between linear and non-linear evolution\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m fig, axs \u001b[38;5;241m=\u001b[39m \u001b[43mplt\u001b[49m\u001b[38;5;241m.\u001b[39msubplots(\u001b[38;5;241m3\u001b[39m, \u001b[38;5;241m3\u001b[39m, figsize\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m15\u001b[39m, \u001b[38;5;241m10\u001b[39m), dpi\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m80\u001b[39m, constrained_layout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m      3\u001b[0m axs_flat \u001b[38;5;241m=\u001b[39m axs\u001b[38;5;241m.\u001b[39mflat\n\u001b[1;32m      5\u001b[0m axs_flat[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mplot(history_times, history_o_points[:, \u001b[38;5;241m0\u001b[39m],\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mk+\u001b[39m\u001b[38;5;124m'\u001b[39m, label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlinear\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
-      "\u001b[0;31mNameError\u001b[0m: name 'plt' is not defined"
-     ]
-    }
-   ],
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "# Plot evolution of tracked values and compare between linear and non-linear evolution\n",
     "fig, axs = plt.subplots(3, 3, figsize=(15, 10), dpi=80, constrained_layout=True)\n",
@@ -553,7 +579,7 @@
    "source": [
     "## Limited plasma example\n",
     "\n",
-    "In this example, we examine a limiter configuration plasma and show how to use time-varying plasma profile parameters and plasma resisitivity."
+    "In this example, we examine a limiter configuration plasma and show how to use time-varying plasma profile parameters and plasma resisitivity. Note there is also the option to specify time-dependent active coil resistances. For example, if you want to switch a coil \"off\" during a simulation you can set its resistance to an arbitrarily high number. See the optional parameter `custom_active_coil_resistances`. "
    ]
   },
   {
@@ -648,7 +674,7 @@
    "outputs": [],
    "source": [
     "# determine the (constant) active voltages to be applied at each time step\n",
-    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.R)[:stepping.evol_metal_curr.n_active_coils] "
+    "voltages = (stepping.vessel_currents_vec*stepping.evol_metal_curr.coil_resist)[:stepping.evol_metal_curr.n_active_coils] "
    ]
   },
   {
@@ -782,6 +808,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# recalculate the active voltages using the new currents\n",
+    "U_active = (stepping.vessel_currents_vec*stepping.evol_metal_curr.coil_resist)[:stepping.evol_metal_curr.n_active_coils] \n",
+    "\n",
+    "# number of time steps to simulate\n",
+    "max_count = 50\n",
+    "\n",
+    "# reset the solver object by resetting the initial conditions\n",
+    "stepping.initialize_from_ICs(eq, profiles)\n",
+    "\n",
     "# initialising some variables for iteration and logging\n",
     "counter = 0\n",
     "t = 0\n",

freegsnke/circuit_eq_metal.py

@@ -98,7 +98,8 @@ def __init__(
             self.coil_resist = eq.tokamak.coil_resist
 
         self.Rm1 = 1.0 / self.coil_resist
-        self.R = self.coil_resist
+        # self.R = np.copy(self.coil_resist)
+        self.active_coil_resistances = np.copy(self.coil_resist[: self.n_active_coils])
 
         # prepare inductance data
         if coil_self_ind is not None:
@@ -110,6 +111,8 @@ def __init__(
         else:
             self.coil_self_ind = eq.tokamak.coil_self_ind
 
+        self.build_rm1l()
+
         self.flag_vessel_eig = flag_vessel_eig
         self.flag_plasma = flag_plasma
 
@@ -138,6 +141,12 @@ def __init__(
         # Dummy voltage vector
         self.empty_U = np.zeros(self.n_coils)
 
+    def build_rm1l(
+        self,
+    ):
+        # active + passive
+        self.rm1l_non_symm = np.diag(self.coil_resist**-1.0) @ self.coil_self_ind
+
     def make_selected_mode_mask(self, mode_coupling_masks, verbose):
         """Creates a mask for the vessel normal modes to include in the circuit
         equations, based on the maximum frequency of the selected modes.
@@ -243,7 +252,7 @@ def initialize_for_eig(
         # diagonalised separately from the active coils. The modes of used for the passive structures
         # diagonalise the isolated dynamics of the walls.
         # Equation is Lambda**(-1)Iddot + I = F
-        self.Lambdam1 = self.Pm1 @ (self.normal_modes.rm1l_non_symm @ self.P)
+        self.Lambdam1 = self.Pm1 @ (self.rm1l_non_symm @ self.P)
         # self.RP = np.diag(self.coil_resist) @ self.P
         # self.RP_inv = np.linalg.solve(self.RP.T @ self.RP, self.RP.T)
         # self.Lambdam1 = (self.RP_inv @ self.coil_self_ind) @ self.P
@@ -260,6 +269,15 @@ def initialize_for_eig(
         else:
             self.forcing_term = self.forcing_term_eig_no_plasma
 
+    def reset_active_coil_resistances(self, active_coil_resistances):
+        self.coil_resist = np.concatenate(
+            (active_coil_resistances, self.coil_resist[self.n_active_coils :])
+        )
+        self.active_coil_resistances = np.copy(self.coil_resist[: self.n_active_coils])
+        self.Rm1 = 1 / self.coil_resist
+        self.build_rm1l()
+        self.Lambdam1 = self.Pm1 @ (self.rm1l_non_symm @ self.P)
+
     def initialize_for_no_eig(self):
         """Initializes the metal currents object for the case where vessel
         eigenmodes are not used."""

freegsnke/linear_solve.py

@@ -34,7 +34,7 @@ class linear_solver:
 
     def __init__(
         self,
-        eq,
+        coil_numbers,
         Lambdam1,
         P,
         Pm1,
@@ -93,10 +93,10 @@ def __init__(
         self.Pm1Rm1Mey = np.matmul(self.Pm1Rm1, Mey)
         self.MyeP = np.matmul(Mey.T, P).T
 
-        if Lambdam1 is None:
-            self.Lambdam1 = Pm1 @ (Rm1 @ (eq.tokamak.coil_self_ind @ P))
-        else:
-            self.Lambdam1 = Lambdam1
+        # if Lambdam1 is None:
+        #     self.Lambdam1 = Pm1 @ (Rm1 @ (eq.tokamak.coil_self_ind @ P))
+        # else:
+        self.Lambdam1 = Lambdam1
         self.n_independent_vars = np.shape(self.Lambdam1)[0]
 
         self.Mmatrix = np.zeros(
@@ -109,7 +109,7 @@ def __init__(
             (self.n_independent_vars + 1, self.n_independent_vars + 1)
         )
 
-        self.n_active_coils = eq.tokamak.n_active_coils
+        self.n_active_coils, self.n_coils = coil_numbers
 
         self.solver = implicit_euler_solver(
             Mmatrix=np.eye(self.n_independent_vars + 1),