Commits

Lisandro Dalcin committed ef1adca

Add support for PCBDDC (thanks to Stefano Zampini)

  • Participants
  • Parent commits 10b99ff

Comments (0)

Files changed (7)

File include/petiga.h

 
 /* ---------------------------------------------------------------- */
 
+PETSC_EXTERN PetscErrorCode IGAPreparePCBDDC(IGA iga,PC pc);
+PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerPC(PC pc);
+PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerKSP(KSP ksp);
+PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerSNES(SNES snes);
+PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerTS(TS ts);
+
+/* ---------------------------------------------------------------- */
+
 PETSC_EXTERN PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAPointEval(IGA iga,IGAPoint point);
 

File src/makefile

 #FFLAGS_STD = ${FFLAGS_STD_F03}
 #FFLAGS = ${FFLAGS_STD} -pedantic -Wall -Wextra -fcheck=all
 
-SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigaksp.c petigasnes.c petigats.c petigats2.c petigaio.c petigagrid.c petigapart.c snesfdcolor.c tsalpha2.c
 SOURCEH  = ../include/petiga.h petigabl.h petigagrid.h petigapart.h
+SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigapc.c petigaksp.c petigasnes.c petigats.c petigats2.c petigaio.c petigagrid.c petigapart.c snesfdcolor.c tsalpha2.c
 SOURCEF1 = petigaftn.F90 petigaval.F90
 SOURCEF2 = petigabsp.f90 petigaqdr.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
 SOURCEF  = ${SOURCEF1} ${SOURCEF2}

File src/petigaksp.c

     ierr = IGAElementAssembleVec(element,B,vecB);CHKERRQ(ierr);
   }
   ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
-  
+
   ierr = PetscLogEventEnd(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
 
   ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
   ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
   ierr = KSPCreate(comm,ksp);CHKERRQ(ierr);
   ierr = PetscObjectCompose((PetscObject)*ksp,"IGA",(PetscObject)iga);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerKSP(*ksp);CHKERRQ(ierr);
   /*ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);*/
   /*ierr = KSPSetOperators(*ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);*/
   /*ierr = MatDestroy(&A);CHKERRQ(ierr);*/
   PetscFunctionReturn(0);
 }
+
+/*
+#undef  __FUNCT__
+#define __FUNCT__ "IGA_OptionsHandler_KSP"
+static PetscErrorCode IGA_OptionsHandler_KSP(PetscObject obj,void *ctx)
+{
+  KSP            ksp = (KSP)obj;
+  IGA            iga;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
+  if (PetscOptionsPublishCount != 1) PetscFunctionReturn(0);
+  ierr = PetscObjectQuery((PetscObject)ksp,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
+  if (!iga) PetscFunctionReturn(0);
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+
+  PetscFunctionReturn(0);
+}
+*/
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetOptionsHandlerKSP"
+PetscErrorCode IGASetOptionsHandlerKSP(KSP ksp)
+{
+  PC             pc;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
+  /*ierr = PetscObjectAddOptionsHandler((PetscObject)ksp,IGA_OptionsHandler_KSP,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);*/
+  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerPC(pc);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}

File src/petigapc.c

