Strange values when adding cell data to grid

In Situ Support
1

I’m trying to add cell data to a grid. This is what my FECxx.cxx looks like:

    // Adaptor for getting Fortran simulation code into ParaView CoProcessor.

    // CoProcessor specific headers
    #include "vtkCPDataDescription.h"
    #include "vtkCPInputDataDescription.h"
    #include "vtkCPProcessor.h"
    #include "vtkCPPythonScriptPipeline.h"
    #include "vtkDoubleArray.h"
    #include "vtkCellData.h"
    #include "vtkSmartPointer.h"
    #include "vtkRectilinearGrid.h"

    // Fortran specific header
    #include "vtkCPPythonAdaptorAPI.h"

    // These will be called from the Fortran "glue" code"
    // Completely dependent on data layout, structured vs. unstructured, etc.
    // since VTK/ParaView uses different internal layouts for each.

    // Creates the data container for the CoProcessor.
    extern "C" void createcpdata_(int* is, int* ie, int* js, int* je, int* ks, int* ke, int* nx, int* ny, int* nz, double* x, double* y, double* z)
    {
        if (!vtkCPPythonAdaptorAPI::GetCoProcessorData())
        {
            vtkGenericWarningMacro("Unable to access CoProcessorData.");
            return;
        }
        
        // Set whole extent
        vtkCPPythonAdaptorAPI::GetCoProcessorData()->
            GetInputDescriptionByName("input")->SetWholeExtent(0, *nx, 0, *ny, 0, *nz);

        // create a Rectilinear grid
        cout << "FECxx.cxx : Creating Rectilinear Grid\n";
        auto grid = vtkSmartPointer<vtkRectilinearGrid>::New();
        
        grid->SetExtent(*is-1, *ie, *js-1, *je, *ks-1, *ke);
        
        vtkSmartPointer<vtkDoubleArray> xCoords = vtkSmartPointer<vtkDoubleArray>::New();
        xCoords->SetNumberOfComponents(1);
        vtkSmartPointer<vtkDoubleArray> yCoords = vtkSmartPointer<vtkDoubleArray>::New();
        yCoords->SetNumberOfComponents(1);
        vtkSmartPointer<vtkDoubleArray> zCoords = vtkSmartPointer<vtkDoubleArray>::New();
        zCoords->SetNumberOfComponents(1);

        for(int i = *is-1; i <= *ie; i++)
        {
          xCoords->InsertNextValue( x [i] );
        }
        for(int i = *js-1; i <= *je; i++)
        {
          yCoords->InsertNextValue( y [i] );
        }
        for(int i = *ks-1; i <= *ke; i++)
        {
          zCoords->InsertNextValue( z [i] );
        }

        grid->SetXCoordinates(xCoords);
        grid->SetYCoordinates(yCoords);
        grid->SetZCoordinates(yCoords);


        // Name should be consistent between here, Fortran and Python client script.q
        vtkCPPythonAdaptorAPI::GetCoProcessorData()->
            GetInputDescriptionByName("input")->SetGrid(grid);
    }

    // Add field(s) to the data container.
    // Separate from above because this will be dynamic, grid is static.
    // By hand name mangling for fortran.
    extern "C" void addfield_(double* scalars, char* name)
    {
        auto idd = vtkCPPythonAdaptorAPI::GetCoProcessorData()->
                        GetInputDescriptionByName("input");

        auto grid = vtkRectilinearGrid::SafeDownCast(idd->GetGrid());

        if (!grid)
        {
          vtkGenericWarningMacro("No adaptor grid to attach field data to.");
          return;
        }

        // field name must match that in the fortran code.
        if (idd->IsFieldNeeded(name))
        { 
            auto field = vtkSmartPointer<vtkDoubleArray>::New();
            field->SetName(name);
            field->SetNumberOfComponents(1);
            //field->SetNumberOfTuples(grid->GetNumberOfCells());
            field->SetArray(scalars, 1 * grid->GetNumberOfCells(), 1);
            // for debugging =================================================
            for(int i = 0; i < grid->GetNumberOfCells(); i++)
            {
              cout << scalars[i] << " ";
            }
            cout << "\n";
            // for debugging =================================================
            grid->GetCellData()->AddArray(field);
        }
    }

