Errors in core and edge profile interpolation and extrapolation

[anchal-physics]

I ran a scan of core and edge profiles of electron density on our sample file to see how the interpolation functions are doing. While for some values of Z and R, they are working exactly as we would like them to, there are weird problems with them on some values of R, Z. See these plots:

---

Scan Z at different values of R. The X-point is near R=5.5, Z=-3.5

1. At R=4.5m, 5.0m, and 5.5m, we expect the core profile to not have any values outside of core. But for some reason, it is returning this constant value at Z< -4m which is very similar to core value.
2. In almost all the plots, there is some kind of peak in edge_profile below Z=-4m. This is actually because of the points' proximity to separatix. So even though this region is outside the SOLPS mesh, it is smoothly decreasing the electron density away from separatix. Is this ok?
3. There is a sudden drop in core profile value near the center at R=6.5. Seems like a hole.

---

Scan R at different values of Z. The X-point is near R=5.5, Z=-3.5

1. At Z=-4.0m, -4.5m, again we expect the core_profile to not have any zero value as the core does not exist in this region. But the core_profile_2d returned non-zero values in some region.
2. At Z = -4.0, the behavior of edge_profile shows that the simple inverse distance weighing for interpolation has weird effects when closest known value can suddenly change like it happened a near R=5.5, we see a sudden drop. Maybe we should just not use this function outside the SOLPS boundary or make it output flat zero everywhere outside or throw an error if outside point is provided.

3. Same as above, one hole can be seen at Z=0.5m scan near the center. The core profile suddenly drops to zero (which is not plotted in this logscale plot)

[eldond]

Scan Z at different values of R. The X-point is near R=5.5, Z=-3.5, # 1

The double valued core profile thing is because of the interpolation through psi. psi_N has similar values in the core plasma and private flux region. So any function that only considers psi_N will not distinguish between PFR and core. This happens all the time and is pretty easy to get around by checking against X-point height.

By the way, we need X-point height. I'm working on it (https://github.com/ProjectTorreyPines/SD4SOLPS.jl/issues/22).

Scan Z at different values of R. The X-point is near R=5.5, Z=-3.5, # 2

I don't know how edge profiles is being extrapolated that far out beyond the SOLPS mesh. I didn't finish that thing, right? What does the green highlight mean?

The core profiles assume all quantities are constant on flux surfaces (the flux function approximatin), and the edge profiles do not. The flux function approximation is usually good deep in the core, and is useless in the SOL. Somewhere in between, the approximation can be seen breaking down. I based the core profile extrapolation on the outboard midplane, so I would expect core and edge profiles to match exactly at the outboard midplane and have errors elsewhere. See SOLPS-only ne(R,Z) (or Te) 2D plots to help clarify this.

Okay, but as for the specific variations that are seen here... some look weird. The large spike in the slice at Z=-2.5 m, for instance, is unexpected. That sort of thing is sometimes seen on the high field (smaller R) side. It's called the high-field-side-high-density. It's a very creative name.

[anchal-physics]

• The edge_profile interpolation/extrapolation works on this simple algorithm:

  • Find four closest grid points
  • Take their mean with weights as inverse of distance to those points.

  Thus this edge_profile extrapolation gives values even outside the SOLPS mesh, gradually decreasing the value as 1/d if d is the distance from SOL edge.

• The green region is the SOLPS mesh region. The red region is the core region (that is anything inside the separatix). These regions have a small overlap as SOLPS mesh has some core region as well.

• Just FYI, in the SOLPS mesh, X-point is defined and grid_subset with index 6 holds the X-point value. But you probably are looking for a general tool that uses magnetic field to find those in https://github.com/ProjectTorreyPines/SD4SOLPS.jl/issues/22

• Yeah the peak at around R=6.5, Z=-2.5 is bit weird. It is the separatix region. We see this peak in electron density just outside of separatix in teh heatmap as well.

• What did you think about the holes in the core profile near the center. See at R=6.5, Z = 0.5.

[eldond]

Dude that is really not bad for a simple first try at edge extrapolation. This significantly lowers the urgency of https://github.com/ProjectTorreyPines/SD4SOLPS.jl/issues/4 or maybe even solves it, depending on how far we want to take that thing. Anyway, it gives time to do a higher quality edge extrapolation while we use what you have for testing. Good work. Very nice.

Can we have a look at your extrapolated edge profile in 2d? I'm curious now.

The gap in the core should not be there. Could there have been some rounding error or use of the wrong value for psi axis (such as from another time or something?) such that psi_N < 0 deep in the core?

[eldond]

And yes, I need a more general X-point finder that operates on equilibrium data so that I can find a secondary X-point that's outside of the SOLPS grid.

[anchal-physics]

"Can we have a look at your extrapolated edge profile in 2d?"

Here is the plot made by only using the interp function from GGDUtils.jl at efd8468 To reiterate, this interpolation/extrapolation function simply looks for the closest 4 (which is default, can be made any number at cost of run time) grid points in SOLPS grid where the value is defined and takes their weighted mean weighted by inverse distance.

Minimum code to generate this:

