ProjectTorreyPines / SOLPS2ctrl.jl

Synthetic diagnostic workflow manager for use with SOLPS models
Apache License 2.0
2 stars 0 forks source link

core_profile_2d only works inside equilibrium grid #20

Closed anchal-physics closed 11 months ago

anchal-physics commented 1 year ago

I found out that core_profile_2d is really only interpolating within the equilibrium profile grid. It results in an error if a point outside the grid is provided. I think it should return 0 outside the grid instead of failing, wasn't that what we wanted?

Reproduce

cd SD4SOLPS.jl
git pull
git checkout issue_20
julia --project=./ test/test_core_profile_2d.jl

Looking deeper

SD4SOLPS is using Interpolations.jl at v0.13.6 while the latest version is v0.14.7. We are using the convenience constructor LinearInterpolation which is not mentioned in the documentation anymore. The documentation suggests using linear_interpolation constructor that has an argument for setting extrapolation boundary condition using extrapolation_bc keyword. This can be set to 0 to make the value 0 outside everywhere. But the reason out Interpolations.jl is stuck at v0.13.6 is because the package NumericalIntegration.jl pins it to that. This package has unknown build status and was updated 2 years ago. We are only using it for a cumulative integration, maybe we should remove dependency on this package and use something newer and up to date, and we can update Interpolations.jl to use extrapolation feature after that.

orso82 commented 1 year ago

for numerical integration take a look at QuadGK.jl https://juliamath.github.io/QuadGK.jl/stable/

anchal-physics commented 1 year ago

Yeah I'm using QuadGK.jl too in SynthDiag.jl as most people online suggest it. But it requires a function input which can compute the value at arbitrary points. I think that wasn't the use case here that's why NumericalIntergration.jl was used. But I think creating such a function should not be hard.

eldond commented 1 year ago

Technically, core profiles are only defined within the plasma boundary, so this isn't wrong, exactly, but yes it's limiting. 0 isn't necessarily the correct value outside of the boundary. We do have some SOLPS data in the SOL that should be used when a point is within the SOLPS mesh but outside the main plasma boundary. The edge extrapolation solution should be used beyond that. Since edge extrapolation isn't done yet, it's safe to estimate electron density = 0 outside the SOLPS mesh.

The core_profile_2d() function is implementing a simple feature, that while inconveniently narrow in scope, isn't technically incorrect in what it's doing. The 1D profiles shouldn't be swept around into 2D on the open flux surfaces because the approximation that plasma parameters like T_e and n_e are constant on a flux surface only holds on the closed flux surfaces; we know there's significant poloidal variation in the SOL. Perhaps core_profile_2d() can be called by a different, more convenient function that gets a value anywhere in R,Z, sometimes calling core_profile_2d() and sometimes doing something else.