Adding Reconstruction

.gitignore

@@ -4,10 +4,9 @@
 *.geqdsk
 *.pdf
 *.h5
+*.png
 *.svg
 *.jpg
 FreeGS.egg-info/
 build/
 dist/
-*.iml
-*.xml

15-TestMachineReconstruction.py

@@ -1,5 +1,5 @@
 """
-This script simuates an equilibrium using the forward solver in freegs,
+This script simulates an equilibrium using the forward solver in freegs,
 then uses the sensor measurments and coil currents to reconstruct the data
 """
 from freegs import machine, reconstruction
@@ -18,9 +18,9 @@
 pprime_order = 3
 ffprime_order = 3
 x_z1 = -0.6
-x_z2 = 0.6
-x_r1 = 1.1
-x_r2 = 1.1
+x_z2 = 0.7
+x_r1 = 1
+x_r2 = 1
 
 tokamak = machine.EfitTestMachine()
 
@@ -30,6 +30,6 @@
 print('Starting Reconstruction')
 
 eq_setup = {'tokamak': tokamak, 'Rmin': Rmin, 'Rmax': Rmax, 'Zmin': Zmin, 'Zmax': Zmax, 'nx': nx, 'ny': ny}
-Recon = reconstruction.Reconstruction(pprime_order, ffprime_order,use_VslCrnt=True, test=True, tolerance=1e-10,  **eq_setup)
+Recon = reconstruction.Reconstruction(pprime_order, ffprime_order, test=True, tolerance=1e-10,  **eq_setup)
 Recon.solve_from_tokamak()
 Recon.plot()

15-sensorMeasurements.py

@@ -0,0 +1,21 @@
+from freegs import machine, equilibrium, control, jtor, boundary, solve
+
+tokamak = machine.TestTokamakSensor()
+
+eq = equilibrium.Equilibrium(tokamak=tokamak,
+                        Rmin=0.1, Rmax=2.0,
+                        Zmin=-1.0, Zmax=1.0,
+                        nx=65, ny=65,
+                        boundary=boundary.freeBoundaryHagenow)
+
+profiles = jtor.ConstrainPaxisIp(eq, 1e3, 2e5, 2.0)
+
+xpoints = [(1.1, -0.6), (1.1, 0.8)]
+isoflux = [(1.1,-0.6, 1.1,0.6)]
+
+constrain = control.constrain(xpoints=xpoints, isoflux=isoflux)
+
+solve(eq, profiles, constrain, show=True)
+
+tokamak.printCurrents()
+tokamak.printMeasurements(eq=eq)

README.md