using SOLPS2IMAS
using GGDUtils
b2gmtry = "GGDUtils.jl/samples/b2fgmtry"
b2output = "GGDUtils.jl/samples/b2time_red.nc"
gsdesc = "GGDUtils.jl/samples/gridspacedesc.yml"
b2mn = "GGDUtils.jl/samples/b2mn.dat"
ids = solps2imas(b2gmtry, b2output, gsdesc, b2mn)
grid_ggd = ids.edge_profiles.grid_ggd[1]
n_e  =  GGDUtils.interp(ids.edge_profiles.ggd[1].electrons.density, grid_ggd, 5)
R = 3:.01:9
Z=  -5:0.01:5
rz = [(r, z) for r in R, z in Z]
n_e_rz = n_e.(rz)
heatmap(R, Z, n_e_rz', c=:inferno, xlabel="R / m",
ylabel="Z / m", title=L"n_e / m^{-3}", colorbar_scale=:log10)
p = plot!(space)

Inferences

• This interpolation is doing an ok job of interpolating inside the core as well.
• Outside the SOL, the behavior is what you would expect for a 1/d decaying profile to look like where d is the distance from SOL.
• Near the targets, the behavior is clearly weird and probably wrong.
• Since the grid i highly irregular, maybe it makes more sense to use more than 4 points

Using more nearest neighbors

As I was writing this, I figures I can quickly try using more nearest points. Here are the plots:

Using 8 nearest neighbors

Using 16 nearest neighbors

Using 32 nearest neighbors

• So the halo gains more prominent features on increasing the use of nearest neighbors.
• The distribution inside separatix is essentially the same for all cases.
• There is actually no time penalty on using more points. Seems like kdtree is very efficient in finding any number of nearest neighbors for our cases.
• But it seems it is best to use just 4 points. The behavior near divertor targets is weird because of the choice of this extrapolation curve. If you have something else in mind, I can quickly try it out or add an argument in interp function to provide what extrapolation curve to use.

[eldond]

Well, some of these look like they could be a glorious and heroic logo or insignia, so that's nice. Very pretty. But these things are producing artifacts even within the mesh.

That's no good.

1/r doesn't fall off fast enough.

The behavior at the targets is just nuts.

The fall off length scale will vary poloidally and needs information about the flux surfaces.

To try to improve this for now, since it does seem pretty useful for testing, the behavior within the mesh has to not do this jagged pattern, and the rays of glory from the targets have to be suppressed. Outside the mesh, you can probably get some improvement by changing to exp(-r/lambda). Ideally, one would define lambda from the last few rows of the mesh, but you can also just pick something in the range 1 mm to 1 cm or so and play with it until it looks reasonable-ish (I think we're just doing test samples here).

In order for values within the mesh to be lower than what I think must be the values at the mesh centers, I suspect there's a flaw in how the weights are normalized, but I'm just speculating.

[anchal-physics]

It suffices to conclude then that outside SOLPS mesh, the extrapolation needs to be more careful and follow a different approach. I realized one error in my understanding of what is happening outside the mesh. Since the algorithm takes weighted mean of inverse distance, it actually just picks the closest cell value when it is outside the mesh. That's why even if I use an exponentially decaying weighing function with decay length of just 1 cm, the extrapolated profile looks exactly the same as 1/d weighing function.

Now that I think about it, I remember knowing this at the time when I wrote this function. That outside the grid, it will just return the value of the closest grid cell. That's why the density does not actually drop with distance outside the mesh.

This also explain the nutty behavior at the target. It is essentially just copying the target value to all nearest positions.

Artifacts within the mesh

The artifacts inside the mesh are coming because this is a 2D interpolation using values from the four closest points. Thus a "X" like feature will appear if the values at nearest points are very different from each other.

The other extreme would be just giving the whole cell the same value everywhere inside which will create discontinuities at the boundary.

So I am not sure how we can not have such artifacts in places where the gradient is large. My guesses are:

• Use an analytical model for the distribution and fitting the data to that model.
• Using gradients in nearby cells to do a local fiting, something similar to cubic spline in 1D, to get a more general way of interpolation.

Possible flaw in weighed mean

I cross-checked the code and there does not seem to be an error in it, atleast to me. But I might be blind to an error I committed myself, so if you think that the reproduced mesh values are too high, maybe you should check closely once.

In the core region, edge_profile extrapolation agreed with the core_profile values which again makes the case stronger for edge_profile interpolation to be correct.

[anchal-physics]

Thin-plate spline

Well, well, well...

GGDUtils_interp_from_edge_profiles_thin_plates.pdf

I used the thin plate spline method generously defined in this document (distributed under Creative Commons license): http://www.geometrictools.com/Documentation/ThinPlateSplines.pdf

It is basically a surface solution that reduces the bending energy while remaining on all the specified grid points. The above plot has enforced color bar limits to show the variation properly in mesh and core area. Without such limits though, we can see how the extrapolation goes outside the grid:

GGDUtils_interp_from_edge_profiles_thin_plates_complete.pdf

P.S. Use the pdf if you want to zoom into some area.

[eldond]

Good work, this is great for a simple approximation. This will take us pretty far with testing. It's still a little excitable near the targets, but that's okay for now. When I finish off the extended mesh and my fancy witchcraft for extrapolating into the SOL, this same interpolation scheme can draw from that extension.