Jed Brown avatar Jed Brown committed 792fabe

TS ex18.c: new unstructured FEM heat equation

Comments (0)

Files changed (2)

src/ts/examples/tutorials/ex18.c

+static char help[] = "Unstructured heat equation solver on simplicial meshes.\n";
+/*
+  The mesh is read in from an ExodusII file generated by Cubit.
+*/
+#include <petscts.h>
+#include <petscdmplex.h>
+#if defined(PETSC_HAVE_EXODUSII)
+#  include <exodusII.h>
+#else
+#  error This example requires ExodusII support. Reconfigure using --download-exodusii
+#endif
+
+/* Generate element tabulation using:
+ *
+ *  python2 ./bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 1 1 laplacian src/ts/examples/tutorials/ex18.h
+ *
+ * This example needs a consistent mass matrix so PetscGenerateFEMQuadrature.py needs to be modified to raise the
+ * quadrature order by one.
+ */
+
+#include "ex18.h"
+
+#define DIM SPATIAL_DIM_0
+#define ALEN(a) (sizeof(a)/sizeof((a)[0]))
+
+/* 'User' implements a discretization of a continuous model. */
+typedef struct _n_User *User;
+
+struct _n_User {
+  PetscInt vtkInterval;
+};
+
+PETSC_STATIC_INLINE PetscScalar DotDIM(const PetscScalar *x,const PetscScalar *y)
+{
+  PetscInt    i;
+  PetscScalar prod=0.0;
+
+  for (i=0; i<DIM; i++) prod += x[i]*y[i];
+  return prod;
+}
+PETSC_STATIC_INLINE PetscReal NormDIM(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(DotDIM(x,x))); }
+PETSC_STATIC_INLINE void axDIM(const PetscScalar a,PetscScalar *x)
+{
+  PetscInt i;
+  for (i=0; i<DIM; i++) x[i] *= a;
+}
+PETSC_STATIC_INLINE void waxDIM(const PetscScalar a,const PetscScalar *x,PetscScalar *w)
+{
+  PetscInt i;
+  for (i=0; i<DIM; i++) w[i] = x[i]*a;
+}
+PETSC_STATIC_INLINE void NormalSplitDIM(const PetscReal *n,const PetscScalar *x,PetscScalar *xn,PetscScalar *xt)
+{                               /* Split x into normal and tangential components */
+  PetscInt    i;
+  PetscScalar c;
+  c = DotDIM(x,n)/DotDIM(n,n);
+  for (i=0; i<DIM; i++) {
+    xn[i] = c*n[i];
+    xt[i] = x[i]-xn[i];
+  }
+}
+
+PETSC_STATIC_INLINE void Mult(PetscInt m,PetscInt n,const PetscReal *A,const PetscReal *x,PetscReal *y) {
+  PetscInt i,j;
+  for (i=0; i<m; i++) {
+    y[i] = 0;
+    for (j=0; j<n; j++) y[i] += A[i*n+j]*x[j];
+  }
+}
+
+PETSC_STATIC_INLINE void MultTrans(PetscInt m,PetscInt n,const PetscReal *A,const PetscReal *x,PetscReal *y) {
+  PetscInt i,j;
+  for (j=0; j<n; j++) {
+    y[j] = 0;
+    for (i=0; i<m; i++) {
+      y[j] += A[i*n+j]*x[i];
+    }
+  }
+}
+
+PETSC_STATIC_INLINE void MultTransAdd(PetscInt m,PetscInt n,const PetscReal *A,const PetscReal *x,PetscReal *y) {
+  PetscInt i,j;
+  for (j=0; j<n; j++) {
+    for (i=0; i<m; i++) {
+      y[j] += A[i*n+j]*x[i];
+    }
+  }
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "UserSetupElem"
+/* Map to physical element */
+PETSC_STATIC_INLINE PetscErrorCode UserSetupElem(DM dm,PetscInt cell,const PetscReal **Weight,const PetscReal **Basis,const PetscReal **Deriv)
+{
+  static PetscReal _weights[NUM_QUADRATURE_POINTS_0];
+  static PetscReal _deriv[NUM_QUADRATURE_POINTS_0][DIM][NUM_BASIS_FUNCTIONS_0];
+  PetscReal        J[DIM][DIM],Jinv[DIM][DIM],Jdet,v0[DIM];
+  PetscInt         q,p,d,d2;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexComputeCellGeometry(dm,cell,v0,&J[0][0],&Jinv[0][0],&Jdet);CHKERRQ(ierr);
+  if (Jdet < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Inverted element %D Jdet %G",cell,Jdet);
+  for (q=0; q<NUM_QUADRATURE_POINTS_0; q++) {
+    _weights[q] = weights_0[q] * Jdet;
+    for (p=0; p<NUM_BASIS_FUNCTIONS_0; p++) {
+      for (d=0; d<DIM; d++) {
+        PetscReal tmp = 0;
+        for (d2=0; d2<DIM; d2++) {
+          /* Jinv[i][j] is derivative of reference coordinate i with respect to physical coordinate j */
+          tmp += Jinv[d2][d] * BasisDerivatives_0[(q*NUM_BASIS_FUNCTIONS_0+p)*DIM + d2];
+        }
+        _deriv[q][d][p] = tmp;
+      }
+    }
+  }
+  *Weight = _weights;
+  *Basis  = Basis_0;
+  *Deriv  = &_deriv[0][0][0];
+
+   /* Squelch unused compiler warning */
+  (void)points_0;
+  (void)numDof_0;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "UserGetCoefficients"
+static PetscErrorCode UserGetCoefficients(User user,PetscInt cellset,PetscReal *Kheat,PetscReal *HeatSource)
+{
+
+  PetscFunctionBegin;
+  if (cellset == 10) {          /* Interior */
+    *Kheat = 5.;
+    *HeatSource = 1.;
+  } else {
+    *Kheat = 1.;
+    *HeatSource = 0.;
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "UserSetupSection"
+static PetscErrorCode UserSetupSection(User user,DM dm)
+{
+  PetscErrorCode ierr;
+  PetscInt numComp[1] = {1},numDof[4] = {0,0,0,0};
+  PetscSection section;
+
+  PetscFunctionBegin;
+  numDof[0] = 1;
+  ierr = DMPlexCreateSection(dm,DIM,1,numComp,numDof,0,NULL,NULL,&section);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dm,section);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
+
+  {
+    PetscInt nlabels,i;
+    ierr = DMPlexGetNumLabels(dm,&nlabels);CHKERRQ(ierr);
+    for (i=0; i<nlabels; i++) {
+      const char *name;
+      ierr = DMPlexGetLabelName(dm,i,&name);CHKERRQ(ierr);
+      ierr = PetscPrintf(PetscObjectComm((PetscObject)dm),"Label % 2D: %s\n",i,name);CHKERRQ(ierr);
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetInitialCondition"
+static PetscErrorCode SetInitialCondition(DM dm,Vec U,User user)
+{
+  DM                cdm;
+  const PetscScalar *x;
+  PetscScalar       *u;
+  PetscInt          v,vStart,vEnd;
+  Vec               X;
+  PetscErrorCode    ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);CHKERRQ(ierr);
+  ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
+  ierr = DMGetCoordinates(dm,&X);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(U,&u);CHKERRQ(ierr);
+
+  /* Vertex loop, setting solution state */
+  for (v=vStart; v<vEnd; v++) {
+    const PetscScalar *vx;
+    PetscScalar *vu;
+    ierr = DMPlexPointGlobalRead(cdm,v,x,&vx);CHKERRQ(ierr);
+    ierr = DMPlexPointGlobalRef(dm,v,u,&vu);CHKERRQ(ierr);
+    if (!vu) continue;
+    vu[0] = DotDIM(vx,vx);
+  }
+  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "IFunction"
+static PetscErrorCode IFunction(TS ts,PetscReal time,Vec U,Vec Udot,Vec F,void *ctx)
+{
+  User              user = (User)ctx;
+  DM                dm;
+  Vec               locU,locUdot,locF;
+  PetscErrorCode    ierr;
+  PetscSection      section;
+  PetscInt          c,cStart,cEnd;
+
+  PetscFunctionBeginUser;
+  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(dm,&locU);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(dm,&locUdot);CHKERRQ(ierr);
+
+  ierr = DMGlobalToLocalBegin(dm,U,INSERT_VALUES,locU);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalEnd(dm,U,INSERT_VALUES,locU);CHKERRQ(ierr);
+
+  ierr = DMGlobalToLocalBegin(dm,Udot,INSERT_VALUES,locUdot);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalEnd(dm,Udot,INSERT_VALUES,locUdot);CHKERRQ(ierr);
+
+  ierr = DMGetLocalVector(dm,&locF);CHKERRQ(ierr);
+  ierr = VecZeroEntries(locF);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
+  ierr = DMGetDefaultSection(dm,&section);CHKERRQ(ierr);
+  for (c=cStart; c<cEnd; c++) {
+    PetscScalar     *cu,*cudot;
+    PetscInt        csize,cellset,q,d;
+    const PetscReal *Weight,*Basis,*Deriv;
+    PetscReal       Kheat,HeatSource;
+    PetscScalar     qu[NUM_QUADRATURE_POINTS_0],qudot[NUM_QUADRATURE_POINTS_0],qdu[NUM_QUADRATURE_POINTS_0][DIM];
+    PetscScalar     qf[NUM_QUADRATURE_POINTS_0],qdf[NUM_QUADRATURE_POINTS_0][DIM],cf[NUM_BASIS_FUNCTIONS_0];
+
+    ierr = UserSetupElem(dm,c,&Weight,&Basis,&Deriv);CHKERRQ(ierr);
+    ierr = DMPlexGetLabelValue(dm,"Cell Sets",c,&cellset);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosure(dm,section,locU,c,&csize,&cu);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosure(dm,section,locUdot,c,&csize,&cudot);CHKERRQ(ierr);
+    if (csize != NUM_BASIS_FUNCTIONS_0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for dim %D csize %D",DIM,csize);CHKERRQ(ierr);
+
+    Mult(NUM_QUADRATURE_POINTS_0,NUM_BASIS_FUNCTIONS_0,Basis,cu,qu);
+    Mult(NUM_QUADRATURE_POINTS_0,NUM_BASIS_FUNCTIONS_0,Basis,cudot,qudot);
+    Mult(NUM_QUADRATURE_POINTS_0*DIM,NUM_BASIS_FUNCTIONS_0,Deriv,cu,&qdu[0][0]);
+
+    ierr = UserGetCoefficients(user,cellset,&Kheat,&HeatSource);CHKERRQ(ierr);
+    for (q=0; q<NUM_QUADRATURE_POINTS_0; q++) {
+      qf[q] = Weight[q] * (qudot[q] - HeatSource);
+      for (d=0; d<DIM; d++) qdf[q][d] = Weight[q] * Kheat * qdu[q][d];
+    }
+    MultTrans(NUM_QUADRATURE_POINTS_0,NUM_BASIS_FUNCTIONS_0,Basis,qf,cf);
+    MultTransAdd(NUM_QUADRATURE_POINTS_0*DIM,NUM_BASIS_FUNCTIONS_0,Deriv,&qdf[0][0],cf);
+
+    ierr = DMPlexVecRestoreClosure(dm,section,locU,c,&csize,&cu);CHKERRQ(ierr);
+    ierr = DMPlexVecRestoreClosure(dm,section,locUdot,c,&csize,&cudot);CHKERRQ(ierr);
+    ierr = DMPlexVecSetClosure(dm,section,locF,c,cf,ADD_VALUES);CHKERRQ(ierr);
+  }
+
+  ierr = VecZeroEntries(F);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalBegin(dm,locF,ADD_VALUES,F);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalEnd(dm,locF,ADD_VALUES,F);CHKERRQ(ierr);
+
+  ierr = DMRestoreLocalVector(dm,&locU);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(dm,&locUdot);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(dm,&locF);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "OutputVTK"
+static PetscErrorCode OutputVTK(DM dm,const char *filename,PetscViewer *viewer)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm),viewer);CHKERRQ(ierr);
+  ierr = PetscViewerSetType(*viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
+  ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MonitorVTK"
+static PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec U,void *ctx)
+{
+  User           user = (User)ctx;
+  DM             dm;
+  PetscViewer    viewer;
+  char           filename[PETSC_MAX_PATH_LEN];
+  PetscReal      unorm;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = PetscObjectSetName((PetscObject)U,"Temperature");CHKERRQ(ierr);
+  ierr = VecGetDM(U,&dm);CHKERRQ(ierr);
+  ierr = VecNorm(U,NORM_INFINITY,&unorm);CHKERRQ(ierr);
+
+  if (user->vtkInterval < 1) PetscFunctionReturn(0);
+  if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
+    if (stepnum == -1) {        /* Final time is not multiple of normal time interval,write it anyway */
+      ierr = TSGetTimeStepNumber(ts,&stepnum);CHKERRQ(ierr);
+    }
+    ierr = PetscSNPrintf(filename,sizeof filename,"ex18-%03D.vtu",stepnum);CHKERRQ(ierr);
+    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
+    ierr = VecView(U,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc,char **argv)
+{
+  MPI_Comm          comm;
+  User              user;
+  DM                dm,dmDist;
+  PetscReal         ftime,dt;
+  PetscInt          dim,nsteps;
+  int               CPU_word_size = 0,IO_word_size = 0,exoid;
+  float             version;
+  TS                ts;
+  TSConvergedReason reason;
+  Vec               U;
+  PetscMPIInt       rank;
+  char              filename[PETSC_MAX_PATH_LEN] = "sevenside.exo";
+  PetscErrorCode    ierr;
+
+  ierr = PetscInitialize(&argc,&argv,(char*) 0,help);CHKERRQ(ierr);
+  comm = PETSC_COMM_WORLD;
+  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
+
+  ierr = PetscNew(struct _n_User,&user);CHKERRQ(ierr);
+
+  ierr = PetscOptionsBegin(comm,NULL,"Unstructured Finite Element Heat Options","");CHKERRQ(ierr);
+  {
+    ierr = PetscOptionsString("-f","Exodus.II filename to read","",filename,filename,sizeof(filename),NULL);CHKERRQ(ierr);
+    user->vtkInterval = 1;
+    ierr = PetscOptionsInt("-ufe_vtk_interval","VTK output interval (0 to disable)","",user->vtkInterval,&user->vtkInterval,NULL);CHKERRQ(ierr);
+  }
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  if (!rank) {
+    exoid = ex_open(filename,EX_READ,&CPU_word_size,&IO_word_size,&version);
+    if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ex_open(\"%s\",...) did not return a valid file ID",filename);
+  } else exoid = -1;                 /* Not used */
+  ierr = DMPlexCreateExodus(comm,exoid,PETSC_FALSE,&dm);CHKERRQ(ierr);
+  if (!rank) {ierr = ex_close(exoid);CHKERRQ(ierr);}
+  /* Distribute mesh */
+  ierr = DMPlexDistribute(dm,"chaco",0,&dmDist);CHKERRQ(ierr);
+  if (dmDist) {
+    ierr = DMDestroy(&dm);CHKERRQ(ierr);
+    dm   = dmDist;
+  }
+  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
+
+  ierr = DMPlexGetDimension(dm,&dim);CHKERRQ(ierr);
+  ierr = DMPlexSetPreallocationCenterDimension(dm,dim);CHKERRQ(ierr);
+  ierr = UserSetupSection(user,dm);CHKERRQ(ierr);
+
+  ierr = DMCreateGlobalVector(dm,&U);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject)U,"Temperature");CHKERRQ(ierr);
+  ierr = SetInitialCondition(dm,U,user);CHKERRQ(ierr);
+
+  ierr = TSCreate(comm,&ts);CHKERRQ(ierr);
+  ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
+  ierr = TSSetDM(ts,dm);CHKERRQ(ierr);
+  ierr = TSMonitorSet(ts,MonitorVTK,user,NULL);CHKERRQ(ierr);
+  ierr = TSSetIFunction(ts,NULL,IFunction,user);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,1000,1.0);CHKERRQ(ierr);
+  dt = 0.1;
+  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  ierr = TSSolve(ts,U);CHKERRQ(ierr);
+  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
+  ierr = TSGetTimeStepNumber(ts,&nsteps);CHKERRQ(ierr);
+  ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,nsteps);CHKERRQ(ierr);
+
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = PetscFree(user);CHKERRQ(ierr);
+  ierr = VecDestroy(&U);CHKERRQ(ierr);
+  ierr = DMDestroy(&dm);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+  return(0);
+}

src/ts/examples/tutorials/makefile

 FPPFLAGS        =
 LOCDIR          = src/ts/examples/tutorials/
 EXAMPLESC       = ex1.c ex2.c ex3.c ex4.c ex5.c ex6.c ex7.c ex8.c \
-                ex9.c ex10.c ex12.c ex13.c ex14.c ex15.c ex16.c ex17.c \
+                ex9.c ex10.c ex12.c ex13.c ex14.c ex15.c ex16.c ex17.c ex18.c \
                 ex19.c ex20.c ex21.c ex22.c ex23.c ex24.c ex25.c ex26.c ex28.c
 EXAMPLESF       = ex1f.F ex2f.F ex22f.F # ex22f_mf.F90
 EXAMPLESFH      = ex2f.h
 	-${CLINKER} -o ex17 ex17.o  ${PETSC_TS_LIB}
 	${RM} ex17.o
 
+ex18: ex18.o  chkopts
+	-${CLINKER} -o ex18 ex18.o  ${PETSC_TS_LIB}
+	${RM} ex18.o
+
 ex19: ex19.o  chkopts
 	-${CLINKER} -o ex19 ex19.o  ${PETSC_TS_LIB}
 	${RM} ex19.o
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.