Commits

Stefano Zampini  committed 892d802

PCBDDC: More robust method to compute local boundary nodes

Using counter vector to infer local boundary nodes is wrong in case a
subdomain has local periodicity.

  • Participants
  • Parent commits 93fb5fd

Comments (0)

Files changed (1)

File src/ksp/pc/impls/is/pcis.c

 {
   PC_IS          *pcis  = (PC_IS*)(pc->data);
   Mat_IS         *matis = (Mat_IS*)pc->mat->data;
-  PetscInt       i;
   PetscErrorCode ierr;
   PetscBool      flg;
   Vec            counter;
 
   pcis->pure_neumann = matis->pure_neumann;
 
-  /*
-    Creating the local vector vec1_N, containing the inverse of the number
-    of subdomains to which each local node (either owned or ghost)
-    pertains. To accomplish that, we scatter local vectors of 1's to
-    a global vector (adding the values); scatter the result back to
-    local vectors and finally invert the result.
-  */
-  ierr = VecDuplicate(matis->x,&pcis->vec1_N);CHKERRQ(ierr);
-  ierr = MatGetVecs(pc->pmat,&counter,0);CHKERRQ(ierr); /* temporary auxiliar vector */
-  ierr = VecSet(counter,0.0);CHKERRQ(ierr);
-  ierr = VecSet(pcis->vec1_N,1.0);CHKERRQ(ierr);
-  ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterBegin(matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd  (matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  /*
-    Creating local and global index sets for interior and
-    inteface nodes. Notice that interior nodes have D[i]==1.0.
-  */
+  /* get info on mapping */
+  ierr = ISLocalToGlobalMappingGetInfo(((Mat_IS*)(pc->mat->data))->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));CHKERRQ(ierr);
+  pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_TRUE;
+
+  /* Creating local and global index sets for interior and inteface nodes. */
   {
     PetscInt    n_I;
     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
-    PetscScalar *array;
+    PetscInt    *array;
+    PetscInt    i,j;
+
     /* Identifying interior and interface nodes, in local numbering */
-    ierr = VecGetSize(pcis->vec1_N,&pcis->n);CHKERRQ(ierr);
-    ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
+    ierr = VecGetSize(matis->x,&pcis->n);CHKERRQ(ierr);
+    ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&array);CHKERRQ(ierr);
+    ierr = PetscMemzero(array,pcis->n*sizeof(PetscInt));CHKERRQ(ierr);
+    for (i=1;i<pcis->n_neigh;i++)
+      for (j=0;j<pcis->n_shared[i];j++)
+          array[pcis->shared[i][j]] += 1;
+ 
     ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&idx_I_local);CHKERRQ(ierr);
     ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&idx_B_local);CHKERRQ(ierr);
     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
-      if (array[i] == 1.0) {
+      if (!array[i]) {
         idx_I_local[n_I] = i;
         n_I++;
       } else {
     /* Freeing memory and restoring arrays */
     ierr = PetscFree(idx_B_local);CHKERRQ(ierr);
     ierr = PetscFree(idx_I_local);CHKERRQ(ierr);
-    ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
+    ierr = PetscFree(array);CHKERRQ(ierr);
   }
 
   /*
   /*
     Creating work vectors and arrays
   */
-  /* pcis->vec1_N has already been created */
+  ierr = VecDuplicate(matis->x,&pcis->vec1_N);CHKERRQ(ierr);
   ierr = VecDuplicate(pcis->vec1_N,&pcis->vec2_N);CHKERRQ(ierr);
   ierr = VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);CHKERRQ(ierr);
   ierr = VecDuplicate(pcis->vec1_D,&pcis->vec2_D);CHKERRQ(ierr);
     }
   }
   ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
+  ierr = MatGetVecs(pc->pmat,&counter,0);CHKERRQ(ierr); /* temporary auxiliar vector */
   ierr = VecSet(counter,0.0);CHKERRQ(ierr);
   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
   ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
   ierr = VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   ierr = VecScatterEnd  (pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
+  ierr = VecDestroy(&counter);CHKERRQ(ierr);
 
   /* See historical note 01, at the bottom of this file. */
 
     ierr = KSPSetUp(pcis->ksp_N);CHKERRQ(ierr);
   }
 
-  ierr = ISLocalToGlobalMappingGetInfo(((Mat_IS*)(pc->mat->data))->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));CHKERRQ(ierr);
-
-  pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_TRUE;
-
-  ierr = VecDestroy(&counter);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }