Sensei2 ghost libsim

Jacobi/Sensei/C/Bridge.cxx

@@ -21,29 +21,26 @@ namespace BridgeInternals
 }
 
 //-----------------------------------------------------------------------------
-void bridge_initialize(MPI_Comm comm,
-  int m, int rankx, int ranky, int bx, int by, int ng,
-  const char* config_file)
+void bridge_initialize(MPI_Comm comm, simulation_data *sim, const char *config_file)
 {
   BridgeInternals::comm = comm;
 
   // initialize the data adaptor
   BridgeInternals::GlobalDataAdaptor =
     vtkSmartPointer<pjacobi::JacobiDataAdaptor>::New();
 
-  BridgeInternals::GlobalDataAdaptor->Initialize(m, rankx, ranky, bx, by, ng);
+  BridgeInternals::GlobalDataAdaptor->Initialize(sim);
 
   // initialize the analysis adaptor
   BridgeInternals::GlobalAnalysisAdaptor = vtkSmartPointer<sensei::ConfigurableAnalysis>::New();
   BridgeInternals::GlobalAnalysisAdaptor->Initialize(comm, config_file);
 }
 
 //-----------------------------------------------------------------------------
-void bridge_update(int tstep, double time, double* temperature)
+void bridge_update(simulation_data *sim)
 {
-  BridgeInternals::GlobalDataAdaptor->SetDataTime(time);
-  BridgeInternals::GlobalDataAdaptor->SetDataTimeStep(tstep);
-  BridgeInternals::GlobalDataAdaptor->AddArray("temperature", temperature);
+  BridgeInternals::GlobalDataAdaptor->Update(sim);
+
   if (!BridgeInternals::GlobalAnalysisAdaptor->Execute(BridgeInternals::GlobalDataAdaptor))
     {
     cerr << "ERROR: Failed to execute the analysis" << endl;

Jacobi/Sensei/C/Bridge.h

@@ -3,22 +3,21 @@
 
 #include <mpi.h>
 #include <stdint.h>
+#include <solvers.h>
 
 #ifdef __cplusplus
 extern "C" {
 #endif
   /// This defines the analysis bridge for the pjacobi miniapp.
 
   /// Called before simulation loop
-  void bridge_initialize(MPI_Comm comm,
-    int m, int rankx, int ranky, int bx, int by, int ng,
-    const char* config_file);
+  void bridge_initialize(MPI_Comm comm, simulation_data *sim, const char *config);
 
   /// Called per timestep in the simulation loop
-  void bridge_update(int tstep, double time, double* temperature);
+  void bridge_update(simulation_data *sim);
 
   /// Called just before simulation terminates.
-  void bridge_finalize();
+  void bridge_finalize(void);
 
 #ifdef __cplusplus
 } // extern "C"

Jacobi/Sensei/C/jacobi_sensei_2.xml

@@ -10,12 +10,13 @@
 
   <!-- CATALYST -->
   <analysis type="catalyst" pipeline="pythonscript"
-    filename="jacobi_catalyst_sensei_2.py" enabled="1" />
+    filename="jacobi_catalyst_sensei_2.py" enabled="0" />
 
   <!-- LIBSIM -->
-  <analysis type="libsim" plots="Pseudocolor,Mesh" plotvars="temperature,mesh"
+  <analysis type="libsim" visitdir="/Users/bjw/Development/MAIN/trunk/install" mode="interactive,paused" options="-debug 5 -clobber_vlogs" trace="trace"
+    plots="Pseudocolor,Mesh" plotvars="temperature,mesh" 
     image-filename="jacobi_%ts" image-width="800" image-height="800"
-    slice-project="1" image-format="png" enabled="0"/>
+    slice-project="1" image-format="png" enabled="1"/>
 
   <!-- ADIOS -->
   <analysis type="adios" filename="jacobi.bp" method="MPI" enabled="0" />

Jacobi/Sensei/C/pjacobi.c

@@ -53,16 +53,15 @@ int main(int argc, char *argv[])
 
   const char *config = argv[1];
 
-  bridge_initialize(MPI_COMM_WORLD, sim.m,
-    sim.rankx, sim.ranky, sim.bx, sim.by, 1, config);
+  bridge_initialize(MPI_COMM_WORLD, &sim, config);
 #endif
 
 
   while ((sim.gdel > TOL) && (sim.iter <= MAXSTEPS))
     {  // iterate until error below threshold
     simulate_one_timestep(&sim);
 #ifdef ENABLE_SENSEI
-   bridge_update(sim.iter, sim.iter*1.0, sim.Temp);
+   bridge_update(&sim);
 #endif
     }
 

Jacobi/Sensei/C/solution/JacobiDataAdaptor.cxx

@@ -7,6 +7,9 @@
 #include <vtkInformation.h>
 #include <vtkObjectFactory.h>
 #include <vtkPointData.h>
+#include <vtkFloatArray.h>
+#include <vtkRectilinearGrid.h>
+#include <vtkUnsignedCharArray.h>
 
 #include <Error.h>
 #if defined(SENSEI_2)
@@ -23,6 +26,7 @@ senseiNewMacro(JacobiDataAdaptor);
 //-----------------------------------------------------------------------------
 JacobiDataAdaptor::JacobiDataAdaptor()
 {
+    this->sim = NULL;
 }
 
 //-----------------------------------------------------------------------------
@@ -31,24 +35,19 @@ JacobiDataAdaptor::~JacobiDataAdaptor()
 }
 
 //-----------------------------------------------------------------------------
-void JacobiDataAdaptor::Initialize(int m, int rankx, int ranky, int bx, int by, int ng)
+void JacobiDataAdaptor::Initialize(simulation_data *s)
 {
-  (void)m; // m could be used for whole extents [0 m+2 0 m+2 0 0]
-
-  this->Origin[0] = -double(ng);
-  this->Origin[1] = -double(ng);
-  this->Origin[2] = 0.0;
-
-  this->Spacing[0] = 1.0;
-  this->Spacing[1] = 1.0;
-  this->Spacing[2] = 1.0;
-
-  this->Extent[0] = rankx*(bx - 1);
-  this->Extent[1] = this->Extent[0] + bx + 2*ng - 1;
-  this->Extent[2] = ranky*by;
-  this->Extent[3] = this->Extent[2] + by + 2*ng - 1;
-  this->Extent[4] = 0;
-  this->Extent[5] = 0;
+    this->sim = s;
+}
+
+//-----------------------------------------------------------------------------
+void JacobiDataAdaptor::Update(simulation_data *s)
+{
+    this->sim = s;
+
+    this->SetDataTime(sim->iter*1.);
+    this->SetDataTimeStep(sim->iter);
+    this->AddArray("temperature", sim->Temp);
 }
 
 //-----------------------------------------------------------------------------
@@ -111,14 +110,30 @@ int JacobiDataAdaptor::GetMesh(const std::string &meshName, bool structureOnly,
     int n_ranks = 1;
     MPI_Comm_size(MPI_COMM_WORLD, &n_ranks);
 
-    vtkImageData *block = vtkImageData::New();
-    block->SetExtent(this->Extent);
-    block->SetOrigin(this->Origin);
-    block->SetSpacing(this->Spacing);
-
+    vtkRectilinearGrid *block = vtkRectilinearGrid::New();
+    int dims[3];
+    dims[0] = this->sim->bx+2;
+    dims[1] = this->sim->by+2;
+    dims[2] = 1;
+    block->SetDimensions(dims);
+
+    vtkFloatArray *x = vtkFloatArray::New();
+    vtkFloatArray *y = vtkFloatArray::New();
+    vtkFloatArray *z = vtkFloatArray::New();
+    int dontDelete = 1;
+    x->SetArray(this->sim->cx, this->sim->bx+2, dontDelete);
+    y->SetArray(this->sim->cy, this->sim->by+2, dontDelete);
+    z->SetNumberOfTuples(1);
+    z->SetTuple1(0,0.);
+    block->SetXCoordinates(x);
+    block->SetYCoordinates(y);
+    block->SetZCoordinates(z);
+    x->Delete();
+    y->Delete();
+    z->Delete();
     this->Mesh = vtkSmartPointer<vtkMultiBlockDataSet>::New();
-    this->Mesh->SetNumberOfBlocks(n_ranks);
-    this->Mesh->SetBlock(rank, block);
+    this->Mesh->SetNumberOfBlocks(1);
+    this->Mesh->SetBlock(0, block);
 
     block->Delete();
     }
@@ -135,14 +150,14 @@ int JacobiDataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName
   vtkMultiBlockDataSet *mb = dynamic_cast<vtkMultiBlockDataSet*>(mesh);
   if (!mb)
     {
-    SENSEI_ERROR("Invlaid mesh type " << mesh->GetClassName())
+    SENSEI_ERROR("Invalid mesh type " << mesh->GetClassName())
     return -1;
     }
 
   int rank = 0;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
-  vtkImageData *block = dynamic_cast<vtkImageData*>(mb->GetBlock(rank));
+  vtkRectilinearGrid *block = dynamic_cast<vtkRectilinearGrid*>(mb->GetBlock(0));
   if (!block)
     return 0;
 
@@ -172,10 +187,7 @@ int JacobiDataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName
     vtkSmartPointer<vtkDoubleArray>& vtkarray = this->Arrays[iterV->first];
     vtkarray = vtkSmartPointer<vtkDoubleArray>::New();
     vtkarray->SetName(arrayName.c_str());
-
-    vtkIdType size = (this->Extent[1] - this->Extent[0] + 1) *
-      (this->Extent[3] - this->Extent[2] + 1) * (this->Extent[5] - this->Extent[4] + 1);
-
+    vtkIdType size = (this->sim->bx+2) * (this->sim->by+2);
     vtkarray->SetArray(iterV->second, size, 1);
 
     block->GetPointData()->SetScalars(vtkarray);
@@ -232,6 +244,83 @@ int JacobiDataAdaptor::GetArrayName(const std::string &meshName, int association
   return 0;
 }
 
+int JacobiDataAdaptor::GetMeshHasGhostNodes(const std::string &meshName, 
+  bool &hasGhostNodes, int &nLayers)
+{
+  if (meshName != "mesh")
+    {
+    hasGhostNodes = false;
+    nLayers = 0;
+    SENSEI_ERROR("No mesh named " << meshName)
+    return -1;
+    }
+
+  hasGhostNodes = true;
+  nLayers = 1;
+  return 0;
+}
+
+int JacobiDataAdaptor::AddGhostNodesArray(vtkDataObject* mesh, const std::string &meshName)
+{
+  if (meshName != "mesh")
+    {
+    SENSEI_ERROR("No mesh named " << meshName)
+    return -1;
+    }
+
+  vtkMultiBlockDataSet *mb = dynamic_cast<vtkMultiBlockDataSet*>(mesh);
+  if (!mb)
+    {
+    SENSEI_ERROR("Invalid mesh type " << mesh->GetClassName())
+    return -1;
+    }
+
+  vtkRectilinearGrid *block = dynamic_cast<vtkRectilinearGrid*>(mb->GetBlock(0));
+  if (!block)
+    return -1;
+
+  int nx = this->sim->bx+2;
+  int ny = this->sim->by+2;
+  vtkUnsignedCharArray *gn = vtkUnsignedCharArray::New();
+  gn->SetNumberOfTuples(nx*ny);
+  gn->SetName(GHOST_NODE_ARRAY_NAME().c_str());
+  unsigned char *gptr = (unsigned char *)gn->GetVoidPointer(0);
+  memset(gn->GetVoidPointer(0), 0, nx*nx*sizeof(unsigned char));
+  unsigned char ghost = 1;
+  // Left column
+  if(this->sim->rankx > 0)
+    {
+    int i = 0, j = 0;
+    for( ; j < ny; ++j)
+        gptr[j*nx+i] = ghost;
+    }
+  // Right column
+  if(this->sim->rankx < this->sim->cart_dims[0]-1)
+    {
+    int i = nx-1, j = 0;
+    for( ; j < ny; ++j)
+        gptr[j*nx+i] = ghost;
+    }
+  // Bottom row
+  if(this->sim->ranky > 0)
+    {
+    int i = 0, j = 0;
+    for( ; i < nx; ++i)
+        gptr[j*nx+i] = ghost;
+    }
+  // Top row
+  if(this->sim->ranky < this->sim->cart_dims[1]-1)
+    {
+    int i = 0, j = ny-1;
+    for( ; i < nx; ++i)
+        gptr[j*nx+i] = ghost;
+    }
+  block->GetPointData()->SetScalars(gn);
+  gn->Delete();
+
+  return 0;
+}
+
 //-----------------------------------------------------------------------------
 int JacobiDataAdaptor::ReleaseData()
 {
@@ -252,17 +341,34 @@ vtkDataObject* JacobiDataAdaptor::GetMesh(bool vtkNotUsed(structure_only))
     int n_ranks = 1;
     MPI_Comm_size(MPI_COMM_WORLD, &n_ranks);
 
-    vtkImageData *block = vtkImageData::New();
-    block->SetExtent(this->Extent);
-    block->SetOrigin(this->Origin);
-    block->SetSpacing(this->Spacing);
-
+    vtkRectilinearGrid *block = vtkRectilinearGrid::New();
+    int dims[3];
+    dims[0] = this->sim->bx+2;
+    dims[1] = this->sim->by+2;
+    dims[2] = 1;
+    block->SetDimensions(dims);
+
+    vtkFloatArray *x = vtkFloatArray::New();
+    vtkFloatArray *y = vtkFloatArray::New();
+    vtkFloatArray *z = vtkFloatArray::New();
+    int dontDelete = 1;
+    x->SetArray(this->sim->cx, this->sim->bx+2, dontDelete);
+    y->SetArray(this->sim->cy, this->sim->by+2, dontDelete);
+    z->SetNumberOfTuples(1);
+    z->SetTuple1(0,0.);
+    block->SetXCoordinates(x);
+    block->SetYCoordinates(y);
+    block->SetZCoordinates(z);
+    x->Delete();
+    y->Delete();
+    z->Delete();
     this->Mesh = vtkSmartPointer<vtkMultiBlockDataSet>::New();
-    this->Mesh->SetNumberOfBlocks(n_ranks);
-    this->Mesh->SetBlock(rank, block);
+    this->Mesh->SetNumberOfBlocks(1);
+    this->Mesh->SetBlock(0, block);
 
     block->Delete();
     }
+
   return this->Mesh;
 }
 