@@ -6,8 +6,6 @@ FreeGS: Free boundary Grad-Shafranov solver
 [![Build Status](https://github.com/bendudson/freegs/workflows/Tests/badge.svg)](https://github.com/bendudson/freegs/workflows/Tests/badge.svg)
 [![codecov](https://codecov.io/gh/bendudson/freegs/branch/master/graph/badge.svg?token=4dc6aHbu7K)](https://codecov.io/gh/bendudson/freegs)
 
-**NOTICE: This is no longer the main FreeGS repository** The official repository has moved to [https://github.com/freegs-plasma/freegs](https://github.com/freegs-plasma/freegs).
-
 This Python module calculates plasma equilibria for tokamak fusion experiments,
 by solving the Grad-Shafranov equation with free boundaries. Given a set of coils,
 plasma profiles and shape, FreeGS finds the currents in the coils which produce

freegs/coil.py

@@ -27,7 +27,7 @@
 from .gradshafranov import Greens, GreensBr, GreensBz, mu0
 import numpy as np
 import numbers
-from shapely.geometry import Point, Polygon
+from shapely.geometry import Point
 
 
 class AreaCurrentLimit:
@@ -206,7 +206,7 @@ def getForces(self, equilibrium):
             ]
         )  # Jphi x Br = - Fz
 
-    def inShape(self,polygon):
+    def inShape(self, polygon):
         if polygon.contains(Point(self.R, self.Z)):
             return 1
         else:

freegs/machine.py

@@ -179,11 +179,7 @@ def getForces(self, equilibrium):
         return forces
 
     def inShape(self, polygon):
-        counter = 0
-        for label, coil, multiplier in self.coils:
-            counter += coil.inShape(polygon) * multiplier
-
-        return counter
+        return sum(coil.inShape(polygon)*multiplier for label, coil, multiplier in self.coils)
 
     def __repr__(self):
         result = "Circuit(["

freegs/picard.py

@@ -189,10 +189,10 @@ def solve(
             # The user wants to wait for a limited plasma. The plasma is not limited.
             ok_to_break = False
 
-        # if show:
-        #     print('psi_relchange: '+str(psi_relchange))
-        #     print('bndry_relchange: '+str(bndry_relchange))
-        #     print('\n')
+        if show:
+            print('psi_relchange: '+str(psi_relchange))
+            print('bndry_relchange: '+str(bndry_relchange))
+            print('\n')
 
         # Check if the changes in psi are small enough and that it is ok to start checking for convergence
         if(((psi_maxchange < atol) or (psi_relchange < rtol)) and bndry_relchange < rtol and ok_to_break):

freegs/plotting.py

@@ -22,7 +22,6 @@
 
 from numpy import linspace, amin, amax
 from . import critical
-from . import machine
 
 
 def plotCoils(coils, axis=None):
@@ -136,4 +135,3 @@ def plotEquilibrium(eq, axis=None, show=True, oxpoints=True, wall=True, plot_sen
         plt.show()
 
     return axis
-

freegs/reconstruction.py

@@ -57,7 +57,7 @@ def __init__(self, pprime_order, ffprime_order, tolerance=1e-8, coil_weight=50,
 
         # Reconstruction attributes
         self.check_limited = True
-        self.pause = 0.0001
+        self.pause = 0.1
 
     # Methods for building measurment and uncertainty matrices
     def take_Measurements_from_dictionary(self, measurement_dict, sigma_dict):
@@ -157,8 +157,15 @@ def initialise_plasma_current(self):
         M_plasma = self.get_M_plasma()
 
         #Uses least squares solver to find coefficients and subsequently jtor
-        coefs = lstsq(self.Gplasma @ self.FEmatrix, M_plasma)[0]
-        jtor = self.FEmatrix @ coefs
+        if self.use_VslCrnt:
+            coefs = lstsq(np.block([self.Gplasma @ self.FEmatrix, self.Gvessel]),M_plasma)[0]
+            jtor = self.FEmatrix @ coefs[:-self.tokamak.n_basis]
+            self.tokamak.updateVesselCurrents(self.tokamak.eigenbasis@coefs[-self.tokamak.n_basis:])
+        else:
+            coefs = lstsq(self.Gplasma @ self.FEmatrix, M_plasma)[0]
+            jtor = self.FEmatrix @ coefs
+
+        print(coefs)
         return jtor
 
     # Calculating the plasma contirbution to the measurements
@@ -471,7 +478,7 @@ def get_fittingWeightVector(self):
         return np.diag([val**(-1) for val in self.sigma])
 
     # Calculate Finite Elements Grid
-    def get_FiniteElements(self, m=5, n=5):
+    def get_FiniteElements(self, m=3, n=3):
         """
         Parameters
         ----------
@@ -481,21 +488,20 @@ def get_FiniteElements(self, m=5, n=5):
         -------
         T - finite element basis matrix for current initialisation
         """
-        m += 2
-        n += 2
         R = self.R.flatten(order='F')
         Z = self.Z.flatten(order='F')
-        a, b, c, d = self.Rmin, self.Rmax, self.Zmin, self.Zmax
-        dR = (b - a) / (m - 1)
-        dZ = (d - c) / (n - 1)
-        r_FE = np.linspace(a, b, m)
-        z_FE = np.linspace(c, d, n)
-
-        T = np.zeros((self.nx * self.ny, (m - 2) * (n - 2)))
+        Rdif, Zdif = (self.Rmax-self.Rmin)/3, (self.Zmax-self.Zmin)/3
+        a, b, c, d = self.Rmin+Rdif, self.Rmin+2*Rdif, self.Zmin+Zdif, self.Zmin+2*Zdif
+        dR = (b - a) / (m + 1)
+        dZ = (d - c) / (n + 1)
+        r_FE = np.linspace(a, b, m+2)
+        z_FE = np.linspace(c, d, n+2)
+
+        T = np.zeros((self.nx * self.ny, m * n))
         row_num = 0
-        for i in range(1, m - 1):
+        for i in range(1, m +1):
             r_h = r_FE[i]
-            for j in range(1, n - 1):
+            for j in range(1, n +1):
                 z_h = z_FE[j]
 
                 for h in range(T.shape[0]):
@@ -504,7 +510,6 @@ def get_FiniteElements(self, m=5, n=5):
                             T[h, row_num] = (1 - abs(R[h] - r_h) / dR) * (
                                         1 - abs(Z[h] - z_h) / dZ)
                 row_num += 1
-
         return T
 
     # Indexing function
@@ -563,6 +568,7 @@ def reconstruct(self):
 
                 # Performs plasma current density initialisation
                 if self.use_CrntInit:
+                    print('yes')
                     jtor_1d = self.initialise_plasma_current()
                     self.Jtor = np.reshape(jtor_1d, (self.nx, self.ny),order='F')
 
@@ -718,7 +724,6 @@ def check_all_measurements_used(self, measurement_dict):
                 print(key)
                 print('Not all inputted Measurement values used.')
 
-
     # Methods for solving
     def solve_from_tokamak(self):
         self.take_Measurements_from_tokamak()

freegs/shaped_coil.py

@@ -26,7 +26,7 @@
 from . import quadrature
 from . import polygons
 from .gradshafranov import Greens, GreensBr, GreensBz
-from shapely.geometry import Point, Polygon
+from shapely.geometry import Polygon
 
 import numpy as np
 
@@ -118,15 +118,15 @@ def controlBz(self, R, Z):
             result += GreensBz(R_fil, Z_fil, R, Z) * weight
         return result
 
+    def inShape(self,polygon):
+        Shaped_Coil = Polygon([shape for shape in self.shape])
+        return (polygon.intersection(Shaped_Coil).area) / (self._area)
+
     def __repr__(self):
         return "ShapedCoil({0}, current={1:.1f}, turns={2}, control={3})".format(
             self.shape, self.current, self.turns, self.control
         )
 
-    def inShape(self,polygon):
-        Shaped_Coil = Polygon([shape for shape in self.shape])
-        return (polygon.intersection(Shaped_Coil).area) / (self._area)
-
     @property
     def R(self):
         """

freegs/test_reconstruction.py

@@ -1,4 +1,4 @@
-from freegs import machine, equilibrium, reconstruction, plotting, critical
+from freegs import machine, reconstruction
 import numpy as np
 import pickle
 
@@ -143,4 +143,5 @@ def test_vessel_eigenmode():
     tokamak = machine.EfitTestMachine()
     Recon = reconstruction.Reconstruction(pprime_order, ffprime_order, show=show, **eq_setup)
     Recon.solve_from_dictionary(measurement_dict=measurement_dict, sigma_dict=sigma_dict)
-    assert abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i]) > 50*abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i+1]) and abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i]) > 50*abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i-1])
+    assert abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i]) > 50*abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i+1]) and abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i]) > 50*abs(Recon.coefs[-tokamak.eigenbasis.shape[1]+i-1])
+