I have the coprocessor write the data to a pvtr file which I’m inspecting in Paraview. My grid is fine, but there is an error in the file addfield_ function that I just can’t find. For some reason the cell data that is written to the pvtr file is strange (holding zeroes and really,really,really huge values). As you can see I am printing out the values of scalars and there they are fine.

Any idea what I’m doing wrong?

2

Maybe instead of:

should just do:

My guess is that the referencing counting is getting messed up with auto.

3

Unfortunately that does not resolve the issue. Still getting the same strange values.

4

You need to make sure that grid->SetExent() is called before addfield_() is called. Check to see what grid->GetNumberOfCells() is set during your addfield_() call.

5

That is the case and grid->GetNumberOfCells() returns the proper number of cells.

6

I’m basically calling such a subroutine in my Fortran code

  subroutine post_catalyst_process( [...], p )

    use CoProcessor

    implicit none
    
    PetscScalar, dimension(gis:gie,gjs:gje,gks:gke) :: p
    
    [...]

    call runcoprocessor( p(is:ie,js:je,ks:ke) )

end subroutine

where p includes ghost cells that I intended not to pass on to the runcoprocessor subroutine. (gis=is-1, gie=ie+1, …)

 module CoProcessor
   use iso_c_binding
   implicit none
   public
   interface CoProcessor_adaptor
       module procedure runcoprocessor
   end interface 

 contains

   subroutine runcoprocessor(is,ie,js,je,ks,ke,nx,ny,nz,x,y,z,step,time,pressure)
     use iso_c_binding
     implicit none
     integer, intent(in) :: is,ie,js,je,ks,ke,nx,ny,nz,step
     real(kind=8), dimension(:), intent(in) :: x,y,z
     real(kind=8), intent(in) :: time
     real(kind=8), dimension(:,:,:), intent (in) :: pressure
     integer :: flag
     call requestdatadescription(step,time,flag)
     if (flag .ne. 0) then
       call needtocreategrid(flag)
       if (flag .ne. 0) then
         call createcpdata(is,ie,js,je,ks,ke,nx,ny,nz,x,y,z)
       end if
       ! adding //char(0) appends the C++ terminating character to the Fortran array
       call addfield(pressure,"pressure"//char(0))
       call coprocess()
     end if
     
     return
     
   end subroutine runcoprocessor

 end module CoProcessor

The trouble appears to be related to the fact that p is an array that is handled by PETSc. If instead I create my own array p inside post_catalyst_process everything is fine.

Not sure how to resolve my issue.

7

I still think it’s puzzling though that

            for(int i = 0; i < grid->GetNumberOfCells(); i++)
            {
              cout << scalars[i] << " ";
            }

prints out the proper values, whereas

field->SetArray(scalars, 1 * grid->GetNumberOfCells(), 1);

introduces strange values.

8

The following needs the memory in scalars to not be freed before Catalyst uses the memory since it doesn’t copy the data but just keeps the pointer to the memory:

Blockquote
field->SetArray(scalars, 1 * grid->GetNumberOfCells(), 1);

9

I’m pretty confident, the memory is not being altered in the meantime. I’m even calling coprocess() just after calling addfield(..)

10

Maybe check what’s in the array after coprocess() is called. This appears to be a memory access issue where I’m suspecting the memory is deleted before Catalyst is using it. Beyond the array itself not being freed you can also look into reference counting of the VTK/ParaView objects to see if something is getting freed before you think it should be.

11

Just the right values.

In my Fortran code I call runcoprocessor
call runcoprocessor(is,ie,js,je,ks,ke,nx,ny,nz,xst(0:nx),yst(0:ny),zst(0:nz),nt,t,p(is:ie,js:je,ks:ke))

 module CoProcessor
   use iso_c_binding
   implicit none
   public
   interface CoProcessor_adaptor
       module procedure runcoprocessor
   end interface 

 contains

   subroutine runcoprocessor(is,ie,js,je,ks,ke,nx,ny,nz,x,y,z,step,time,pressure)
     use iso_c_binding
     implicit none
     integer, intent(in) :: is,ie,js,je,ks,ke,nx,ny,nz,step
     real(kind=8), dimension(:), intent(in) :: x,y,z
     real(kind=8), intent(in) :: time
     real(kind=8), dimension(:,:,:), intent (in) :: pressure
     integer :: flag
     call requestdatadescription(step,time,flag)
     if (flag .ne. 0) then
       call needtocreategrid(flag)
       if (flag .ne. 0) then
         call createcpdata(is,ie,js,je,ks,ke,nx,ny,nz,x,y,z)
       end if
       ! adding //char(0) appends the C++ terminating character to the Fortran array
       write(*,*) "before addfield  ",pressure
       call addfield(pressure,"pressure"//char(0))
       write(*,*) "before coprocess ",pressure
       call coprocess()
       write(*,*) "after  coprocess ",pressure
     end if
     
     return
     
   end subroutine runcoprocessor

 end module CoProcessor

as you can see I’m printing out the values of pressure just before and after calling addfield, and before and after calling coprocess

I’m also printing out the values inside the addfield function before and after calling AddArray

extern "C" void addfield_(double* scalars, char* name)
{
    auto idd = vtkCPPythonAdaptorAPI::GetCoProcessorData()->
                    GetInputDescriptionByName("input");

    auto grid = vtkRectilinearGrid::SafeDownCast(idd->GetGrid());

    if (!grid)
    {
      vtkGenericWarningMacro("No adaptor grid to attach field data to.");
      return;
    }

    // field name must match that in the fortran code.
    if (idd->IsFieldNeeded(name))
    { 
        auto field = vtkSmartPointer<vtkDoubleArray>::New();
        field->SetName(name);
        field->SetNumberOfComponents(1);
        //field->SetNumberOfTuples(grid->GetNumberOfCells());
        field->SetArray(scalars, 1 * grid->GetNumberOfCells(), 1);
        cout << "in addfield before AddArray(field) "; for(int i = 0; i < grid->GetNumberOfCells(); i++) { cout << scalars[i] << " "; } cout << "\n"; //debug
        grid->GetCellData()->AddArray(field);
        cout << "in addfield after  AddArray(field) "; for(int i = 0; i < grid->GetNumberOfCells(); i++) { cout << scalars[i] << " "; } cout << "\n"; //debug
    }
}

These values are always printed out correctly.
And still, the pvtr file created by the processing script contains only really strange values.

This is my coproc.py:

from paraview.simple import *
from paraview import coprocessing

#--------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.


# ----------------------- CoProcessor definition -----------------------

def CreateCoProcessor():
    def _CreatePipeline(coprocessor, datadescription):
        class Pipeline:
            adaptorinput = coprocessor.CreateProducer(datadescription, "input")

            writer = servermanager.writers.XMLPRectilinearGridWriter(Input = adaptorinput)
            coprocessor.RegisterWriter(writer, "pressure_%t.pvtr", 1)
  
        return Pipeline()

    class CoProcessor(coprocessing.CoProcessor):
        def CreatePipeline(self, datadescription):
            self.Pipeline = _CreatePipeline(self, datadescription)

    coprocessor = CoProcessor()
    freqs = {'input': [1]}
    coprocessor.SetUpdateFrequencies(freqs)
    return coprocessor

#--------------------------------------------------------------
# Global variables that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()

#--------------------------------------------------------------
# Enable Live-Visualizaton with ParaView
coprocessor.EnableLiveVisualization(False)


# ---------------------- Data Selection method ----------------------

def RequestDataDescription(datadescription):
    "Callback to populate the request for current timestep"
    global coprocessor
    if datadescription.GetForceOutput() == True:
        # We are just going to request all fields and meshes from the simulation
        # code/adaptor.
        for i in range(datadescription.GetNumberOfInputDescriptions()):
            datadescription.GetInputDescription(i).AllFieldsOn()
            datadescription.GetInputDescription(i).GenerateMeshOn()
        return

    # setup requests for all inputs based on the requirements of the
    # pipeline.
    coprocessor.LoadRequestedData(datadescription)

# ------------------------ Processing method ------------------------

def DoCoProcessing(datadescription):
    "Callback to do co-processing for current timestep"
    global coprocessor

    # Update the coprocessor by providing it the newly generated simulation data.
    # If the pipeline hasn't been setup yet, this will setup the pipeline.
    coprocessor.UpdateProducers(datadescription)

    # Write output data, if appropriate.
    coprocessor.WriteData(datadescription);

    # Write image capture (Last arg: rescale lookup table), if appropriate.
    coprocessor.WriteImages(datadescription, rescale_lookuptable=False)

    # Live Visualization, if enabled.
    coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)

I don’t know how to do this.

12

Ah I think I detected the origin of the problem!

The following fails:

    PetscScalar, dimension(gis:gie,gjs:gje,gks:gke) :: p
    call runcoprocessor(is,ie,js,je,ks,ke,nx,ny,nz,xst(0:nx),yst(0:ny),zst(0:nz),nt,t,p(is:ie,js:je,ks:ke))

while this works fine:

    PetscScalar, dimension(is:ie,js:je,ks:ke) :: p
    call runcoprocessor(is,ie,js,je,ks,ke,nx,ny,nz,xst(0:nx),yst(0:ny),zst(0:nz),nt,t,p(is:ie,js:je,ks:ke))

gis=is-1, gie=ie+1 and so on.

The array p includes one layer of ghost cells so that the volumes corresponding to the local p-fields of each processor overlap. Can I just pass them on to Catalyst? I think that this would solve my issue.

13

Use vtkObjectBase::GetReferenceCount() to get the reference count. Or maybe use vtkObject::DebugOn to check the status of the object.

14

@Andy_Bauer thank you so much for following my posts and the tips you gave.

I now finally have a specific question:

Inside my addfield-function I’m currently adding data like this:

auto field = vtkSmartPointer<vtkDoubleArray>::New();
field->SetName(name);
field->SetNumberOfComponents(1);
field->SetArray(scalars, 1 * grid->GetNumberOfCells(), 1);
grid->GetCellData()->AddArray(field);

The key to solve my issue probably is to use the entire original array (scalars here) with the layer of ghost cells it includes on all 6 sides. So scalars is two entries larger in each dimension than the number of cells that the processor owns.

I need to do something like (thinking pythonic…)

field->SetArray(scalars[1:-1][1:-1][1:-1], 1 * grid->GetNumberOfCells(), 1);

How would one do that?

15

You can’t do that but what you can do is to expand the rectilinear grid to include the ghost cells and then mark the cells as ghosts. Otherwise you’ll need to create a full copy of the data without the values you need.

16

Fantastic, this should work. I am now passing on the entire array including ghost cells and it indeed works fine when running on a single core. When using multiple cores, the corresponding subdomains overlap by two cells, since my subdomains are extended by a layer of ghost cells.

My search for guidance on how to mark ghost cells took me here where I Identified the following commands: (Some pseudo-code in the second line!)

vtkUnsignedCharArray* ghosts = grid->GetCellGhostArray();
for ( vtkIdType cellId in cellIdsOfGhostCells ) {
    ghosts->SetValue(cellId, ghosts->GetValue(cellId) | vtkDataSetAttributes::DUPLICATECELL);
}
vtkDataSetAttributes::GenerateGhostArray

How do I get the corresponding cell ids though (second line)? According to the Catalyst Guide V5.6.0 on page 26/27 I assume that the cells are simply numbered, starting with 0 in the bottom/left/back corner, and then incremented in the logical i direction first, next in j, then k.
If that’s the case I would know how to loop over ghost cells on the six sides of my subdomain.

If there is a simpler way of telling Catalyst that I have a single layer of ghost cells on all six sides (including the corners) I would be happy to use that method.

17

The simplest way I can think of is just looping over all of the cells and marking them. I don’t know of any method that does that automatically based on the number of ghost levels.

18

To conclude this discussion …

The problem here was that I tried to pass only a (center) piece of an array to Catalyst which is not possible, and which leads to memory issues.
The motivation behind my attempt to do so was that my ranks each have overlapping parts of the overall information due to layers of ghost cells.
The key is to pass all information to Catalyst but to then also mark ghost cells as such so that Catalyst knows how to deal with the overlapping information.

I still don’t know how to mark cells as ghost cells, but that’s another question which I’m asking here.