@@ -283,7 +389,7 @@ bool JacobiDataAdaptor::AddArray(vtkDataObject* mesh, int association, const std
   int rank = 0;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   vtkMultiBlockDataSet *mb = dynamic_cast<vtkMultiBlockDataSet*>(mesh);
-  vtkImageData *block = mb ? dynamic_cast<vtkImageData*>(mb->GetBlock(rank)) : nullptr;
+  vtkRectilinearGrid *block = mb ? dynamic_cast<vtkRectilinearGrid*>(mb->GetBlock(0)) : nullptr;
   if (!block)
     return false;
 
@@ -293,8 +399,7 @@ bool JacobiDataAdaptor::AddArray(vtkDataObject* mesh, int association, const std
     vtkSmartPointer<vtkDoubleArray>& vtkarray = this->Arrays[iterV->first];
     vtkarray = vtkSmartPointer<vtkDoubleArray>::New();
     vtkarray->SetName(name.c_str());
-    vtkIdType size = (this->Extent[1] - this->Extent[0] + 1) *
-      (this->Extent[3] - this->Extent[2] + 1) * (this->Extent[5] - this->Extent[4] + 1);
+    vtkIdType size = (this->sim->bx+2) * (this->sim->by+2);
     vtkarray->SetArray(iterV->second, size, 1);
     block->GetPointData()->SetScalars(vtkarray);
     return true;