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