Small patches

freegs4e/critical.py

@@ -410,9 +410,11 @@ def scan_for_crit(R, Z, psi):
                     )  # + 0.5*(fRR*delta_R**2 + fZZ*delta_Z**2 + fRZ*delta_R*delta_Z)
                     crpoint = (R0 + delta_R, Z0 + delta_Z, est_psi)
                     if det > 0.0:
-                        opoint = [crpoint] + opoint
+                        # opoint = [crpoint] + opoint
+                        opoint.append(crpoint)
                     else:
-                        xpoint = [crpoint] + xpoint
+                        # xpoint = [crpoint] + xpoint
+                        xpoint.append(crpoint)
 
     xpoint = np.array(xpoint)
     opoint = np.array(opoint)
@@ -859,6 +861,16 @@ def inside_mask(
     if use_geom:
         # cure flooding
         mask = mask * geom_inside_mask(R, Z, opoint, xpoint)
+        # apply geometric masking criterion to second Xpoint if close to double null
+        if len(xpoint > 1):
+            if (
+                np.abs(
+                    (xpoint[0, 2] - xpoint[1, 2])
+                    / (opoint[0, 2] - xpoint[0, 2])
+                )
+                < 0.1
+            ):
+                mask = mask * geom_inside_mask(R, Z, opoint, xpoint[1:])
     return mask
 
 

freegs4e/equilibrium.py

@@ -564,7 +564,7 @@ def psi(self):
             The total poloidal flux due to the plasma and the coils [Webers/2pi].
         """
 
-        return self.plasma_psi + self.tokamak.calcPsiFromGreens(self._pgreen)
+        return self.plasma_psi + self.tokamak.getPsitokamak(self._vgreen)
 
     def psiRZ(self, R, Z):
         """
@@ -2272,7 +2272,7 @@ def normalised_total_Beta(self):
         """
         Calculates the total beta from the following definition:
 
-            normalised_total_Beta = ( (1 / poloidalBeta2) + (1/toroidalBeta) )^(-1).
+            normalised_total_Beta = ( (1 / poloidalBeta1) + (1/toroidalBeta1) )^(-1).
 
         Parameters
         ----------
@@ -2285,7 +2285,7 @@ def normalised_total_Beta(self):
         """
 
         return 1.0 / (
-            (1.0 / self.poloidalBeta()) + (1.0 / self.toroidalBeta())
+            (1.0 / self.poloidalBeta1()) + (1.0 / self.toroidalBeta1())
         )
 
     def strikepoints(

freegs4e/jtor.py

@@ -184,12 +184,14 @@ def Jtor_part1(self, R, Z, psi, psi_bndry=None, mask_outside_limiter=None):
             if hasattr(self, "diverted_core_mask"):
                 if self.diverted_core_mask is not None:
                     previous_core_size = np.sum(self.diverted_core_mask)
+                    skipped_xpts = 0
                     # check size change
                     check = (
                         np.sum(diverted_core_mask) / previous_core_size < 0.5
                     )
+                    # check there's more candidates
                     check *= len(xpt) > 1
-                    if check:
+                    while check:
                         # try using second xpt as primary xpt
                         alt_diverted_core_mask = critical.inside_mask(
                             R,
@@ -205,12 +207,18 @@ def Jtor_part1(self, R, Z, psi, psi_bndry=None, mask_outside_limiter=None):
                             self.edge_mask * alt_diverted_core_mask
                         )
                         if edge_pixels == 0:
-                            # print(
-                            #     "Discarding 'primary' Xpoint! Please check final result"
-                            # )
+                            # the candidate is valid
                             xpt = xpt[1:]
                             psi_bndry = xpt[1, 2]
                             diverted_core_mask = alt_diverted_core_mask.copy()
+
+                            # check if there could be better candidates
+                            check = (
+                                np.sum(diverted_core_mask) / previous_core_size
+                                < 0.5
+                            )
+                            # check there's more candidates
+                            check *= len(xpt) > 1
         else:
             # No X-points
             psi_bndry = psi[0, 0]
@@ -353,6 +361,8 @@ def Jtor_part2(self, R, Z, psi, psi_axis, psi_bndry, mask):
         IR = (
             np.sum(jtorshape * R / self.Raxis) * dR * dZ
         )  # romb(romb(jtorshape * R / self.Raxis)) * dR * dZ
+        if IR == 0:
+            raise ValueError("No core mask!")
         I_R = (
             np.sum(jtorshape * self.Raxis / R) * dR * dZ
         )  # romb(romb(jtorshape * self.Raxis / R)) * dR * dZ
@@ -392,6 +402,11 @@ def pprime(self, pn):
         """
 
         shape = (1.0 - np.clip(pn, 0.0, 1.0) ** self.alpha_m) ** self.alpha_n
+        if hasattr(self, "Beta0") is False:
+            raise ValueError(
+                "The pprime profile cannot be normalised. "
+                "Please first calculate Jtor for this profile. "
+            )
 
         return self.L * self.Beta0 / self.Raxis * shape
 
@@ -412,6 +427,11 @@ def ffprime(self, pn):
         """
 
         shape = (1.0 - np.clip(pn, 0.0, 1.0) ** self.alpha_m) ** self.alpha_n
+        if hasattr(self, "Beta0") is False:
+            raise ValueError(
+                "The ffprime profile cannot be normalised. "
+                "Please first calculate Jtor for this profile. "
+            )
 
         return mu0 * self.L * (1 - self.Beta0) * self.Raxis * shape
 
@@ -579,11 +599,15 @@ def pprime(self, pn):
         """
 
         shape = (1.0 - np.clip(pn, 0.0, 1.0) ** self.alpha_m) ** self.alpha_n
-        if hasattr(self, "L") is False:
-            self.L = 1
-            print(
-                "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
+        if hasattr(self, "Beta0") is False:
+            raise ValueError(
+                "The pprime profile cannot be normalised. "
+                "Please first calculate Jtor for this profile. "
             )
+            # self.L = 1
+            # print(
+            #     "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
+            # )
         return self.L * self.Beta0 / self.Raxis * shape
 
     def ffprime(self, pn):
@@ -603,10 +627,10 @@ def ffprime(self, pn):
         """
 
         shape = (1.0 - np.clip(pn, 0.0, 1.0) ** self.alpha_m) ** self.alpha_n
-        if hasattr(self, "L") is False:
-            self.L = 1
-            print(
-                "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
+        if hasattr(self, "Beta0") is False:
+            raise ValueError(
+                "The ffprime profile cannot be normalised. "
+                "Please first calculate Jtor for this profile. "
             )
         return mu0 * self.L * (1 - self.Beta0) * self.Raxis * shape
 
@@ -753,10 +777,10 @@ def pprime(self, pn):
         """
 
         shape = (1.0 - np.clip(pn, 0.0, 1.0) ** self.alpha_m) ** self.alpha_n
-        if hasattr(self, "L") is False:
-            self.L = 1
-            print(
-                "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
+        if hasattr(self, "Beta0") is False:
+            raise ValueError(
+                "The pprime profile cannot be normalised. "
+                "Please first calculate Jtor for this profile. "
             )
         return self.L * self.Beta0 / self.Raxis * shape
 
@@ -777,10 +801,10 @@ def ffprime(self, pn):
         """
 
         shape = (1.0 - np.clip(pn, 0.0, 1.0) ** self.alpha_m) ** self.alpha_n
-        if hasattr(self, "L") is False:
-            self.L = 1
-            print(
-                "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
+        if hasattr(self, "Beta0") is False:
+            raise ValueError(
+                "The ffprime profile cannot be normalised. "
+                "Please first calculate Jtor for this profile. "
             )
         return mu0 * self.L * (1 - self.Beta0) * self.Raxis * shape
 
@@ -881,6 +905,27 @@ def initialize_profile(self):
             self.beta = np.concatenate((self.beta, [-np.sum(self.beta)]))
         self.beta_exp = np.arange(0, len(self.beta))
 
+    def build_dJtorpsin1(
+        self,
+    ):
+        # calculate dJ/dpsi_n at psi_n=1
+        pprime_term = (
+            self.alpha[1:, np.newaxis, np.newaxis]
+            * self.alpha_exp[1:, np.newaxis, np.newaxis]
+        )
+        pprime_term = np.sum(pprime_term, axis=0)
+        pprime_term *= self.eqR / self.Raxis
+
+        ffprime_term = (
+            self.beta[1:, np.newaxis, np.newaxis]
+            * self.beta_exp[1:, np.newaxis, np.newaxis]
+        )
+        ffprime_term = np.sum(ffprime_term, axis=0)
+        ffprime_term *= self.Raxis / self.eqR
+        ffprime_term /= mu0
+
+        self.dJtorpsin1 = pprime_term + ffprime_term
+
     def Jtor_part2(
         self,
         R,
@@ -966,16 +1011,44 @@ def Jtor_part2(
         # sum together
         Jtor = pprime_term + ffprime_term
 
+        # calculate dJ/dpsi_n
+        pprime_term = (
+            psi_norm[np.newaxis, :, :]
+            ** self.alpha_exp[:-1, np.newaxis, np.newaxis]
+        )
+        pprime_term *= (
+            self.alpha[1:, np.newaxis, np.newaxis]
+            * self.alpha_exp[1:, np.newaxis, np.newaxis]
+        )
+        pprime_term = np.sum(pprime_term, axis=0)
+        pprime_term *= R / self.Raxis
+
+        ffprime_term = (
+            psi_norm[np.newaxis, :, :]
+            ** self.beta_exp[:-1, np.newaxis, np.newaxis]
+        )
+        ffprime_term *= (
+            self.beta[1:, np.newaxis, np.newaxis]
+            * self.beta_exp[1:, np.newaxis, np.newaxis]
+        )
+        ffprime_term = np.sum(ffprime_term, axis=0)
+        ffprime_term *= self.Raxis / R
+        ffprime_term /= mu0
+
+        dJtordpsin = pprime_term + ffprime_term
+        self.dJtordpsi = dJtordpsin / (psi_axis - psi_bndry)
+
         # put to zero all current outside the LCFS
-        Jtor *= psi > psi_bndry
+        # Jtor *= psi > psi_bndry
 
-        Jtor *= self.Ip * Jtor > 0
+        # Jtor *= self.Ip * Jtor > 0
 
         if torefine:
             return Jtor
 
         if mask is not None:
             Jtor *= mask
+            self.dJtordpsi *= mask
 
         # if Ip normalisation is required, do it
         if self.Ip_logic:
@@ -987,6 +1060,8 @@ def Jtor_part2(
                 )
             L = self.Ip / (jtorIp * dR * dZ)
             Jtor = L * Jtor
+            self.dJtordpsi *= L
+
         else:
             L = 1.0
 
@@ -1022,11 +1097,12 @@ def pprime(self, pn):
             list(np.shape(self.alpha)) + [1] * len(shape_pn)
         )
         shape = np.sum(shape, axis=0)
-        if hasattr(self, "L") is False:
-            self.L = 1
-            print(
-                "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
-            )
+        if self.Ip_logic is True:
+            if hasattr(self, "L") is False:
+                raise ValueError(
+                    "The pprime profile cannot be normalised. "
+                    "Please first calculate Jtor for this profile. "
+                )
         return self.L * shape / self.Raxis
 
     def ffprime(self, pn):
@@ -1055,11 +1131,12 @@ def ffprime(self, pn):
             list(np.shape(self.beta)) + [1] * len(shape_pn)
         )
         shape = np.sum(shape, axis=0)
-        if hasattr(self, "L") is False:
-            self.L = 1
-            print(
-                "This is using self.L=1, which is likely not appropriate. Please calculate Jtor first to ensure the correct normalization."
-            )
+        if self.Ip_logic is True:
+            if hasattr(self, "L") is False:
+                raise ValueError(
+                    "The ffprime profile cannot be normalised. "
+                    "Please first calculate Jtor for this profile. "
+                )
         return self.L * shape * self.Raxis
 
     def pressure(self, pn):