+#include "petiga.h"
+/* ---------------------------------------------------------------- */
+
+PETSC_STATIC_INLINE
+PetscInt Index(const PetscInt shape[],
+               PetscInt i,PetscInt j,PetscInt k)
+{
+  return i + j * shape[0] + k * shape[0] * shape[1];
+}
+
+PETSC_STATIC_INLINE
+PetscInt Color(const PetscInt shape[3],
+               const PetscInt start[3],
+               const PetscInt width[3],
+               PetscInt i,PetscInt j,PetscInt k)
+{
+  PetscInt r=0,g=0,b=0;
+  PetscInt L[3],R[3],C[2];
+
+  if (i<0||i>=shape[0]) return PETSC_MIN_INT;
+  if (j<0||j>=shape[1]) return PETSC_MIN_INT;
+  if (k<0||k>=shape[2]) return PETSC_MIN_INT;
+
+  L[0] = start[0]; R[0] = start[0] + width[0] - 1;
+  L[1] = start[1]; R[1] = start[1] + width[1] - 1;
+  L[2] = start[2]; R[2] = start[2] + width[2] - 1;
+  if (i<L[0]) r = i - L[0]; if (i>R[0]) r = i - R[0];
+  if (j<L[1]) g = j - L[1]; if (j>R[1]) g = j - R[1];
+  if (k<L[2]) b = k - L[2]; if (k>R[2]) b = k - R[2];
+
+  C[0] = shape[0] - width[0] + 1;
+  C[1] = shape[1] - width[1] + 1;
+  return Index(C,r,g,b);
+}
+
+static const
+PetscInt STENCIL[7][3] = {{ 0,  0, -1},
+                          { 0, -1,  0},
+                          {-1,  0,  0},
+                          { 0,  0,  0},
+                          { 0,  0, +1},
+                          { 0, +1,  0},
+                          {+1,  0,  0}};
+
+PETSC_STATIC_INLINE
+#undef  __FUNCT__
+#define __FUNCT__ "ComputeBDDCGraph"
+PetscErrorCode ComputeBDDCGraph(PetscInt bs,const PetscInt shape[3],
+                                const PetscInt start[3],const PetscInt width[3],
+                                PetscInt *_nvtx,PetscInt *_xadj[],PetscInt *_adjy[])
+{
+  PetscInt       c,i,j,k,s,v,pos;
+  PetscInt       nvtx=0,*xadj;
+  PetscInt       nadj=0,*adjy;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidIntPointer(_nvtx,5);
+  PetscValidPointer(_xadj,6);
+  PetscValidPointer(_adjy,7);
+
+  /* Compute the number of vertices and adjacencies */
+  nvtx = shape[0]*shape[1]*shape[2];
+  for (k=0; k<shape[2]; k++)
+    for (j=0; j<shape[1]; j++)
+      for (i=0; i<shape[0]; i++) {
+        PetscInt color = Color(shape,start,width,i,j,k);
+        for (s=0; s<7; s++) {
+          PetscInt ii = i + STENCIL[s][0];
+          PetscInt jj = j + STENCIL[s][1];
+          PetscInt kk = k + STENCIL[s][2];
+          PetscInt cc = Color(shape,start,width,ii,jj,kk);
+          if (cc == color) nadj++;
+        }
+      }
+
+  /* Allocate arrays to store the adjacency graph */
+  nvtx *= bs; nadj *= bs; /* adjust for block size */
+  ierr = PetscMalloc((nvtx+1)*sizeof(PetscInt),&xadj);CHKERRQ(ierr);
+  ierr = PetscMalloc(nadj*sizeof(PetscInt),&adjy);CHKERRQ(ierr);
+
+  /* Fill the adjacency graph */
+  pos = 0; xadj[pos++] = 0;
+  for (k=0; k<shape[2]; k++)
+    for (j=0; j<shape[1]; j++)
+      for (i=0; i<shape[0]; i++) {
+        /* Compute the color of this vertex */
+        PetscInt color = Color(shape,start,width,i,j,k);
+        /* Compute the list of neighbor vertices
+           having the same color */
+        PetscInt nv = 0, vertices[7];
+        for (s=0; s<7; s++) { /* loop over neighbors */
+          PetscInt ii = i + STENCIL[s][0];
+          PetscInt jj = j + STENCIL[s][1];
+          PetscInt kk = k + STENCIL[s][2];
+          PetscInt cc = Color(shape,start,width,ii,jj,kk);
+          if (cc == color)
+            vertices[nv++] = Index(shape,ii,jj,kk);
+        }
+        for (c=0; c<bs; c++) {
+          xadj[pos] = xadj[pos-1];
+          for (v=0; v<nv; v++)
+            adjy[xadj[pos]++] = c + bs*vertices[v];
+          pos++;
+        }
+      }
+
+  *_nvtx = nvtx;
+  *_xadj = xadj;
+  *_adjy = adjy;
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPreparePCBDDC"
+PetscErrorCode IGAPreparePCBDDC(IGA iga,PC pc)
+{
+  Mat            mat;
+  void           (*f)(void);
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(pc,PC_CLASSID,2);
+  IGACheckSetUp(iga,1);
+
+  ierr = PCGetOperators(pc,PETSC_NULL,&mat,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscObjectQueryFunction((PetscObject)mat,"MatISGetLocalMat_C",&f);CHKERRQ(ierr);
+  if (!f) PetscFunctionReturn(0);
+  ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr);
+  if (!f) PetscFunctionReturn(0);
+
+#if defined(PETSC_HAVE_PCBDDC)  /* XXX */
+  {
+    PetscInt i,dim,dof;
+    PetscInt shape[3] = {1,1,1};
+    PetscInt start[3] = {0,0,0};
+    PetscInt width[3] = {1,1,1};
+    PetscInt nvtx=0,*xadj=0,*adjy=0;
+    ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+    ierr = IGAGetDof(iga,&dof);CHKERRQ(ierr);
+    for (i=0; i<dim; i++) {
+      shape[i] = iga->node_gwidth[i];
+      start[i] = iga->node_lstart[i] - iga->node_gstart[i];
+      width[i] = iga->node_lwidth[i];
+    }
+    ierr = ComputeBDDCGraph(dof,shape,start,width,&nvtx,&xadj,&adjy);CHKERRQ(ierr);
+    ierr = PCBDDCSetLocalAdjacencyGraph(pc,nvtx,xadj,adjy,PETSC_OWN_POINTER);CHKERRQ(ierr);
+  }
+#endif
+
+  PetscFunctionReturn(0);
+}
+
+/* ---------------------------------------------------------------- */
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGA_OptionsHandler_PC"
+static PetscErrorCode IGA_OptionsHandler_PC(PetscObject obj,void *ctx)
+{
+  PC             pc = (PC)obj;
+  Mat            mat;
+  IGA            iga;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  if (PetscOptionsPublishCount != 1) PetscFunctionReturn(0);
+  ierr = PCGetOperators(pc,PETSC_NULL,&mat,PETSC_NULL);CHKERRQ(ierr);
+  if (!mat) PetscFunctionReturn(0);
+  PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
+  ierr = PetscObjectQuery((PetscObject)mat,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
+  if (!iga) PetscFunctionReturn(0);
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  /* */
+  ierr = IGAPreparePCBDDC(iga,pc);CHKERRQ(ierr);
+  /* */
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetOptionsHandlerPC"
+PetscErrorCode IGASetOptionsHandlerPC(PC pc)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
+  ierr = PetscObjectAddOptionsHandler((PetscObject)pc,IGA_OptionsHandler_PC,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/* ---------------------------------------------------------------- */

File src/petigasnes.c

   ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
   ierr = SNESCreate(comm,snes);CHKERRQ(ierr);
   ierr = PetscObjectCompose((PetscObject)*snes,"IGA",(PetscObject)iga);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerSNES(*snes);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&F);CHKERRQ(ierr);
   ierr = SNESSetFunction(*snes,F,IGASNESFormFunction,iga);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
+
+/*
+#undef  __FUNCT__
+#define __FUNCT__ "IGA_OptionsHandler_SNES"
+static PetscErrorCode IGA_OptionsHandler_SNES(PetscObject obj,void *ctx)
+{
+  SNES            snes = (SNES)obj;
+  IGA            iga;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
+  if (PetscOptionsPublishCount != 1) PetscFunctionReturn(0);
+  ierr = PetscObjectQuery((PetscObject)snes,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
+  if (!iga) PetscFunctionReturn(0);
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+
+  PetscFunctionReturn(0);
+}
+*/
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetOptionsHandlerSNES"
+PetscErrorCode IGASetOptionsHandlerSNES(SNES snes)
+{
+  KSP            ksp;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
+  /*ierr = PetscObjectAddOptionsHandler((PetscObject)snes,IGA_OptionsHandler_SNES,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);*/
+  ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerKSP(ksp);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}

File src/petigats.c

   ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
   ierr = TSCreate(comm,ts);CHKERRQ(ierr);
   ierr = PetscObjectCompose((PetscObject)*ts,"IGA",(PetscObject)iga);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerTS(*ts);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
   ierr = TSSetSolution(*ts,U);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
+
+/*
+#undef  __FUNCT__
+#define __FUNCT__ "IGA_OptionsHandler_TS"
+static PetscErrorCode IGA_OptionsHandler_TS(PetscObject obj,void *ctx)
+{
+  TS             ts = (TS)obj;
+  IGA            iga;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(ts,TS_CLASSID,1);
+  if (PetscOptionsPublishCount != 1) PetscFunctionReturn(0);
+  ierr = PetscObjectQuery((PetscObject)ts,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
+  if (!iga) PetscFunctionReturn(0);
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+
+  PetscFunctionReturn(0);
+}
+*/
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetOptionsHandlerTS"
+PetscErrorCode IGASetOptionsHandlerTS(TS ts)
+{
+  SNES           snes;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(ts,TS_CLASSID,1);
+  /*ierr = PetscObjectAddOptionsHandler((PetscObject)ts,IGA_OptionsHandler_TS,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);*/
+  ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerSNES(snes);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}

File src/petigats2.c

   ierr = TSCreate(comm,ts);CHKERRQ(ierr);
   ierr = TSSetType(*ts,TSALPHA2);CHKERRQ(ierr);
   ierr = PetscObjectCompose((PetscObject)*ts,"IGA",(PetscObject)iga);CHKERRQ(ierr);
+  ierr = IGASetOptionsHandlerTS(*ts);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&V);CHKERRQ(ierr);