# petsc / src / ksp / ksp / examples / tutorials / ex45.c

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153``` ```/* Laplacian in 3D. Modeled by the partial differential equation - Laplacian u = 1,0 < x,y,z < 1, with boundary conditions u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. This uses multigrid to solve the linear system See src/snes/examples/tutorials/ex50.c Can also be run with -pc_type exotic -ksp_type fgmres */ static char help[] = "Solves 3D Laplacian using multigrid.\n\n"; #include #include #include extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*); extern PetscErrorCode ComputeRHS(KSP,Vec,void*); extern PetscErrorCode ComputeInitialGuess(KSP,Vec,void*); #undef __FUNCT__ #define __FUNCT__ "main" int main(int argc,char **argv) { PetscErrorCode ierr; KSP ksp; PetscReal norm; DM da; Vec x,b,r; Mat A; PetscInitialize(&argc,&argv,(char*)0,help); ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-7,-7,-7,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);CHKERRQ(ierr); ierr = KSPSetDM(ksp,da);CHKERRQ(ierr); ierr = KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,NULL);CHKERRQ(ierr); ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr); ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr); ierr = DMDestroy(&da);CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr); ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr); ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr); ierr = VecDuplicate(b,&r);CHKERRQ(ierr); ierr = KSPGetOperators(ksp,&A,NULL,NULL);CHKERRQ(ierr); ierr = MatMult(A,x,r);CHKERRQ(ierr); ierr = VecAXPY(r,-1.0,b);CHKERRQ(ierr); ierr = VecNorm(r,NORM_2,&norm);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); ierr = KSPDestroy(&ksp);CHKERRQ(ierr); ierr = PetscFinalize(); return 0; } #undef __FUNCT__ #define __FUNCT__ "ComputeRHS" PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx) { PetscErrorCode ierr; PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs; DM dm; PetscScalar Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy; PetscScalar ***barray; PetscFunctionBeginUser; ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr); ierr = DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1); HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx; ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); ierr = DMDAVecGetArray(dm,b,&barray);CHKERRQ(ierr); for (k=zs; k
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.