Commits

Lisandro Dalcin committed 43b015a

Some work to get collocation working in parallel

Comments (0)

Files changed (5)

 runex6c_4:
 	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
 runex6d_1:
-	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady -iga_collocation
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_2:
+	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_4:
+	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_8:
+	-@${MPIEXEC} -n 8 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_9:
+	-@${MPIEXEC} -n 9 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
 runex7a_1:
 	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -dim 1
 runex7a_4:
 Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
 Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
 Neumann := Neumann.PETSc runex7a_1 runex7a_4 runex7b_1 runex7b_4 runex7c_4 Neumann.rm
-Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 runex6d_1 Bratu.rm
+Bratu   := \
+Bratu.PETSc \
+runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 \
+runex6d_1 runex6d_2 runex6d_4 runex6d_8 runex6d_9 \
+Bratu.rm
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
 PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 ElasticRod := ElasticRod.PETSc runex8_1 runex8_4 ElasticRod.rm
   PetscInt  width[3];
   PetscInt  sizes[3];
   PetscInt  ID[3];
+  IGABasis  BD[3];
   PetscBool atboundary;
   PetscInt  boundary_id;
   /**/
   PetscInt nsd;
   PetscInt npd;
 
-  IGABasis *BD;
-
   PetscInt    *mapping;   /*[nen]      */
   PetscBool   geometry;
   PetscReal   *geometryX; /*[nen][nsd] */
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveBool(iga,collocation,2);
-  if (collocation) {
-    PetscMPIInt size = 1;
-    ierr = MPI_Comm_size(((PetscObject)iga)->comm,&size);CHKERRQ(ierr);
-    if (size > 1) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
-                          "Collocation not supported in parallel");
-  }
   if (collocation && iga->setup) {
     PetscInt i, dim = iga->dim;
     PetscBool periodic = PETSC_FALSE;
   PetscFunctionReturn(0);
 }
 
