From 562f9a65b20cfa72eff0b6f42db919fa67ecbdd7 Mon Sep 17 00:00:00 2001 From: Nico Amorisco Date: Wed, 4 Jun 2025 10:15:19 +0100 Subject: [PATCH 1/7] added geometrical mask exclusion to second xpoint for cases close to double null --- freegs4e/critical.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/freegs4e/critical.py b/freegs4e/critical.py index c840617..7af27c3 100644 --- a/freegs4e/critical.py +++ b/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,10 @@ 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]))<.1: + mask = mask * geom_inside_mask(R, Z, opoint, xpoint[1:]) return mask From d01809d40839c09ee9746972d2550fe28e719ac3 Mon Sep 17 00:00:00 2001 From: Nico Amorisco Date: Wed, 4 Jun 2025 10:16:54 +0100 Subject: [PATCH 2/7] added valueerror for cases with no plasma domain --- freegs4e/jtor.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/freegs4e/jtor.py b/freegs4e/jtor.py index 96dff25..b939ddd 100644 --- a/freegs4e/jtor.py +++ b/freegs4e/jtor.py @@ -353,6 +353,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 From 3696082cfa6521f8b00e26e57dd9e725ca8bf24a Mon Sep 17 00:00:00 2001 From: Nico Amorisco Date: Wed, 4 Jun 2025 10:23:05 +0100 Subject: [PATCH 3/7] black --- freegs4e/critical.py | 10 ++++++++-- freegs4e/jtor.py | 2 +- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/freegs4e/critical.py b/freegs4e/critical.py index 7af27c3..0228d92 100644 --- a/freegs4e/critical.py +++ b/freegs4e/critical.py @@ -862,8 +862,14 @@ def inside_mask( # 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]))<.1: + 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 diff --git a/freegs4e/jtor.py b/freegs4e/jtor.py index 5a7a074..90c5244 100644 --- a/freegs4e/jtor.py +++ b/freegs4e/jtor.py @@ -353,7 +353,7 @@ 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: + if IR == 0: raise ValueError("No core mask!") I_R = ( np.sum(jtorshape * self.Raxis / R) * dR * dZ From 3212297881e91a4fc6c9073a4dc5729edce26392 Mon Sep 17 00:00:00 2001 From: kpentland Date: Sat, 21 Jun 2025 13:24:44 +0100 Subject: [PATCH 4/7] minor typo in normalised total beta --- freegs4e/equilibrium.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index 1a3fbf7..83d1486 100644 --- a/freegs4e/equilibrium.py +++ b/freegs4e/equilibrium.py @@ -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( From 85e007c72f2ae76f9e5a96396928e9d16a814942 Mon Sep 17 00:00:00 2001 From: Nico Amorisco Date: Wed, 9 Jul 2025 11:33:21 +0100 Subject: [PATCH 5/7] using vectorised psi_tokamak in eq.psi() --- freegs4e/equilibrium.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index ee5c2b5..0fdebc2 100644 --- a/freegs4e/equilibrium.py +++ b/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): """ From cdfa4ad7190b69c281da8e1ab8aadc76eaf176d6 Mon Sep 17 00:00:00 2001 From: Nico Amorisco Date: Wed, 9 Jul 2025 11:35:01 +0100 Subject: [PATCH 6/7] introduced dJdpsi calc and checks for normalization of profiles --- freegs4e/jtor.py | 119 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 83 insertions(+), 36 deletions(-) diff --git a/freegs4e/jtor.py b/freegs4e/jtor.py index 90c5244..67a4860 100644 --- a/freegs4e/jtor.py +++ b/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,17 @@ 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] @@ -394,6 +401,9 @@ 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 @@ -414,6 +424,9 @@ 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 @@ -581,11 +594,13 @@ 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): @@ -605,11 +620,9 @@ 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 def fvac(self): @@ -755,11 +768,9 @@ 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 def ffprime(self, pn): @@ -779,11 +790,9 @@ 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 def fvac(self): @@ -883,6 +892,22 @@ 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, @@ -968,16 +993,38 @@ 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: @@ -989,6 +1036,8 @@ def Jtor_part2( ) L = self.Ip / (jtorIp * dR * dZ) Jtor = L * Jtor + self.dJtordpsi *= L + else: L = 1.0 @@ -1024,11 +1073,10 @@ 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): @@ -1057,11 +1105,10 @@ 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): From 903f99c3efa2deaf7db4148eae975864fea220e5 Mon Sep 17 00:00:00 2001 From: Nico Amorisco Date: Wed, 9 Jul 2025 11:38:23 +0100 Subject: [PATCH 7/7] black --- freegs4e/jtor.py | 82 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 55 insertions(+), 27 deletions(-) diff --git a/freegs4e/jtor.py b/freegs4e/jtor.py index 67a4860..e2837a8 100644 --- a/freegs4e/jtor.py +++ b/freegs4e/jtor.py @@ -214,7 +214,8 @@ def Jtor_part1(self, R, Z, psi, psi_bndry=None, mask_outside_limiter=None): # check if there could be better candidates check = ( - np.sum(diverted_core_mask) / previous_core_size < 0.5 + np.sum(diverted_core_mask) / previous_core_size + < 0.5 ) # check there's more candidates check *= len(xpt) > 1 @@ -402,8 +403,10 @@ 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. ") + raise ValueError( + "The pprime profile cannot be normalised. " + "Please first calculate Jtor for this profile. " + ) return self.L * self.Beta0 / self.Raxis * shape @@ -425,8 +428,10 @@ 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. ") + 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 @@ -595,8 +600,10 @@ 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. ") + 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." @@ -621,8 +628,10 @@ 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. ") + 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 def fvac(self): @@ -769,8 +778,10 @@ 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. ") + raise ValueError( + "The pprime profile cannot be normalised. " + "Please first calculate Jtor for this profile. " + ) return self.L * self.Beta0 / self.Raxis * shape def ffprime(self, pn): @@ -791,8 +802,10 @@ 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. ") + 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 def fvac(self): @@ -892,22 +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,): + 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 = ( + 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 = ( + 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, @@ -998,7 +1016,10 @@ def Jtor_part2( 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 *= ( + 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 @@ -1006,14 +1027,17 @@ def Jtor_part2( 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 *= ( + 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) - + self.dJtordpsi = dJtordpsin / (psi_axis - psi_bndry) + # put to zero all current outside the LCFS # Jtor *= psi > psi_bndry @@ -1075,8 +1099,10 @@ def pprime(self, pn): shape = np.sum(shape, axis=0) 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. ") + 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): @@ -1107,8 +1133,10 @@ def ffprime(self, pn): shape = np.sum(shape, axis=0) 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. ") + 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):