xref: /petsc/src/snes/tutorials/ex25.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
1*c4762a1bSJed Brown static const char help[] ="Minimum surface problem in 2D.\n\
2*c4762a1bSJed Brown Uses 2-dimensional distributed arrays.\n\
3*c4762a1bSJed Brown \n\
4*c4762a1bSJed Brown   Solves the linear systems via multilevel methods \n\
5*c4762a1bSJed Brown \n\n";
6*c4762a1bSJed Brown 
7*c4762a1bSJed Brown /*T
8*c4762a1bSJed Brown    Concepts: SNES^solving a system of nonlinear equations
9*c4762a1bSJed Brown    Concepts: DMDA^using distributed arrays
10*c4762a1bSJed Brown    Concepts: multigrid;
11*c4762a1bSJed Brown    Processors: n
12*c4762a1bSJed Brown T*/
13*c4762a1bSJed Brown 
14*c4762a1bSJed Brown /*
15*c4762a1bSJed Brown 
16*c4762a1bSJed Brown     This example models the partial differential equation
17*c4762a1bSJed Brown 
18*c4762a1bSJed Brown          - Div((1 + ||GRAD T||^2)^(1/2) (GRAD T)) = 0.
19*c4762a1bSJed Brown 
20*c4762a1bSJed Brown 
21*c4762a1bSJed Brown     in the unit square, which is uniformly discretized in each of x and
22*c4762a1bSJed Brown     y in this simple encoding.  The degrees of freedom are vertex centered
23*c4762a1bSJed Brown 
24*c4762a1bSJed Brown     A finite difference approximation with the usual 5-point stencil
25*c4762a1bSJed Brown     is used to discretize the boundary value problem to obtain a
26*c4762a1bSJed Brown     nonlinear system of equations.
27*c4762a1bSJed Brown 
28*c4762a1bSJed Brown */
29*c4762a1bSJed Brown 
30*c4762a1bSJed Brown #include <petscsnes.h>
31*c4762a1bSJed Brown #include <petscdm.h>
32*c4762a1bSJed Brown #include <petscdmda.h>
33*c4762a1bSJed Brown 
34*c4762a1bSJed Brown extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,void*);
35*c4762a1bSJed Brown 
36*c4762a1bSJed Brown int main(int argc,char **argv)
37*c4762a1bSJed Brown {
38*c4762a1bSJed Brown   SNES           snes;
39*c4762a1bSJed Brown   PetscErrorCode ierr;
40*c4762a1bSJed Brown   PetscInt       its,lits;
41*c4762a1bSJed Brown   PetscReal      litspit;
42*c4762a1bSJed Brown   DM             da;
43*c4762a1bSJed Brown 
44*c4762a1bSJed Brown   ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
45*c4762a1bSJed Brown   /*
46*c4762a1bSJed Brown       Set the DMDA (grid structure) for the grids.
47*c4762a1bSJed Brown   */
48*c4762a1bSJed Brown   ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);
49*c4762a1bSJed Brown   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
50*c4762a1bSJed Brown   ierr = DMSetUp(da);CHKERRQ(ierr);
51*c4762a1bSJed Brown   ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,NULL);CHKERRQ(ierr);
52*c4762a1bSJed Brown   ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
53*c4762a1bSJed Brown   ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
54*c4762a1bSJed Brown   ierr = DMDestroy(&da);CHKERRQ(ierr);
55*c4762a1bSJed Brown 
56*c4762a1bSJed Brown   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
57*c4762a1bSJed Brown 
58*c4762a1bSJed Brown   ierr    = SNESSolve(snes,0,0);CHKERRQ(ierr);
59*c4762a1bSJed Brown   ierr    = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
60*c4762a1bSJed Brown   ierr    = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr);
61*c4762a1bSJed Brown   litspit = ((PetscReal)lits)/((PetscReal)its);
62*c4762a1bSJed Brown   ierr    = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr);
63*c4762a1bSJed Brown   ierr    = PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);CHKERRQ(ierr);
64*c4762a1bSJed Brown   ierr    = PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);CHKERRQ(ierr);
65*c4762a1bSJed Brown 
66*c4762a1bSJed Brown   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
67*c4762a1bSJed Brown   ierr = PetscFinalize();
68*c4762a1bSJed Brown   return ierr;
69*c4762a1bSJed Brown }
70*c4762a1bSJed Brown 
71*c4762a1bSJed Brown PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr)
72*c4762a1bSJed Brown {
73*c4762a1bSJed Brown   PetscInt    i,j;
74*c4762a1bSJed Brown   PetscScalar hx,hy;
75*c4762a1bSJed Brown   PetscScalar gradup,graddown,gradleft,gradright,gradx,grady;
76*c4762a1bSJed Brown   PetscScalar coeffup,coeffdown,coeffleft,coeffright;
77*c4762a1bSJed Brown 
78*c4762a1bSJed Brown   PetscFunctionBeginUser;
79*c4762a1bSJed Brown   hx = 1.0/(PetscReal)(info->mx-1);  hy    = 1.0/(PetscReal)(info->my-1);
80*c4762a1bSJed Brown 
81*c4762a1bSJed Brown   /* Evaluate function */
82*c4762a1bSJed Brown   for (j=info->ys; j<info->ys+info->ym; j++) {
83*c4762a1bSJed Brown     for (i=info->xs; i<info->xs+info->xm; i++) {
84*c4762a1bSJed Brown 
85*c4762a1bSJed Brown       if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {
86*c4762a1bSJed Brown         f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
87*c4762a1bSJed Brown       } else {
88*c4762a1bSJed Brown 
89*c4762a1bSJed Brown         gradup    = (t[j+1][i] - t[j][i])/hy;
90*c4762a1bSJed Brown         graddown  = (t[j][i] - t[j-1][i])/hy;
91*c4762a1bSJed Brown         gradright = (t[j][i+1] - t[j][i])/hx;
92*c4762a1bSJed Brown         gradleft  = (t[j][i] - t[j][i-1])/hx;
93*c4762a1bSJed Brown 
94*c4762a1bSJed Brown         gradx = .5*(t[j][i+1] - t[j][i-1])/hx;
95*c4762a1bSJed Brown         grady = .5*(t[j+1][i] - t[j-1][i])/hy;
96*c4762a1bSJed Brown 
97*c4762a1bSJed Brown         coeffup   = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
98*c4762a1bSJed Brown         coeffdown = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);
99*c4762a1bSJed Brown 
100*c4762a1bSJed Brown         coeffleft  = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
101*c4762a1bSJed Brown         coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);
102*c4762a1bSJed Brown 
103*c4762a1bSJed Brown         f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
104*c4762a1bSJed Brown       }
105*c4762a1bSJed Brown 
106*c4762a1bSJed Brown     }
107*c4762a1bSJed Brown   }
108*c4762a1bSJed Brown   PetscFunctionReturn(0);
109*c4762a1bSJed Brown }
110*c4762a1bSJed Brown 
111*c4762a1bSJed Brown 
112*c4762a1bSJed Brown /*TEST
113*c4762a1bSJed Brown 
114*c4762a1bSJed Brown    test:
115*c4762a1bSJed Brown       args: -pc_type mg -da_refine 1 -ksp_type fgmres
116*c4762a1bSJed Brown 
117*c4762a1bSJed Brown    test:
118*c4762a1bSJed Brown       suffix: 2
119*c4762a1bSJed Brown       nsize: 2
120*c4762a1bSJed Brown       args: -pc_type mg -da_refine 1 -ksp_type fgmres
121*c4762a1bSJed Brown 
122*c4762a1bSJed Brown TEST*/
123