conversion of boundary_separatrix and boundary_secondary_separatrix data in dd 3.42 to contour_tree data in dd 4.0

[jfartaud]

The conversion of the data from substructure time_slice/boundary_separatrix and time_slice/boundary_secondary_separatrix of IDS equilibrium in DD 3.42 to the substructure time_slice/contour_tree of IDS equilibrium in DD 4.x is missing. These information are lost, in particular the outline contour, the poloidal flux, also data relative on special points like O and X points and strike points. The moments (R,a, K, ...) for close part of the contour are not included in DD 4.x as it is relevant only for the boundary.

Can-you add it ?

[olivhoenen]

@imbeauf @Simon-McIntosh can we summarize what information from DD3 in boundary_separatrix can have a well defined transfer into contour_tree?

[jfartaud]

For both boundary_separatrix and boundary_secondary_separatrix ?

[jfartaud]

From the adaptation of equilibrium code from dd 3.42 to dd 4.0, most of the information for the separatrix and second separatrix can be transferred. In addition the contour_tree can also include the magnetic axis information.

[imbeauf]

Here is the information that we could retrieve from DDv3 datasets.
We understood that the first element of countour_tree node should be the magnetic axis. If this is correct, then the algorithm would be as follows:
If v3 boundary_second_separatrix/psi is not empty, then contour_tree.node.resize(3)
Elseif v3 boundary_separatrix/psi is not empty, then contour_tree.node.resize(2)
Else contour_tree.node.resize(1)

Node [0]:

• critical_type = 0
• r,z to be taken from global_quantities/magnetic_axis/r,z
• psi to be taken from global_quantities/psi_axis WITH A CHANGE OF SIGN

Node[1]: (if v3 boundary_separatrix/psi is not empty)

• If boundary_separatrix/type = 0 (limiter) critical_type =2, (r,z) = active_limiter_point/r,z, psi = - psi NOTE THE CHANGE OF SIGN
• else critical_type = 1; (r,z) = x_point[0]/r, z; psi = -psi NOTE THE CHANGE OF SIGN
• levelset/r,z = outline/r,z

Node[2]: (If v3 boundary_second_separatrix/psi is not empty)

• critical_type = 1; (r,z) = x_point[0]/r,z; psi = -psi NOTE THE CHANGE OF SIGN
• levelset/r,z = outline/r,z

@maarten-ic , can you add this logic in IMPAS-Python DDv3 to DDv4 specific conversion of the Equilibrium IDS ?

[maarten-ic]

Thanks @imbeauf, I can pick this up. I do have a couple of questions though:

1. The magnetic axis is an O-point, but not necessarily a local minimum of the flux map (critical_type=0). It could also be a maximum (depending on the sign of Ip), right?
2. Same question applies for using critical_type=2 for limiter boundaries.
3. What to do when there are multiple x-points defined in boundary_(second_)separatrix/x_point?
4. There's still quite a number of quantities not converted: geometric parameters, strike points, closest wall point, dr_dz_zero point and gaps. Should we emit a warning (if these are filled) that this data is lost?

[maarten-ic]

And some more questions:
5. Should the contour_tree be filled from boundary when boundary/psi_norm == 1 and both boundary_separatrix and boundary_second_separatrix are not filled?
6. Do you have some example IDSs that we can run the conversion on to test if it's all working properly?

[jfartaud]

answer to 6: You can try to use /home/ITER/hoeneno/public/dina3420 (MAS-AL-Matlab/5.3.0) that has been generated by Olivier Hoenen Olivier for the test of PDS. This is the test case that have trigger the current work.

[jfartaud]

1. and 2. : Yes you are right.
2. I think as the structure is "point oriented", we should have an entry per X-point.
3. the quantities closest wall point, dr_dz_zero point and gaps are relevant only for the boundary used for plasma position control and are already filled if I am right. The strike point information is lost as a standalone information but can be retrieved from levelset points information in conjunction with information in wall or divertor. In this case the edge information will be useful but I think this information is not available in DD3. Maybe it will be useful to develop a tool to compute the strike point?

[maarten-ic]

Thanks @jfartaud.

1. and 2. : Yes you are right.

I've added some logic to detect if the magnetic axis is a (local) minimum or maximum, by comparing global_quantities.psi_axis with global_quantities.psi_boundary (taking into account the COCOS sign flip)

I think as the structure is "point oriented", we should have an entry per X-point.

Right! I've updated the logic to add a node (at the end) for every extra x-point.

the quantities closest wall point, dr_dz_zero point and gaps are relevant only for the boundary used for plasma position control and are already filled if I am right. The strike point information is lost as a standalone information but can be retrieved from levelset points information in conjunction with information in wall or divertor.

Okay clear. I'm not putting in any warnings for these.

answer to 6: You can try to use /home/ITER/hoeneno/public/dina3420 (MAS-AL-Matlab/5.3.0) that has been generated by Olivier Hoenen Olivier for the test of PDS. This is the test case that have trigger the current work.

