Skip to content

Commit

Permalink
Merge pull request #138 from boutproject/better-initial-nonorth-separ…
Browse files Browse the repository at this point in the history 
…atrix-spacing

Better orthogonal separatrix spacing for nonorthogonal grids
  • Loading branch information
@johnomotani
johnomotani authored Nov 28, 2022
2 parents 4ec834f + 1376977 commit fc72bd7
Show file tree
Hide file tree
Showing 3 changed files with 156 additions and 120 deletions.
12 changes: 8 additions & 4 deletions .flake8
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,14 @@
[flake8]
max-line-length = 89
#E741 - 'ambiguous variable names' forbids using 'I', 'O' or 'l'
#W503 - 'line break before binary operator', but this is allowed and useful inside brackets
#E203 - 'whitespace before ':'', but black formats some slice expressions with space before ':'
#E231 - missing whitespace after ',', but black formats some expressions without space after ','
ignore =
E741, # 'ambiguous variable names' forbids using 'I', 'O' or 'l'
W503, # 'line break before binary operator', but this is allowed and useful inside brackets
E203, # 'whitespace before ':'', but black formats some slice expressions with space before ':'
E231, # missing whitespace after ',', but black formats some expressions without space after ','
E741,
W503,
E203,
E231,
exclude =
hypnotoad/gui/hypnotoad_mainWindow.py,
hypnotoad/gui/hypnotoad_preferences.py,
Expand Down
22 changes: 21 additions & 1 deletion doc/whats-new.md
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,29 @@
What's new
==========

0.4.5 (unreleased)
0.5.0 (unreleased)
------------------

### Breaking changes

- Change how initial spacing of points on separatrix for nonorthogonal grids is
calculated (this spacing is used to construct the underlying orthogonal
grid). Now use weights that are just
{`cos(0.5*pi*i/index_length)**2`,`sin(0.5*pi*i/index_length)**2`} rather than
trying to construct a function that will give the same spacings when called
with the 'orthogonal spacing' that came from the initial call passed in as
`sfunc_orthogonal`. This means that the underlying orthogonal grid does not
rely at all on `nonorthogonal_*` settings, so removes an issue where some
grids could only be constructed by generating first with one set of
`nonorthogonal_*` settings, then changing them and regridding. It also
removes the possibility of a "Weight too small. Suggest increasing poloidal
'range' settings" error, since the branch that produced that error no longer
exists. Does change the nonorthogonal grids that will be produced by
hypnotoad, although the integrated tests pass without updating the expected
results, so the changes should be small (less than 1e-9 for the integrated
test cases) (#138)\
By [John Omotani](https://github.com/johnomotani)

### Bug fixes

- Catch error when `psi_sol_inner` is set wrong. Add suggestions to 'gradient
Expand Down
242 changes: 127 additions & 115 deletions hypnotoad/core/equilibrium.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2884,142 +2884,154 @@ def combineSfuncs(
and spacings["nonorthogonal_range_upper"] is not None
):

def new_sfunc(i):
sfixed_lower = sfunc_fixed_lower(i)
if sfunc_orthogonal is None:
# Define new_sfunc in a sensible way to create the initial distribution
# of points on the separatrix that is then used to create the orthogonal
# grid. This shouldn't rely on 'range' parameters as it should be fixed
# when nonorthogonal_* parameters are changed. Also want to avoid a
# sharp transition between sfixed_lower and sfixed_upper which might
# give odd spacings
def new_sfunc(i):
sfixed_lower = sfunc_fixed_lower(i)
sfixed_upper = sfunc_fixed_upper(i)

# define weight_lower so it is 1. at the lower boundary and 0. at the
# upper boundary and the gradient is zero at both boundaries
weight_lower = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
1.0,
0.0,
lambda i: numpy.cos(0.5 * numpy.pi * i / index_length) ** 2,
],
)

sfixed_upper = sfunc_fixed_upper(i)
# define weight_upper so it is 1. at the upper boundary and 0. at the
# lower boundary and the gradient is zero at both boundaries
weight_upper = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
0.0,
1.0,
lambda i: numpy.sin(0.5 * numpy.pi * i / index_length) ** 2,
],
)

if sfunc_orthogonal is None:
sorth = None
else:
# Note weight_lower+weight_upper=1 by construction
return weight_lower * sfixed_lower + weight_upper * sfixed_upper

else:

def new_sfunc(i):
sfixed_lower = sfunc_fixed_lower(i)
sfixed_upper = sfunc_fixed_upper(i)
sorth = sfunc_orthogonal(i)

# define weight_lower so it is 1. at the lower boundary and 0. at the
# upper boundary and the gradient is zero at both boundaries
weight_lower = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
1.0,
0.0,
lambda i: numpy.exp(-((i / N_norm / this_range_lower) ** 2)),
],
)
# define weight_lower so it is 1. at the lower boundary and 0. at the
# upper boundary and the gradient is zero at both boundaries
weight_lower = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
1.0,
0.0,
lambda i: numpy.exp(
-((i / N_norm / this_range_lower) ** 2)
),
],
)

# define weight_upper so it is 1. at the upper boundary and 0. at the
# lower boundary and the gradient is zero at both boundaries
weight_upper = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
0.0,
1.0,
lambda i: numpy.exp(
-(((index_length - i) / N_norm / this_range_upper) ** 2)
),
],
)
# define weight_upper so it is 1. at the upper boundary and 0. at the
# lower boundary and the gradient is zero at both boundaries
weight_upper = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
0.0,
1.0,
lambda i: numpy.exp(
-(((index_length - i) / N_norm / this_range_upper) ** 2)
),
],
)

