Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Integrate on boundary to compute length scale quantities #1094

Draft
wants to merge 31 commits into
base: master
Choose a base branch
from

Conversation

unalmis
Copy link
Collaborator

@unalmis unalmis commented Jul 2, 2024

See #1101, in particular this comment for some of other changes.

Some quantities like the plasma volume V, major radius R0, minor radius a etc. can be computed via stoke's theorem on the last closed flux surface. This pull request modifies the compute functions to default to such an integration on the boundary rather than throughout the volume so that the length scale quantities can be computed accurately on grids other than QuadratureGrid. This lets us avoid building transforms for QuadratureGrid in objectives, saving memory.

Also, might have some improvements in general

  1. Reducing dimensionality of integral makes integration more tractable with faster convergence.
  2. Transforming the integration to one over periodic domains makes it more accurate due to the exponential convergence of the quadrature on periodic domains.

@unalmis unalmis added the speed New feature or request to make the code faster label Jul 2, 2024
@unalmis unalmis changed the base branch from master to fieldline_compute July 2, 2024 19:52
@unalmis unalmis marked this pull request as ready for review July 2, 2024 19:55
@unalmis unalmis marked this pull request as draft July 2, 2024 20:31
Copy link
Contributor

github-actions bot commented Jul 2, 2024

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_midres         |     +1.04 +/- 3.60     | +6.31e-03 +/- 2.17e-02 |  6.11e-01 +/- 1.3e-02  |  6.04e-01 +/- 1.7e-02  |
 test_build_transform_fft_highres        |     +0.68 +/- 2.57     | +6.73e-03 +/- 2.55e-02 |  1.00e+00 +/- 2.3e-02  |  9.93e-01 +/- 1.1e-02  |
 test_equilibrium_init_lowres            |     +0.66 +/- 1.35     | +2.50e-02 +/- 5.14e-02 |  3.84e+00 +/- 4.5e-02  |  3.81e+00 +/- 2.6e-02  |
 test_objective_compile_atf              |     -0.27 +/- 4.23     | -2.12e-02 +/- 3.29e-01 |  7.76e+00 +/- 2.4e-01  |  7.78e+00 +/- 2.3e-01  |
 test_objective_compute_atf              |     +0.40 +/- 1.97     | +4.20e-05 +/- 2.07e-04 |  1.05e-02 +/- 1.7e-04  |  1.05e-02 +/- 1.2e-04  |
 test_objective_jac_atf                  |     +1.05 +/- 3.54     | +1.97e-02 +/- 6.65e-02 |  1.90e+00 +/- 6.1e-02  |  1.88e+00 +/- 2.7e-02  |
 test_perturb_1                          |     -0.21 +/- 2.24     | -2.62e-02 +/- 2.82e-01 |  1.26e+01 +/- 1.7e-01  |  1.26e+01 +/- 2.3e-01  |
 test_proximal_jac_atf                   |     +0.60 +/- 1.00     | +4.84e-02 +/- 8.03e-02 |  8.12e+00 +/- 4.4e-02  |  8.07e+00 +/- 6.7e-02  |
 test_proximal_freeb_compute             |     -1.09 +/- 1.03     | -2.00e-03 +/- 1.90e-03 |  1.82e-01 +/- 1.3e-03  |  1.84e-01 +/- 1.4e-03  |
 test_build_transform_fft_lowres         |     +0.31 +/- 5.38     | +1.63e-03 +/- 2.86e-02 |  5.33e-01 +/- 2.3e-02  |  5.32e-01 +/- 1.7e-02  |
 test_equilibrium_init_medres            |     +0.87 +/- 1.71     | +3.65e-02 +/- 7.21e-02 |  4.24e+00 +/- 5.0e-02  |  4.20e+00 +/- 5.2e-02  |
 test_equilibrium_init_highres           |     +0.15 +/- 1.43     | +8.14e-03 +/- 7.96e-02 |  5.57e+00 +/- 6.7e-02  |  5.56e+00 +/- 4.3e-02  |
 test_objective_compile_dshape_current   |     -0.15 +/- 1.10     | -5.77e-03 +/- 4.29e-02 |  3.89e+00 +/- 3.9e-02  |  3.90e+00 +/- 1.8e-02  |
 test_objective_compute_dshape_current   |     -1.33 +/- 1.94     | -4.84e-05 +/- 7.08e-05 |  3.60e-03 +/- 5.6e-05  |  3.65e-03 +/- 4.4e-05  |
 test_objective_jac_dshape_current       |     +4.00 +/- 8.22     | +1.61e-03 +/- 3.30e-03 |  4.17e-02 +/- 3.1e-03  |  4.01e-02 +/- 1.2e-03  |
 test_perturb_2                          |     +0.06 +/- 1.20     | +1.08e-02 +/- 2.14e-01 |  1.79e+01 +/- 1.9e-01  |  1.79e+01 +/- 1.0e-01  |
 test_proximal_freeb_jac                 |     +0.09 +/- 1.13     | +6.41e-03 +/- 8.47e-02 |  7.52e+00 +/- 2.3e-02  |  7.51e+00 +/- 8.1e-02  |
 test_solve_fixed_iter                   |     -0.58 +/- 59.85    | -2.96e-02 +/- 3.06e+00 |  5.09e+00 +/- 2.1e+00  |  5.12e+00 +/- 2.2e+00  |

