/*T Concepts: KSP^solving a system of linear equations Concepts: KSP^Laplacian, 3d Processors: n T*/ /* Laplacian in 3D. Modeled by the partial differential equation div grad u = f, 0 < x,y,z < 1, with pure Neumann boundary conditions u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. The functions are cell-centered This uses multigrid to solve the linear system Contributed by Jianming Yang */ static char help[] = "Solves 3D Laplacian using multigrid.\n\n"; #include #include #include extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*); extern PetscErrorCode ComputeRHS(KSP,Vec,void*); int main(int argc,char **argv) { KSP ksp; DM da; PetscReal norm; PetscErrorCode ierr; PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,d,dof; PetscScalar Hx,Hy,Hz; PetscScalar ****array; Vec x,b,r; Mat J; ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; dof = 1; ierr = PetscOptionsGetInt(NULL,NULL,"-da_dof",&dof,NULL);CHKERRQ(ierr); 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,12,12,12,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,0,&da);CHKERRQ(ierr); ierr = DMSetFromOptions(da);CHKERRQ(ierr); ierr = DMSetUp(da);CHKERRQ(ierr); ierr = DMDASetInterpolationType(da, DMDA_Q0);CHKERRQ(ierr); ierr = KSPSetDM(ksp,da);CHKERRQ(ierr); ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr); ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);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 = KSPGetOperators(ksp,NULL,&J);CHKERRQ(ierr); ierr = VecDuplicate(b,&r);CHKERRQ(ierr); ierr = MatMult(J,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",(double)norm);CHKERRQ(ierr); ierr = DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); Hx = 1.0 / (PetscReal)(mx); Hy = 1.0 / (PetscReal)(my); Hz = 1.0 / (PetscReal)(mz); ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); ierr = DMDAVecGetArrayDOF(da, x, &array);CHKERRQ(ierr); for (k=zs; k