1: static const char help[] ="Minimum surface problem in 2D.\n\
  2: Uses 2-dimensional distributed arrays.\n\
  3: \n\
  4:   Solves the linear systems via multilevel methods \n\
  5: \n\n";

  7: /*T
  8:    Concepts: SNES^solving a system of nonlinear equations
  9:    Concepts: DMDA^using distributed arrays
 10:    Concepts: multigrid;
 11:    Processors: n
 12: T*/

 14: /*  
 15:   
 16:     This example models the partial differential equation 
 17:    
 18:          - Div((1 + ||GRAD T||^2)^(1/2) (GRAD T)) = 0.
 19:        
 20:     
 21:     in the unit square, which is uniformly discretized in each of x and 
 22:     y in this simple encoding.  The degrees of freedom are vertex centered
 23:  
 24:     A finite difference approximation with the usual 5-point stencil 
 25:     is used to discretize the boundary value problem to obtain a 
 26:     nonlinear system of equations. 
 27:  
 28: */

 30: #include <petscsnes.h>
 31: #include <petscdmda.h>
 32: #include <petscpcmg.h>

 34: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,void*);

 38: int main(int argc,char **argv)
 39: {
 40:   SNES           snes;
 42:   PetscInt       its,lits;
 43:   PetscReal      litspit;
 44:   DM             da;

 46:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

 48:   /*
 49:       Set the DMDA (grid structure) for the grids.
 50:   */
 51:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
 52:   DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
 53:   SNESCreate(PETSC_COMM_WORLD,&snes);
 54:   SNESSetDM(snes,da);
 55:   DMDestroy(&da);

 57:   SNESSetFromOptions(snes);

 59:   SNESSolve(snes,0,0);
 60:   SNESGetIterationNumber(snes,&its);
 61:   SNESGetLinearSolveIterations(snes,&lits);
 62:   litspit = ((PetscReal)lits)/((PetscReal)its);
 63:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
 64:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
 65:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",litspit);

 67:   SNESDestroy(&snes);
 68:   PetscFinalize();

 70:   return 0;
 71: }

 75: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr)
 76: {
 77:   PetscInt     i,j;
 78:   PetscScalar  hx,hy;
 79:   PetscScalar  gradup,graddown,gradleft,gradright,gradx,grady;
 80:   PetscScalar  coeffup,coeffdown,coeffleft,coeffright;

 83:   hx    = 1.0/(PetscReal)(info->mx-1);  hy    = 1.0/(PetscReal)(info->my-1);
 84: 
 85:   /* Evaluate function */
 86:   for (j=info->ys; j<info->ys+info->ym; j++) {
 87:     for (i=info->xs; i<info->xs+info->xm; i++) {

 89:       if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {

 91:         f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
 92: 
 93:       } else {

 95:         gradup     = (t[j+1][i] - t[j][i])/hy;
 96:         graddown   = (t[j][i] - t[j-1][i])/hy;
 97:         gradright  = (t[j][i+1] - t[j][i])/hx;
 98:         gradleft   = (t[j][i] - t[j][i-1])/hx;

100:         gradx      = .5*(t[j][i+1] - t[j][i-1])/hx;
101:         grady      = .5*(t[j+1][i] - t[j-1][i])/hy;

103:         coeffup    = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
104:         coeffdown  = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);

106:         coeffleft  = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
107:         coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);

109:         f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
110: 
111:       }

113:     }
114:   }
115:   return(0);
116: }