VisIt failed to show the LAST vertical layer of some variables from the FVCOM NetCDF results

[SiqiLiOcean]

I am a ocean model FVCOM user and I love using VisIt to check the FVCOM outputs (in NetCDF). But I found that VisIt failed to show the last vertical layer of some variables, when drawing pseudocolor. I wish to solve this problem and hope anyone could give me some suggestion.

There are three kinds of variables in the FVCOM files.
1> in size of (node, siglay, time), such as 'temp'. For the second dimension, the last layer (j=siglay) could not be drawn. (Please try the 'temp' in the attached data)
2> in size of (nele, siglay, time), such as 'u'. All the layers could be drawn.
3> in size of (node, siglev, time), such as 'omega'. For the second dimension, the last layer (j=siglev) could not be drawn. (Please try the 'omega' in the attached data)

Please download the test data from:
https://drive.google.com/file/d/19GYoOWd8Vqss5OkEFXigFYmw2k7qih-b/view?usp=sharing

Any suggestion would be helpful.

Thanks.

Siqi

[JustinPrivitera]

Hello @SiqiLiOcean. I'm not sure what you're referring to when you say VisIt failed to show the last vertical layer of some variables. I downloaded the attached data and opened it in the latest version of VisIt, 3.4.2. I made plots of temp, u, and omega. I see the following data extents:

What should be different about these plots? Are we missing data that exists in the files?

[SiqiLiOcean]

Thank you for your reply.

Take the variable 'temp' for an example. This variable is in size of (5620, 10, 4), where 5620 is the horizontal grid size, 10 is the vertical layer number, 4 is the time number. So there should be 10 layers, but when I use visit to see it, there is only 9 layers.

[JustinPrivitera]

I found these two code snippets in the reader. This first one is where the material set object is created which is what you see when you turn on and off individual sigma layers. It looks like it is using the number of sigma layers minus 1 on purpose:

visit/src/databases/NETCDF/avtFVCOMReader.C

Lines 1875 to 1894 in 457f046

	if(mesh4==true)
	{
	//Add Material MetaData for: SigmaLayer_Mesh
	avtMaterialMetaData *matmd_lay = new avtMaterialMetaData;
	matmd_lay->name = "Sigma_Layers";
	matmd_lay->meshName = "SigmaLayer_Mesh";
	matmd_lay->numMaterials=(nSiglay-1);
	for (int i=0; i < (nSiglay-1); ++i)
	{
	char buffer[50];
	int n;
	n=sprintf(buffer, "Layer %d",(i+1)); (void) n;
	matmd_lay->materialNames.push_back(buffer);
	}
	md->Add(matmd_lay);
	}

This next snippet has a comment explaining why there is one less sigma layer in the sigma layers mesh:

visit/src/databases/NETCDF/avtFVCOMReader.C

Lines 2787 to 2790 in 457f046

	int nzones=(nMats)*nElem; // there are nSigLayers
	// in the levels mesh. There are nSigLayers-1 in the Sigma_Layers
	// mesh because its nodes are at the edge centers of the original
	// mesh.

It isn't clear to me why it is this way. I will confer with other developers. Perhaps that comment would make more sense to you as a user. From your opening comment in this discussion, it sounds like you are saying that the first code block is wrong, and we should not subtract 1 from the number of sigma layers.

[markcmiller86]

@SiqiLiOcean I wrote a ton of stuff here explaining the situation before I realized...your problem is with the SIL controls...not the way VisIt is displaying the variables. But, I don't want to loose what I wrote so am posting this now.

VisIt displays data on meshes. But, a NETCDF file does not define meshes. NETCDF defines multi-dimensional arrays of numbers. Some arrays are #nodes x #layers x #times. Some arrays are #elements x #layers x #times. Some arrays are #nodes x #levels x #times.

For each of these NETCDF arrays, VisIt has a choice to make...is that array defining a node-centered quantity or a zone-centered quantity? For the same array of numbers, treating it as a node-centered quantity invoves constructing a different mesh than treating that same array of numbers as a zone-centerd quantity.

