From f422e8d060d04cbe8275e19df9888339e8aaa120 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 13 Nov 2023 10:49:35 -0500 Subject: [PATCH 1/6] Use numpy in optimise --- freegs/optimise.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/freegs/optimise.py b/freegs/optimise.py index f2e842b..cbe1942 100644 --- a/freegs/optimise.py +++ b/freegs/optimise.py @@ -27,7 +27,7 @@ import matplotlib.pyplot as plt from freegs.plotting import plotEquilibrium -from math import sqrt +from numpy import sqrt, inf # Measures which operate on Equilibrium objects @@ -199,7 +199,7 @@ def solve_and_measure(eq): return measure(eq) except: # Solve failed. - return float("inf") + return float(inf) # Call the generic optimiser, return optimiser.optimise( From f56e34470fda29dc4ca777f258007d264271bb12 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 13 Nov 2023 10:51:13 -0500 Subject: [PATCH 2/6] Escape special chars in docstring in gradshafranov --- freegs/gradshafranov.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/freegs/gradshafranov.py b/freegs/gradshafranov.py index 2b63a56..f8c8c04 100644 --- a/freegs/gradshafranov.py +++ b/freegs/gradshafranov.py @@ -38,7 +38,7 @@ class GSElliptic: Represents the Grad-Shafranov elliptic operator .. math:: - \Delta^* = R^2 \nabla\cdot\frac{1}{R^2}\nabla + \\Delta^* = R^2 \\nabla\\cdot\\frac{1}{R^2}\\nabla which is From 39dad8e1f022230d8aed3d50397d2f6f9d132af8 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 13 Nov 2023 10:51:34 -0500 Subject: [PATCH 3/6] Use numpy in optimiser --- freegs/optimiser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/freegs/optimiser.py b/freegs/optimiser.py index 2c6ea49..1701559 100644 --- a/freegs/optimiser.py +++ b/freegs/optimiser.py @@ -24,7 +24,7 @@ along with FreeGS. If not, see . """ -import random +from numpy import random import copy import bisect import sys From 6c7c253e2d018733e7214bc213132253e37797e5 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 13 Nov 2023 10:51:48 -0500 Subject: [PATCH 4/6] Use inf as initial value in picard --- freegs/picard.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/freegs/picard.py b/freegs/picard.py index 81474d1..94d40ee 100644 --- a/freegs/picard.py +++ b/freegs/picard.py @@ -93,7 +93,7 @@ def solve( bndry = 0.0 # Set an initial value for bndry_change (set high to prevent immediate convergence) - bndry_change = 10.0 + bndry_change = np.inf # Plasma assumed to not be limited at first has_been_limited = False From 24189271da37fbce5af35f9f92e6b3148a8c8a8f Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 13 Nov 2023 10:54:17 -0500 Subject: [PATCH 5/6] Ensure q on axis is smooth. --- freegs/equilibrium.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/freegs/equilibrium.py b/freegs/equilibrium.py index 79f3508..611f93d 100644 --- a/freegs/equilibrium.py +++ b/freegs/equilibrium.py @@ -338,18 +338,30 @@ def q(self, psinorm=None, npsi=100): >>> psinorm, q = eq.q() Note: psinorm = 0 is the magnetic axis, and psinorm = 1 is the separatrix. - Calculating q on either of these flux surfaces is problematic, - and the results will probably not be accurate. + If you request a value close to either of these limits, an extrapolation + based on a 1D grid of values from 0.01 to 0.99 will be used. This gives + smooth and continuous q-profiles, but comes at an increased computational + cost. """ if psinorm is None: # An array which doesn't include psinorm = 0 or 1 psinorm = linspace(1.0 / (npsi + 1), 1.0, npsi, endpoint=False) return psinorm, critical.find_safety(self, psinorm=psinorm) - result = critical.find_safety(self, psinorm=psinorm) + elif np.any((psinorm < 0.01) | (psinorm > 0.99)): + psinorm_inner = np.linspace(0.01, 0.99, num=npsi) + q_inner = critical.find_safety(self, psinorm=psinorm_inner) + + interp = interpolate.interp1d(psinorm_inner, q_inner, + kind = "quadratic", + fill_value = "extrapolate") + result = interp(psinorm) + else: + result = critical.find_safety(self, psinorm=psinorm) + # Convert to a scalar if only one result - if len(result) == 1: - return result[0] + if np.size(result) == 1: + return float(result) return result def tor_flux(self, psi=None): From d2342b1c434b3c3150ed9a44310cf90de1ada164 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 13 Nov 2023 11:07:43 -0500 Subject: [PATCH 6/6] Fix select of a single random element --- freegs/optimiser.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/freegs/optimiser.py b/freegs/optimiser.py index 1701559..4b4626c 100644 --- a/freegs/optimiser.py +++ b/freegs/optimiser.py @@ -65,7 +65,9 @@ def pickUnique(N, m, e): inds = sorted(e) # Sorted list of indices. Used to avoid clashes others = [] # The list of three agents for i in range(m): - newind = random.randint(0, N - 1 - i - len(e)) + high = N - 1 - i - len(e) + newind = random.randint(0, high) if high > 0 else 0 + for ind in inds: if newind == ind: newind += 1