Commits

Matt Knepley committed 6535129

SNES ex12: Fixed output and preparing for Neumann problem

Hg-commit: 80af08fbb63450ee85481f58bc722b97bf4164db

Comments (0)

Files changed (4)

src/snes/examples/tutorials/ex12.c

-static char help[] = "Poisson Problem in 2d and 3d with simplicial and hexahedral finite elements.\n\
+static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements.\n\
 We solve the Poisson problem in a rectangular\n\
 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
 
 
 /*------------------------------------------------------------------------------
   This code can be generated using 'bin/pythonscripts/PetscGenerateFEMQuadrature.py dim order dim 1 laplacian src/snes/examples/tutorials/ex12.h'
-                                   'bin/pythonscripts/PetscGenerateFEMQuadratureTensorProduct.py dim order dim 1 laplacian src/snes/examples/tutorials/ex12.h'
  -----------------------------------------------------------------------------*/
 #include "ex12.h"
 
 }
 
 /*
-  In 2D we use exact solution:
+  In 2D for Dirichlet conditions, we use exact solution:
 
     u = x^2 + y^2
     f = 4
   so that
 
     -\Delta u + f = -4 + 4 = 0
+
+  For Neumann conditions, we have
+
+    \nabla u \cdot -\hat y |_{y=0} = -(2y)|_{y=0} = 0 (bottom)
+    \nabla u \cdot  \hat y |_{y=1} =  (2y)|_{y=1} = 2 (top)
+    \nabla u \cdot -\hat x |_{x=0} = -(2x)|_{x=0} = 0 (left)
+    \nabla u \cdot  \hat x |_{x=1} =  (2x)|_{x=1} = 2 (right)
 */
 PetscScalar quadratic_u_2d(const PetscReal x[])
 {
   IS             bcPoints[1]         = {NULL};
   PetscInt       numDof[NUM_FIELDS*(SPATIAL_DIM_0+1)];
   PetscInt       f, d;
+  PetscBool      has;
   PetscErrorCode ierr;
 
   PetscFunctionBeginUser;
       if ((numDof[f*(dim+1)+d] > 0) && !user->interpolate) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
     }
   }
+  ierr = DMPlexHasLabel(dm, "marker", &has);CHKERRQ(ierr);
+  if (!has) {
+    DMLabel label;
+
+    ierr = DMPlexCreateLabel(dm, "marker");CHKERRQ(ierr);
+    ierr = DMPlexGetLabel(dm, "marker", &label);CHKERRQ(ierr);
+    ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr);
+    if (user->bcType == DIRICHLET) {
+      ierr  = DMPlexLabelComplete(dm, label);CHKERRQ(ierr);
+    }
+  }
   if (user->bcType == DIRICHLET) {
     numBC = 1;
     ierr  = DMPlexGetStratumIS(dm, "marker", 1, &bcPoints[0]);CHKERRQ(ierr);
   SNES           snes;                 /* nonlinear solver */
   Vec            u,r;                  /* solution, residual vectors */
   Mat            A,J;                  /* Jacobian matrix */
-  MatNullSpace   nullSpace;            /* May be necessary for pressure */
+  MatNullSpace   nullSpace;            /* May be necessary for Neumann conditions */
   AppCtx         user;                 /* user-defined work context */
   JacActionCtx   userJ;                /* context for Jacobian MF action */
   PetscInt       its;                  /* iterations for convergence */
   } else {
     A = J;
   }
+#if 0
+  if (user.bcType == NEUMANN) {
+    ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject) user.dm), PETSC_TRUE, 0, NULL, &nullSpace);CHKERRQ(ierr);
+    ierr = MatSetNullSpace(J, nullSpace);CHKERRQ(ierr);
+    if (A != J) {
+      ierr = MatSetNullSpace(A, nullSpace);CHKERRQ(ierr);
+    }
+  }
+#endif
 
   ierr = DMSNESSetFunctionLocal(user.dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);CHKERRQ(ierr);
   ierr = DMSNESSetJacobianLocal(user.dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*))DMPlexComputeJacobianFEM,&user);CHKERRQ(ierr);

src/snes/examples/tutorials/output/ex12_1.out

 Mesh in 3 dimensions:
   0-cells: 8
   3-cells: 6
+Labels:
+  marker: 1 strata of sizes (8)
+  depth: 2 strata of sizes (8, 6)
 Local function
 Vector Object: 1 MPI processes
   type: seq

src/snes/examples/tutorials/output/ex12_2.out

 Mesh in 3 dimensions:
   0-cells: 27
   3-cells: 48
+Labels:
+  marker: 1 strata of sizes (26)
+  depth: 2 strata of sizes (27, 48)
 Local function
 Vector Object: 1 MPI processes
   type: seq

src/snes/examples/tutorials/output/ex12_3.out

-Residual:
-Vector Object: 1 MPI processes
-  type: seq
-0
-0
-0
-0
-0.333333
-0
--0.166667
--2.16667
-Vector Object: 1 MPI processes
-  type: seq
-0
-0.333333
-0
-0
-0
-0
--2.16667
--4.16667
-  0 SNES Function norm 5.217491947500e+00 
-    0 KSP Residual norm 1.289972006587e+00 
-    1 KSP Residual norm 4.648023243078e-01 
-    2 KSP Residual norm 8.105143815092e-02 
-    3 KSP Residual norm 4.765549829478e-03 
-    4 KSP Residual norm 1.806187445385e-04 
-    5 KSP Residual norm 2.470053047512e-16 
-Residual:
-Vector Object: 1 MPI processes
-  type: seq
-0
-0
-0
-0
-0.5
-0
-0
-4.44089e-16
-Vector Object: 1 MPI processes
-  type: seq
-0
--0.5
-0
-0
-0
-0
-8.88178e-16
-1.33227e-15
-  1 SNES Function norm 1.676400004429e-15 
-Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE
-Number of SNES iterations = 1
-L_2 Error: 0.0694444
-Solution
-Vector Object: 2 MPI processes
-  type: mpi
-Process [0]
-0.5
-0.166667
-0.666667
-Process [1]
-0.666667
-1.16667