Lets consider an example using a 3x3 array of numbers. For a node-centered treatment, the mesh that gets constructed is 1:1 with that array. In other words, it too has a 3x3 array of coordinates (3 stesp, 0,1,2 in X and 3 steps, 0,1,2 in Y). For a zone-centered treatment of that same 3x3 array, the mesh that gets constructed actually has 4x4 array of coordinates (4 steps, 0,1,2,3 in X and 4 steps 0,1,2,3 in Y). Worse, in order for the two meshes like this to overlay in an reasonable way, they have to be offset from each other by a half-step in each dimension. In other words, the grid crossings of the 3x3 node-mesh should geometricly line up at the same geometric points as the centers of the zones (cells, elements) in the 4x4 zone-mesh.

In the Image reader plugin, we treat wind up defining every variable in the image twice...once as a node-centered quantity and once as a zone-centered quantity. And we construct the two meshes for each case as described above. When the number of variables is small, and the input database does nothing else to inform down-stream consumers how the data should be treated, that is usually works fine because it gives the user the option of deciding which treatment they want during the session.

Without explicit information in the NETCDF file telling us precisely what to do (and having written the code in the plugin to respond to it), we have to do some guesswork on how to handle those variables. That guesswork is programmed into the FVCOM plugin in this case based on either the dimension names of the underlying dimension of each variable or the variable names themselves.

The plugin defines a SigmaLayer_Mesh, and a SigmaLevel_Mesh...

  else if(strcmp(mesh, "SigmaLayer_Mesh") == 0)
    { // SigmaLayer_Mesh
      // Plot scalars which are calculated at cell edge centers
      // On this grid these scalars will be on cell nodes
                                                          
      debug4<< mName << "Getting SigmaLayer_Mesh" << endl;
      //Sigma layer mesh has one fewer layers and one fewer levels!
      int tNodes = nNode * (nSiglay);                      
      int tCells = nElem * (nSiglay-1);
    
.
.
.
  else if(strcmp(mesh, "SigmaLevel_Mesh") == 0)
    {   // is mesh: SigmaLevel_Mesh
      // plot turbulent quantities which occur at cell nodes (fvcom 2.5+)
      // plot turbulent quantities which occur at cell center (fvcom <2.5)
      // [Note cell center km is not at the right level!] 
      
      // plot velocity at cell centers 
      debug4 << mName << "Getting SigmaLevel Mesh" << endl;
      
      int tNodes = nNode * nSiglev;
      int tCells = nElem * (nSiglev-1);

It then apparently has logic to define variables on either of these meshes as either node or zone-centered. Using File Information from a VisIt session viewing the example file provided, sci_exp.nc, the plugin is defining...

• temp as node-centered on the SigmaLayer_Mesh.
  • ncdump shows that variable is defined in the file as float temp(time, siglay, node)
  • A plot of that variable is shown below. siglay is 10 and it has 10 steps in that dimension. I think it is good.
  • However the SIL is apparently also display controls for each layer and that is showing only the first 9 layers. So, the plugin is somehow interpreting the layers here as layers of cells and not *nodes.
• u as zone-centered on the SigmaLevel-Mesh...which would be a problem but because siglay is siglev-1, it works for zone-centering. However, in general, I think this is a problem
  • ncdump shows that variable is defined in the file as float u(time, siglay, nele) which has 10 steps in that dimension but only 9 cells (gaps between the steps).
• omega as node centered on the SigmaLevel_Mesh
  • ncdump shows that variable is defined in the file as float omega(time, siglev, node) with 11 steps in that dimension

Example omega

Example u

[markcmiller86]

@JustinPrivitera I think the issue is that the dimensions of the mesh representing layers and levels are represented both in the mesh objects themselves and also in the SIL controls. In this example in the input file, siglay (the z dimension of these meshes) is 10 but the number of SIL controls is just 9. Again, I think this has to do with whether the SIL controls are intended to peel off layers of zones or layers of nodes.

[markcmiller86]

I think we should and tests to our test suite using this data set.

[JustinPrivitera]

This is really great context @markcmiller86. Thanks so much for weighing in! @SiqiLiOcean, does what Mark said sound reasonable? Can you let us know, based on what he said, how you expect your variables to display? Can we use the data you provided in our test suite?