+EXTERN_C_BEGIN
+extern PetscReal IGA_Greville(PetscInt i,PetscInt p,const PetscReal U[]);
+extern PetscInt  IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal U[]);
+EXTERN_C_END
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGASetUp_Stage1"
 PetscErrorCode IGASetUp_Stage1(IGA iga)
 {
-  PetscInt       i;
+  PetscInt       i,dim;
+  PetscInt       grid_sizes[3] = {1,1,1};
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   if (iga->setupstage >= 1) PetscFunctionReturn(0);
   iga->setupstage = 1;
 
-  for (i=0; i<iga->dim; i++)
+  dim = iga->dim;
+
+  for (i=0; i<dim; i++)
     {ierr = IGAAxisSetUp(iga->axis[i]);CHKERRQ(ierr);}
-  for (i=iga->dim; i<3; i++)
+  for (i=dim; i<3; i++)
     {ierr = IGAAxisReset(iga->axis[i]);CHKERRQ(ierr);}
 
+  if (!iga->collocation)
+    for (i=0; i<dim; i++)
+      grid_sizes[i] = iga->axis[i]->nel;
+  else
+    for (i=0; i<dim; i++)
+      grid_sizes[i] = iga->axis[i]->nnp;
+
   { /* processor grid and coordinates */
     MPI_Comm    comm = ((PetscObject)iga)->comm;
     PetscMPIInt size,rank;
-    PetscInt    grid_sizes[3] = {1,1,1};
     PetscInt    *proc_sizes = iga->proc_sizes;
     PetscInt    *proc_ranks = iga->proc_ranks;
     ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
     ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
-    for (i=0; i<iga->dim; i++)
-      grid_sizes[i] = iga->axis[i]->nel;
     ierr = IGA_Partition(size,rank,iga->dim,grid_sizes,
                          proc_sizes,proc_ranks);CHKERRQ(ierr);
-    for (i=iga->dim; i<3; i++) {
+    for (i=dim; i<3; i++) {
       proc_sizes[i] = 1;
       proc_ranks[i] = 0;
     }
   }
+
   { /* element partitioning */
     PetscInt *elem_sizes = iga->elem_sizes;
     PetscInt *elem_start = iga->elem_start;
     PetscInt *elem_width = iga->elem_width;
-    for (i=0; i<iga->dim; i++)
-      elem_sizes[i] = iga->axis[i]->nel;
+    for (i=0; i<dim; i++) elem_sizes[i] = grid_sizes[i];
     ierr = IGA_Distribute(iga->dim,iga->proc_sizes,iga->proc_ranks,
                           elem_sizes,elem_width,elem_start);CHKERRQ(ierr);
-    for (i=iga->dim; i<3; i++) {
+    for (i=dim; i<3; i++) {
       elem_sizes[i] = 1;
       elem_start[i] = 0;
       elem_width[i] = 1;
     }
   }
+
+  if (!iga->collocation)
   { /* geometry partitioning */
     PetscInt *geom_sizes  = iga->geom_sizes;
     PetscInt *geom_lstart = iga->geom_lstart;
     PetscInt *geom_lwidth = iga->geom_lwidth;
     PetscInt *geom_gstart = iga->geom_gstart;
     PetscInt *geom_gwidth = iga->geom_gwidth;
-    for (i=0; i<iga->dim; i++) {
+    for (i=0; i<dim; i++) {
       PetscInt nel    = iga->elem_sizes[i];
       PetscInt efirst = iga->elem_start[i];
       PetscInt elast  = iga->elem_start[i] + iga->elem_width[i] - 1;
       geom_gstart[i] = gstart;
       geom_gwidth[i] = gend - gstart;
     }
-    for (i=iga->dim; i<3; i++) {
+    for (i=dim; i<3; i++) {
+      geom_sizes[i]  = 1;
+      geom_lstart[i] = 0;
+      geom_lwidth[i] = 1;
+      geom_gstart[i] = 0;
+      geom_gwidth[i] = 1;
+    }
+  } else
+  {  /* geometry partitioning */
+    PetscInt *geom_sizes  = iga->geom_sizes;
+    PetscInt *geom_lstart = iga->geom_lstart;
+    PetscInt *geom_lwidth = iga->geom_lwidth;
+    PetscInt *geom_gstart = iga->geom_gstart;
+    PetscInt *geom_gwidth = iga->geom_gwidth;
+    for (i=0; i<dim; i++) {
+      PetscInt   p = iga->axis[i]->p;
+      PetscInt   m = iga->axis[i]->m;
+      PetscReal *U = iga->axis[i]->U;
+      PetscInt   n = m - p - 1;
+      geom_sizes[i]  = iga->elem_sizes[i];
+      geom_lstart[i] = iga->elem_start[i];
+      geom_lwidth[i] = iga->elem_width[i];
+      {
+        PetscInt  a = geom_lstart[i];
+        PetscReal u = IGA_Greville(a,p,U);
+        PetscInt  k = IGA_FindSpan(n,p,u,U);
+        geom_gstart[i] = k - p;
+      }
+      {
+        PetscInt  a = geom_lstart[i] + geom_lwidth[i] - 1;
+        PetscReal u = IGA_Greville(a,p,U);
+        PetscInt  k = IGA_FindSpan(n,p,u,U);
+        geom_gwidth[i] = k + 1 - geom_gstart[i];
+      }
+    }
+    for (i=dim; i<3; i++) {
       geom_sizes[i]  = 1;
       geom_lstart[i] = 0;
       geom_lwidth[i] = 1;
       geom_gwidth[i] = 1;
     }
   }
+
   /* element */
   ierr = DMDestroy(&iga->elem_dm);CHKERRQ(ierr);
   /* geometry */
       node_lwidth[i] = iga->geom_lwidth[i];
       node_gstart[i] = iga->geom_gstart[i];
       node_gwidth[i] = iga->geom_gwidth[i];
+      if (rank == 0 && node_lstart[i] < 0) {
+        node_lwidth[i] += node_lstart[i];
+        node_lstart[i] = 0;
+      }
       if (rank == size-1)
         node_lwidth[i] = node_sizes[i] - node_lstart[i];
     }

src/petigabasis.c

 PetscInt IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal U[])
 {
   PetscInt low,high,span;
+  if(u <= U[p])   return p;
   if(u >= U[n+1]) return n;
-  if(u <= U[p])   return p;
   low  = p;
   high = n+1;
   span = (high+low)/2;
   element->count =  0;
   element->index = -1;
 
+  if (element->rowmap != element->mapping)
+    {ierr = PetscFree(element->rowmap);CHKERRQ(ierr);}
+  if (element->colmap != element->mapping)
+    {ierr = PetscFree(element->colmap);CHKERRQ(ierr);}
+  element->rowmap = element->colmap = 0;
   ierr = PetscFree(element->mapping);CHKERRQ(ierr);
+
   ierr = PetscFree(element->geometryX);CHKERRQ(ierr);
   ierr = PetscFree(element->rationalW);CHKERRQ(ierr);
   ierr = PetscFree(element->propertyA);CHKERRQ(ierr);
   PetscInt *start;
   PetscInt *width;
   PetscInt *sizes;
+  IGABasis *BD;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   start = iga->elem_start;
   width = iga->elem_width;
   sizes = iga->elem_sizes;
-  element->BD = iga->basis;
-  if (PetscUnlikely(element->collocation)) { /* collocation */
-    start = iga->node_lstart;
-    width = iga->node_lwidth;
-    sizes = iga->node_sizes;
-    element->BD = iga->node_basis;
-  }
+  BD = element->collocation ? iga->node_basis : iga->basis;
+
   { /* */
     PetscInt i,dim = element->dim;
     PetscInt nel=1,nen=1,nqp=1;
+    for (i=0; i<3; i++)
+      element->BD[i] = BD[i];
     for (i=0; i<dim; i++) {
       element->start[i] = start[i];
       element->width[i] = width[i];
       element->sizes[i] = sizes[i];
-      nel *= element->width[i];
-      nen *= element->BD[i]->nen;
-      nqp *= element->BD[i]->nqp;
+      nel *= width[i];
+      nen *= BD[i]->nen;
+      nqp *= BD[i]->nqp;
     }
     for (i=dim; i<3; i++) {
       element->start[i] = 0;
     element->nen   = nen;
     element->nqp   = nqp;
   }
+  { /**/
+    PetscInt nen = element->nen;
+    ierr = PetscMalloc1(nen,PetscInt,&element->mapping);CHKERRQ(ierr);
+    if (!element->collocation) {
+      element->neq = nen;
+      element->rowmap = element->mapping;
+    } else {
+      element->neq = 1;
+      ierr = PetscMalloc1(1,PetscInt,&element->rowmap);CHKERRQ(ierr);
+    }
+    element->colmap = element->mapping;
+  }
   { /* */
     PetscInt dim = element->dim;
     PetscInt nsd = element->nsd;
 
     /* */
 
-    ierr = PetscMalloc1(nen,PetscInt,&element->mapping);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*nsd,PetscReal,&element->geometryX);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen,PetscReal,&element->rationalW);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*npd,PetscScalar,&element->propertyA);CHKERRQ(ierr);
 
     ierr = PetscMemzero(element->normal,sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
   }
-  element->neq    = element->nen;
-  element->rowmap = element->mapping;
-  element->colmap = element->mapping;
-  if (PetscUnlikely(element->collocation)) { /* collocation */
-    element->neq    = 1;
-    element->rowmap = &element->index;
-    element->colmap = element->mapping;
-  }
   { /* */
     PetscInt nen = element->nen;
     PetscInt dof = element->dof;
         }
       }
     }
+    if (PetscUnlikely(element->collocation)) {
+      PetscInt iA = ID[0] - istart;
+      PetscInt jA = ID[1] - jstart;
+      PetscInt kA = ID[2] - kstart;
+      element->rowmap[0] = iA + jA*jstride + kA*kstride;
+    }
   }
   PetscFunctionReturn(0);
 }
     for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
-      if (element->index == element->colmap[a]) {
+      if (element->rowmap[0] == element->colmap[a]) {
         PetscInt    c = k % element->dof;
         PetscScalar v = element->vfix[f];
         PetscScalar u = element->ufix[f];
     for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
-      if (element->index == element->colmap[a]) {
+      if (element->rowmap[0] == element->colmap[a]) {
         PetscInt c = k % element->dof;
         PetscInt i,j;
-        if (element->index != element->colmap[a]) continue;
         for (i=0; i<M; i++) J[i*N+k] = 0.0;
         for (j=0; j<N; j++) J[c*N+j] = 0.0;
         J[c*N+k] = 1.0;