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 /*T 8 Concepts: SNES^solving a system of nonlinear equations 9 Concepts: DMDA^using distributed arrays 10 Concepts: multigrid; 11 Processors: n 12 T*/ 13 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 */ 29 30 #include <petscsnes.h> 31 #include <petscdm.h> 32 #include <petscdmda.h> 33 34 extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,void*); 35 36 int main(int argc,char **argv) 37 { 38 SNES snes; 39 PetscErrorCode ierr; 40 PetscInt its,lits; 41 PetscReal litspit; 42 DM da; 43 44 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 45 /* 46 Set the DMDA (grid structure) for the grids. 47 */ 48 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 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 50 ierr = DMSetUp(da);CHKERRQ(ierr); 51 ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,NULL);CHKERRQ(ierr); 52 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 53 ierr = SNESSetDM(snes,da);CHKERRQ(ierr); 54 ierr = DMDestroy(&da);CHKERRQ(ierr); 55 56 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 57 58 ierr = SNESSolve(snes,0,0);CHKERRQ(ierr); 59 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 60 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 61 litspit = ((PetscReal)lits)/((PetscReal)its); 62 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr); 63 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);CHKERRQ(ierr); 64 ierr = PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);CHKERRQ(ierr); 65 66 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 67 ierr = PetscFinalize(); 68 return ierr; 69 } 70 71 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr) 72 { 73 PetscInt i,j; 74 PetscScalar hx,hy; 75 PetscScalar gradup,graddown,gradleft,gradright,gradx,grady; 76 PetscScalar coeffup,coeffdown,coeffleft,coeffright; 77 78 PetscFunctionBeginUser; 79 hx = 1.0/(PetscReal)(info->mx-1); hy = 1.0/(PetscReal)(info->my-1); 80 81 /* Evaluate function */ 82 for (j=info->ys; j<info->ys+info->ym; j++) { 83 for (i=info->xs; i<info->xs+info->xm; i++) { 84 85 if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) { 86 f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0)); 87 } else { 88 89 gradup = (t[j+1][i] - t[j][i])/hy; 90 graddown = (t[j][i] - t[j-1][i])/hy; 91 gradright = (t[j][i+1] - t[j][i])/hx; 92 gradleft = (t[j][i] - t[j][i-1])/hx; 93 94 gradx = .5*(t[j][i+1] - t[j][i-1])/hx; 95 grady = .5*(t[j+1][i] - t[j-1][i])/hy; 96 97 coeffup = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx); 98 coeffdown = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx); 99 100 coeffleft = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady); 101 coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady); 102 103 f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy; 104 } 105 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 112 /*TEST 113 114 test: 115 args: -pc_type mg -da_refine 1 -ksp_type fgmres 116 117 test: 118 suffix: 2 119 nsize: 2 120 args: -pc_type mg -da_refine 1 -ksp_type fgmres 121 122 TEST*/ 123