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