diff --git a/freegs4e/critical.py b/freegs4e/critical.py index c840617..0228d92 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,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 diff --git a/freegs4e/equilibrium.py b/freegs4e/equilibrium.py index 1a3fbf7..c39aa5a 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): """ @@ -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( diff --git a/freegs4e/jtor.py b/freegs4e/jtor.py index b332c36..e2837a8 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,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):