These differences are because, prior to the changes in this PR,
these length scale quantities would get computed on a full resolution
quadrature grid in test_compute_everything. Now thare computed on the same
low resolution used in the test that is used for everything else.

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 1.63995582e-05
Max relative difference: 2.03987171e-05
 x: array(0.803967)
 y: array(0.80395)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: S.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.16850603
Max relative difference: 0.00125191
 x: array(134.767312)
 y: array(134.598806)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A(z).
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 0.63797977
Max relative difference: 0.80369303
 x: array([0.763363, 0.763363, 0.763363, 0.763363, 0.763363, 0.763363,
       0.763363, 0.763363, 0.763363, 0.763363, 0.763363, 0.763363,
       0.763363, 0.763363, 0.763363, 0.763363, 0.763363, 0.763363,...
 y: array([0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,
       0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,
       0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: R0.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.00011245
Max relative difference: 2.0398301e-05
 x: array(5.51284)
 y: array(5.512953)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A(r).
Mismatched elements: 594 / 660 (90%)
Max absolute difference: 1.63995582e-05
Max relative difference: 2.03987171e-05
 x: array([0.      , 0.      , 0.      , 0.      , 0.      , 0.      ,
       0.010291, 0.010291, 0.010291, 0.010291, 0.010291, 0.010291,
       0.041077, 0.041077, 0.041077, 0.041077, 0.041077, 0.041077,...
 y: array([0.      , 0.      , 0.      , 0.      , 0.      , 0.      ,
       0.010291, 0.010291, 0.010291, 0.010291, 0.010291, 0.010291,
       0.041077, 0.041077, 0.041077, 0.041077, 0.041077, 0.041077,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: a.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 5.15953253e-06
Max relative difference: 1.01993065e-05
 x: array(0.505876)
 y: array(0.505871)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: R0/a.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.00033345
Max relative difference: 3.05972954e-05
 x: array(10.89761)
 y: array(10.897944)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: a_major/a_minor.
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 10.92082022
Max relative difference: 4.88235767
 x: array([ 4.007305,  4.007305,  4.007305,  4.007305,  4.007305,  4.007305,
        4.007305,  4.007305,  4.007305,  4.007305,  4.007305,  4.007305,
        4.007305,  4.007305,  4.007305,  4.007305,  4.007305,  4.007305,...
 y: array([4.009882, 4.009882, 4.009882, 4.009882, 4.009882, 4.009882,
       4.009882, 4.009882, 4.009882, 4.009882, 4.009882, 4.009882,
       4.009882, 4.009882, 4.009882, 4.009882, 4.009882, 4.009882,...
@unalmis unalmis linked an issue Jul 4, 2024 that may be closed by this pull request
@dpanici
Copy link
Collaborator

dpanici commented Jul 10, 2024

Add error to line_integrals if grid is sym and coord is theta

@dpanici
Copy link
Collaborator

dpanici commented Jul 10, 2024

change defualt grid in AspectRatio and Elongation to sym=False if thing=FourierRZToroidalSurface

Some of the master compute data changes with this commit.
The changes which are not due to floating point differences or simply grid
resolution differences (recall the grid used in test_compute_everything has
L = 9, M=N=5 which is less than required for convergence to true integral
quantities on that equilbrium) are given below:

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A(z).
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 0.6494644
Max relative difference: 0.81816076
 x: array([0.792591, 0.792591, 0.792591, 0.792591, 0.792591, 0.792591,
       0.792591, 0.792591, 0.792591, 0.792591, 0.792591, 0.792591,
       0.792591, 0.792591, 0.792591, 0.792591, 0.792591, 0.792591,...
 y: array([0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,
       0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,
       0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: a_major/a_minor.
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 11.33931768
Max relative difference: 4.99527722
 x: array([ 3.864152,  3.864152,  3.864152,  3.864152,  3.864152,  3.864152,
        3.864152,  3.864152,  3.864152,  3.864152,  3.864152,  3.864152,
        3.864152,  3.864152,  3.864152,  3.864152,  3.864152,  3.864152,...
 y: array([4.057173, 4.057173, 4.057173, 4.057173, 4.057173, 4.057173,
       4.057173, 4.057173, 4.057173, 4.057173, 4.057173, 4.057173,
       4.057173, 4.057173, 4.057173, 4.057173, 4.057173, 4.057173,...

NOTE that both quanties are incorrect, because in general we cannot compute
these quantities accuratly on grids that do not sample the full poloidal
domain.
@unalmis unalmis marked this pull request as ready for review July 11, 2024 22:23
@unalmis unalmis added documentation Add documentation or better warnings etc. bug fix Something was fixed and removed bug fix Something was fixed labels Jul 12, 2024
@unalmis unalmis linked an issue Jul 12, 2024 that may be closed by this pull request
@unalmis unalmis marked this pull request as draft July 12, 2024 22:19
@@ -75,6 +74,10 @@ def _V_FourierRZToroidalSurface(params, transforms, profiles, data, **kwargs):
profiles=[],
coordinates="r",
data=["e_theta", "e_zeta", "Z"],
parameterization=[
"desc.equilibrium.equilibrium.Equilibrium",
"desc.geometry.surface.FourierRZToroidalSurface",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using this for FourierRZToroidalSurface feels a bit weird since in that case its not really a funciton of rho.

)
# To approximate volume at ρ ~ 1, we scale by ρ⁻², assuming the integrand
# varies little from ρ = max_rho to ρ = 1 and a roughly circular cross-section.
data["V"] = jnp.max(data["V(r)"]) / jnp.max(data["rho"]) ** 2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we want to do this scaling here? IE, if the user creates a surface with rho=0.7, do we want to compute the actual volume? or the extrapolated volume?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a good point, if I want an interior volume this would still extrapolate. Maybe there should be a "V_full" key that does extrapolation and a "V" that just gives the volume? I guess that is just an alias for V(r) though

@@ -90,7 +90,7 @@ def test_1d_p(self):
def test_1d_elongation(self):
"""Test plotting 1d elongation as a function of toroidal angle."""
eq = get("precise_QA")
grid = LinearGrid(N=20, NFP=eq.NFP)
grid = LinearGrid(M=eq.M_grid, N=20, NFP=eq.NFP)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would prefer to add logic to eq.compute for that case so that it "just works", especially for plotting stuff. It especially feels weird if you need to create a 2d grid to make a 1d plot.

@@ -405,6 +405,7 @@ class TestZernikeRZToroidalSection:
"""Tests for ZernikeRZToroidalSection class."""

@pytest.mark.unit
@pytest.mark.xfail(reason="GitHub issue 1127.")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why is this failing now? Shouldn't a cross section be able to compute area?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was pushing this off to the PR that resolved issue 1127, but you are right that I should address this here.


@pytest.mark.unit
@pytest.mark.xfail(reason="GitHub issue 1127.")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see other comment, I think this should still work?

vmec.variables["Aminor_p"][:], desc.variables["Aminor_p"][:]
vmec.variables["Aminor_p"][:],
desc.variables["Aminor_p"][:],
rtol=2e-4,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't the new method of computing things be even more accurate?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the limit toward high resolution yes. Looks like at the resolution used in this test, the new method agrees less with vmec.

I could update my sympy script to compute the analytic result of these finite quadrature sampled on this grid for the SOLOVEV's boundary. This would prove whether the new method is less accurate at this resolution for these sample points.

I haven't done that to confirm. However, I had done similar for other surfaces e.g. here, and the numerical quadrature agreed exactly with the analytic quadrature (both of which converge to the true integral).

resolution_requirement="rtz",
)
def _V(params, transforms, profiles, data, **kwargs):
data["V"] = jnp.sum(data["sqrt(g)"] * transforms["grid"].weights)
if isinstance(transforms["grid"], QuadratureGrid):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can we include logic like this for some of the other cases? My main worry is that for higher order quantities that depend on V,A,a etc we're making it so that you need multiple different grids just to compute everything correctly, but i think in theory a quadrature grid should work for everything (though may be overkill in some cases)

@@ -238,6 +240,8 @@ def __init__(
):
if target is None and bounds is None:
target = 1
if grid is not None:
errorif_sym(grid, name)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we need to check this for other objectives?

Copy link
Collaborator Author

@unalmis unalmis Jul 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

aspect ratio can use symmetry. I checked other objectives were also fine. (haven't looked at new coil objectives)

Copy link
Collaborator

@dpanici dpanici left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few questions in the comments, main one is, will "a_major/a_minor" be computed on the correct grid by default when called in eq.compute? how about when passed into plot_1d? I think it should if not

@@ -39,27 +46,19 @@ def _V(params, transforms, profiles, data, **kwargs):
label="V",
units="m^{3}",
units_long="cubic meters",
description="Volume",
description="Volume enclosed by surface, scaled by max(ρ)⁻²",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

idk if this docstring is clear enough, something like "if grid does not contain rho=1" or something could possibly help?

@@ -155,55 +158,32 @@ def _V_rrr_of_r(params, transforms, profiles, data, **kwargs):
label="A(\\zeta)",
units="m^{2}",
units_long="square meters",
description="Cross-sectional area as function of zeta",
description="Area of enclosed cross-section (enclosed constant phi surface), "
"scaled by max(ρ)⁻², as function of zeta",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same point on clarity as above

desc/compute/_geometry.py Show resolved Hide resolved
@@ -275,25 +282,22 @@ def _A_of_r(params, transforms, profiles, data, **kwargs):
label="S",
units="m^{2}",
units_long="square meters",
description="Surface area of outermost flux surface",
description="Surface area of outermost flux surface, scaled by max(ρ)⁻¹",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you referring to the docstring saying scaled by max rho? I think so

grid.sym,
ResolutionWarning,
msg=colored(
"This grid only samples the poloidal domain θ ∈ [0, π], "
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should mention stellarator symmetry in the warning as well (this grid has stellarator sym, which means it only samples...)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe also say that a grid with sym=False should be used instead

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the first note, in a related issue mentioned at top, I mention the goal

Update documentation of Grid class to reflect that grid.sym is unrelated to stellarator symmetry, rather it's a memory optimization to compute certain quantities under symmetry.

So I don't want to refer to grid as having stellarator symmetry.

errorif(
grid.sym,
ResolutionWarning,
msg=colored(
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ditto above comment

desc/objectives/_geometry.py Outdated Show resolved Hide resolved
@@ -90,7 +90,7 @@ def test_1d_p(self):
def test_1d_elongation(self):
"""Test plotting 1d elongation as a function of toroidal angle."""
eq = get("precise_QA")
grid = LinearGrid(N=20, NFP=eq.NFP)
grid = LinearGrid(M=eq.M_grid, N=20, NFP=eq.NFP)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this PR now, will plot_1d("a_major/a_minor") work by default (i.e. will the default grid be sufficient)? If this is a quantity that needs the poloidal resolution, I think the .compute logic should by default compute it on the right grids

@unalmis unalmis marked this pull request as draft July 24, 2024 22:04
@unalmis unalmis added add warning Adds warning to prevent/detect bug and removed documentation Add documentation or better warnings etc. labels Jul 25, 2024
@unalmis unalmis removed the EZ-review label Aug 10, 2024
@unalmis unalmis self-assigned this Aug 17, 2024
Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.00127829
Max relative difference: 0.00159001
 x: array(0.805229)
 y: array(0.80395)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A(z).
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 0.6494644
Max relative difference: 0.81816076
 x: array([0.792591, 0.792591, 0.792591, 0.792591, 0.792591, 0.792591,
       0.792591, 0.792591, 0.792591, 0.792591, 0.792591, 0.792591,
       0.792591, 0.792591, 0.792591, 0.792591, 0.792591, 0.792591,...
 y: array([0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,
       0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,
       0.762978, 0.762978, 0.762978, 0.762978, 0.762978, 0.762978,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: R0.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.00875172
Max relative difference: 0.00158748
 x: array(5.504201)
 y: array(5.512953)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: A(r).
Mismatched elements: 594 / 660 (90%)
Max absolute difference: 0.00127829
Max relative difference: 0.00159001
 x: array([0.      , 0.      , 0.      , 0.      , 0.      , 0.      ,
       0.010307, 0.010307, 0.010307, 0.010307, 0.010307, 0.010307,
       0.041142, 0.041142, 0.041142, 0.041142, 0.041142, 0.041142,...
 y: array([0.      , 0.      , 0.      , 0.      , 0.      , 0.      ,
       0.010291, 0.010291, 0.010291, 0.010291, 0.010291, 0.010291,
       0.041077, 0.041077, 0.041077, 0.041077, 0.041077, 0.041077,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: a.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.00040201
Max relative difference: 0.00079469
 x: array(0.506273)
 y: array(0.505871)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: perimeter(z).
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 0.05135033
Max relative difference: 0.01242666
 x: array([4.250749, 4.250749, 4.250749, 4.250749, 4.250749, 4.250749,
       4.250749, 4.250749, 4.250749, 4.250749, 4.250749, 4.250749,
       4.250749, 4.250749, 4.250749, 4.250749, 4.250749, 4.250749,...
 y: array([4.231189, 4.231189, 4.231189, 4.231189, 4.231189, 4.231189,
       4.231189, 4.231189, 4.231189, 4.231189, 4.231189, 4.231189,
       4.231189, 4.231189, 4.231189, 4.231189, 4.231189, 4.231189,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: a_major/a_minor.
Mismatched elements: 660 / 660 (100%)
Max absolute difference: 11.37253303
Max relative difference: 5.08430436
 x: array([ 3.864152,  3.864152,  3.864152,  3.864152,  3.864152,  3.864152,
        3.864152,  3.864152,  3.864152,  3.864152,  3.864152,  3.864152,
        3.864152,  3.864152,  3.864152,  3.864152,  3.864152,  3.864152,...
 y: array([4.009882, 4.009882, 4.009882, 4.009882, 4.009882, 4.009882,
       4.009882, 4.009882, 4.009882, 4.009882, 4.009882, 4.009882,
       4.009882, 4.009882, 4.009882, 4.009882, 4.009882, 4.009882,...

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: S.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.17114453
Max relative difference: 0.00127152
 x: array(134.76995)
 y: array(134.598806)

Not equal to tolerance rtol=1e-10, atol=1e-10
Parameterization: desc.equilibrium.equilibrium.Equilibrium. Name: R0/a.
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 0.02594014
Max relative difference: 0.00238028
 x: array(10.872004)
 y: array(10.897944)
@unalmis unalmis mentioned this pull request Sep 26, 2024
@unalmis
Copy link
Collaborator Author

unalmis commented Sep 27, 2024

would be less annoying to resolve #1127 before this

# where n is the unit normal such that n dot e_θ|ρ,ϕ = 0 and n dot e_ϕ|R,Z = 0,
# and the labels following | denote those coordinates are fixed.
# Now choose v = [0, 0, Z], and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1].
n = data["n_rho"]
Copy link
Collaborator Author

@unalmis unalmis Oct 2, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

iirc because n_rho should not be defined on zeta cross section (#1127 ) the task is to rewrite this to obtain normal vector purely from e_theta which is defined on zeta cross section alone. basically like #597 (comment).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
add warning Adds warning to prevent/detect bug speed New feature or request to make the code faster
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Change description of S
4 participants