# make sure weight_lower + weight_upper <= 1
weight = weight_lower + weight_upper
weight_over_slice = weight[weight > 1.0]
weight_lower[weight > 1.0] /= weight_over_slice
weight_upper[weight > 1.0] /= weight_over_slice

if sorth is None:
# Fix spacing so that if we call combineSfuncs again for this contour
# with sfunc_orthogonal from self.contourSfunc() then we get the same
# spacing again. This is used to make the contours along the
# separatrix keep the same values when pushing the other contours
# towards orthogonal positions
# s = weight_lower*sfixed_lower
# + weight_upper*sfixed_upper
# + (1. - weight_lower - weight_upper)*s
sorth = (
weight_lower * sfixed_lower + weight_upper * sfixed_upper
) / (weight_lower + weight_upper)

if numpy.any(weight_lower + weight_upper < 1e-200):
print("radial index", ix)
print(weight_lower + weight_upper)
raise ValueError(
"Weight too small. Suggest increasing poloidal 'range' "
"settings"
)
# make sure weight_lower + weight_upper <= 1
weight = weight_lower + weight_upper
weight_over_slice = weight[weight > 1.0]
weight_lower[weight > 1.0] /= weight_over_slice
weight_upper[weight > 1.0] /= weight_over_slice

return (
weight_lower * sfixed_lower
+ weight_upper * sfixed_upper
+ (1.0 - weight_lower - weight_upper) * sorth
)
return (
weight_lower * sfixed_lower
+ weight_upper * sfixed_upper
+ (1.0 - weight_lower - weight_upper) * sorth
)

elif spacings["nonorthogonal_range_lower"] is not None:

def new_sfunc(i):
sfixed_lower = sfunc_fixed_lower(i)
if sfunc_orthogonal is None:
# Fix spacing so that if we call combineSfuncs again for this contour
# with sfunc_orthogonal from self.contourSfunc() then we get the same
# spacing again. This is used to make the contours along the separatrix
# keep the same values when pushing the other contours towards
# orthogonal positions
# s = sorth = weight_lower*sfixed_lower + (1. - weight_lower)*sorth
new_sfunc = sfunc_fixed_lower
else:

if sfunc_orthogonal is None:
sorth = None
else:
def new_sfunc(i):
sfixed_lower = sfunc_fixed_lower(i)
sorth = sfunc_orthogonal(i)

# define weight_lower so it is 1. at the lower boundary and the gradient
# is zero at the lower boundary.
weight_lower = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
1.0,
0.0,
lambda i: numpy.exp(-((i / N_norm / this_range_lower) ** 2)),
],
)

if sorth is None:
# Fix spacing so that if we call combineSfuncs again for this contour
# with sfunc_orthogonal from self.contourSfunc() then we get the same
# spacing again. This is used to make the contours along the
# separatrix keep the same values when pushing the other contours
# towards orthogonal positions
# s = weight_lower*sfixed_lower + (1. - weight_lower)*s
sorth = sfixed_lower
# define weight_lower so it is 1. at the lower boundary and the
# gradient is zero at the lower boundary.
weight_lower = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
1.0,
0.0,
lambda i: numpy.exp(
-((i / N_norm / this_range_lower) ** 2)
),
],
)

return (weight_lower) * sfixed_lower + (1.0 - weight_lower) * sorth
return (weight_lower) * sfixed_lower + (1.0 - weight_lower) * sorth

elif spacings["nonorthogonal_range_upper"] is not None:

def new_sfunc(i):
sfixed_upper = sfunc_fixed_upper(i)
if sfunc_orthogonal is None:
# Fix spacing so that if we call combineSfuncs again for this contour
# with sfunc_orthogonal from self.contourSfunc() then we get the same
# spacing again. This is used to make the contours along the separatrix
# keep the same values when pushing the other contours towards
# orthogonal positions
# s = sorth = weight_upper*sfixed_upper + (1. - weight_upper)*sorth
new_sfunc = sfunc_fixed_upper
else:

if sfunc_orthogonal is None:
sorth = None
else:
def new_sfunc(i):
sfixed_upper = sfunc_fixed_upper(i)
sorth = sfunc_orthogonal(i)

# define weight_upper so it is 1. at the upper boundary and the gradient
# is zero at the upper boundary.
weight_upper = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
0.0,
1.0,
lambda i: numpy.exp(
-(((index_length - i) / N_norm / this_range_upper) ** 2)
),
],
)

if sorth is None:
# Fix spacing so that if we call combineSfuncs again for this contour
# with sfunc_orthogonal from self.contourSfunc() then we get the same
# spacing again. This is used to make the contours along the
# separatrix keep the same values when pushing the other contours
# towards orthogonal positions
# s = weight_upper*sfixed_upper + (1. - weight_upper)*s
sorth = sfixed_upper
# define weight_upper so it is 1. at the upper boundary and the
# gradient is zero at the upper boundary.
weight_upper = numpy.piecewise(
i,
[i < 0.0, i > index_length],
[
0.0,
1.0,
lambda i: numpy.exp(
-(((index_length - i) / N_norm / this_range_upper) ** 2)
),
],
)

return (weight_upper) * sfixed_upper + (1.0 - weight_upper) * sorth
return (weight_upper) * sfixed_upper + (1.0 - weight_upper) * sorth

else:
if sfunc_orthogonal is None:
Expand Down

0 comments on commit fc72bd7

Please sign in to comment.