Why does my EQDSK need lower last_surface_factor?

[theo-brown]

I'm parsing STEP eqdsk files from JETTO to run comparisons in TORAX, and I'm finding that I need to reduce the CONFIG["geometry"]["last_surface_factor"] from its default value of 0.99 to 0.95.
Without this I get strongly reversed q profiles at the edge.

The grid that the ψ profile from the eqdsk is interpolated onto is:

torax/torax/geometry/standard_geometry.py

Lines 686 to 690 in 3ea75f5

	psi_interpolant = np.linspace(
	0,
	(eqfile['psibdry'] - eqfile['psimag']) * last_surface_factor,
	n_surfaces,
	)

So this looks like decreasing last_surface_factor moves the effective LCFS further inside the true LCFS.

I'd like to understand a bit more the reasons behind why this needs to be done in my case. I think it might be because STEP has a thin layer of extremely high current at its edge.

For reference, the eqdsk file that I'm reading looks something like this, from Mitchell et al 2023:

Tagging @fusionby2030 as he is an EQDSK wizard.

[jcitrin]

Is it a diverted geometry? Are quantities like surface_FSA_Bpol_squared going down towards zero near the edge?

the q_LCFS isn't just a function of current, it also depends on the edge psi_grad which goes like Ip_total / g2g3_over_rhon_face[-1], i.e. the q-profile diverges at the LCFS in diverted geometry. So if something weird is happening with q_LCFS, worth seeing how g2g3_over_rhon_face is calculated.

This g2g3_over_rhon_face is formally zero in diverted geometry because part of the flux-surface-average ends up in the infinite toroidal loop where Bpol=0 . I don't recall exactly how that is reflected in EQDSK geometry, but if there is a surface_Bpol~0 somewhere, then that's "bad news" :) , and 0.99 may be too close to that.

But having said all that, in your case q is dropping, not diverging. Maybe it's an issue with

  surface_FSA_Bpol_squared = (
          np.sum(surface_Bpol * surface_dl) / surface_int_dl_over_bpol
      )

?

In fact, is that a correct formula? I'd naively expect surface_Bpol**2 there. @fusionby2030 ?

Edit, the formula is fine (the Bpol in the metric dl/Bpol is cancelled out by the Bpol**2). However, it's still a place worth looking at.

About a thin current layer, if that's really the case (how does that arise?) then it should be seen also in Ip_eqdsk

[theo-brown]

I've plotted the surface_dl and surface_Bpol - I'm not even sure whether plotting them this way makes sense, but it was just to sense-check things vaguely.

Each subplot has the values for a particular surface highlighted in red.

It definitely looks like something funny is happening on the last surface, because the traces are half as long as they are on other surfaces. What's meant to happen here? Is there a better way to visualise it? Again, sorry for potentially nonsensical questions or plots, I'm not familiar with the method or science here and I've struggled to wrap my head around it.

My vague guess would be that this half-length array somehow represents the FS not being fully closed?

[fusionby2030]

Can you plot them as function of theta? Should range from 0:2pi

theta= np.arctan2(Z, R)

Or theta = arctan(Z-ZMID, R-RMID) where Z and R are the flux surface coordinates.

[theo-brown]

Ah yes, that makes more sense as a coordinate!

Just to check, in the code R = x_surface and Z = z_surface?

torax/torax/geometry/standard_geometry.py

Lines 716 to 725 in 0da2ed6

	for _, _psi in enumerate(psi_interpolant[1:]):
	vertices = cg_psi.create_contour(_psi)
	if not vertices:
	raise ValueError(f"""
	Valid contour not found for EQDSK geometry for psi value {_psi}.
	Possible reason is too many surfaces requested.
	Try reducing n_surfaces from the current value of {n_surfaces}.
	""")
	x_surface, z_surface = vertices[0].T[0], vertices[0].T[1]
	surfaces.append((x_surface, z_surface))

[theo-brown]

Plotting just the surfaces (output of the above loop) provides an insight:

While the LCFS from the EQDSK is closed (red), the last contour-generated surface (green) doesn't manage to close properly and is missing its lower half. Will this pose a problem for the rest of the routine? I'm not sure whether this is the source of the problem or a red herring.

I've added a PR that checks for this case, #895.

[jcitrin]

Thanks for looking into this. I'll look at the PR now. Sorry for delays. Great to do a check; I'm also wondering if it's possible to fix the routine, but let's stick with the check for now.