I've converted the equilibrium IDS in this data entry to 4.0.0 (based on the current state of the PR #62 ) and stored it in imas:hdf5?path=/home/ITER/sebregm/public/dina3420-converted400. Could you have a look at this as well?

Some things I noticed while comparing results with the original dataset:

1. DINA apparently fills global_quantities.psi_axis, boundary_separatrix.psi and boundary_secondary_separatrix.psi even when there is no plasma (as in the first time slice). This results in nodes in the contour tree for many time slices with psi = 0, r = 0 and z = 0. I think the conversion logic is working as expected here, and DINA should not fill the psi values if there is no magnetic axis / separatrix / secondary separatrix.
2. In later time slices (e.g. 2000), DINA has filled boundary_separatrix (with type == 1) and boundary_secondary_separatrix. However, the x-points have r,z = (0,0). So this is also what's copied into the contour tree nodes.

[jfartaud]

I try to read path=/home/ITER/sebregm/public/dina3420-converted400 with Viz and it look like there is some problem with the node 3 that should contain the X-point position and flux but also the (R,Z) set of point of the LCFS coming from boundary_separatrix and don't seem to have been copied. If I understand correctly node 1 is the magnetic axis and it look ok and the node 2 is the limiteur contact point and it look also ok.
I don't know if there is a second boundary in the original dataset?

_Some things I noticed while comparing results with the original dataset:

1. DINA apparently fills global_quantities.psi_axis, boundary_separatrix.psi and boundary_secondary_separatrix.psi even when there is no plasma (as in the first time slice). This results in nodes in the contour tree for many time slices with psi = 0, r = 0 and z = 0. I think the conversion logic is working as expected here, and DINA should not fill the psi values if there is no magnetic axis / separatrix / secondary separatrix._

Yes DINA return data when there is not yet a plasma and when there is no more a plasma. I don't know if is is useful for equilibrium IDS. When we use the DINA data we have to filter it. This come from the fact that DINA solve the electromagnetic circuit equation for the pre-magnetization phase and plasma breakdown and also for the return of coil currents to 0 after the end of the plasma.

3. In later time slices (e.g. 2000), DINA has filled boundary_separatrix (with type == 1) and boundary_secondary_separatrix. However, the x-points have r,z = (0,0). So this is also what's copied into the contour tree nodes.

I think it is fine as we don't use these data.

[maarten-ic]

I try to read path=/home/ITER/sebregm/public/dina3420-converted400 with Viz and it look like there is some problem with the node 3 that should contain the X-point position and flux but also the (R,Z) set of point of the LCFS coming from boundary_separatrix and don't seem to have been copied. If I understand correctly node 1 is the magnetic axis and it look ok and the node 2 is the limiteur contact point and it look also ok.
I don't know if there is a second boundary in the original dataset?

I'm not sure if I understand your question correctly, but perhaps you can have a look at time_slice(1963) (with time=168.90 s)?
DINA has filled both boundary_separatrix and boundary_secondary_separatrix for this time slice. This results in the following contour_tree after conversion. The (R,Z) points of the LCFS and secondary separatrix are stored in the levelset/r and levelset/z arrays in the second respectively third node:

>>> imas.util.print_tree(conv["time_slice(1963)/contour_tree"])
contour_tree
├── time_slice[1962]/contour_tree/node[0]
│   ├── critical_type: 0
│   ├── r: 5.764929561824614
│   ├── z: 0.5564499652769378
│   └── psi: 30.132192769171176
├── time_slice[1962]/contour_tree/node[1]
│   ├── critical_type: 1
│   ├── r: 4.403807113972924
│   ├── z: -1.183538925347541
│   ├── psi: 36.66627927264713
│   └── time_slice[1962]/contour_tree/node[1]/levelset
│       ├── r: array([6.5625, 6.5503, 6.4688, ..., 6.75  , 6.6686, 6.6562])
│       └── z: array([-0.4596, -0.4688, -0.5267, ..., -0.3026, -0.375 , -0.3854])
└── time_slice[1962]/contour_tree/node[2]
    ├── critical_type: 1
    ├── r: 4.041587294434299
    ├── z: 2.088147562270086
    ├── psi: 36.92191181422619
    └── time_slice[1962]/contour_tree/node[2]/levelset
        ├── r: array([3.076 , 3.0938, 3.1103, ..., 3.75  , 3.8438, 3.9375])
        └── z: array([6.    , 5.9511, 5.9062, ..., 2.0834, 2.0846, 2.0861])

[olivhoenen]

In order to move forward with this implementation, maybe we should test the conversion on an equilibrium that you know better than this DINA case @jfartaud?

@Simon-McIntosh could you also give it a try with some data and comment, as you were the main designer of this contour_tree structure?

[jfartaud]

What @maarten-ic show is correct. But I can't see it in Viz.
When I read it in Matlab all is OK.
So maybe I don't use the right version of Viz.