Actual source code: ex29.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Processors: n
5: T*/
7: /*
8: Added at the request of Marc Garbey.
10: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
12: div \rho grad u = f, 0 < x,y < 1,
14: with forcing function
16: f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}
18: with Dirichlet boundary conditions
20: u = f(x,y) for x = 0, x = 1, y = 0, y = 1
22: or pure Neumman boundary conditions
24: This uses multigrid to solve the linear system
25: */
27: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
29: #include petscda.h
30: #include petscksp.h
31: #include petscmg.h
32: #include petscdmmg.h
38: typedef enum {DIRICHLET, NEUMANN} BCType;
40: typedef struct {
41: PetscScalar rho;
42: PetscScalar nu;
43: BCType bcType;
44: } UserContext;
48: int main(int argc,char **argv)
49: {
50: DMMG *dmmg;
51: DA da;
52: UserContext user;
53: PetscReal norm;
54: const char *bcTypes[2] = {"dirichlet","neumann"};
56: PetscInt l,bc;
58: PetscInitialize(&argc,&argv,(char *)0,help);
60: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
61: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
62: DMMGSetDM(dmmg,(DM)da);
63: DADestroy(da);
64: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
65: DMMGSetUser(dmmg,l,&user);
66: }
68: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMMG");
69: user.rho = 1.0;
70: PetscOptionsScalar("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, PETSC_NULL);
71: user.nu = 0.1;
72: PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, PETSC_NULL);
73: bc = (PetscInt)DIRICHLET;
74: PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
75: user.bcType = (BCType)bc;
76: PetscOptionsEnd();
78: DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
79: if (user.bcType == NEUMANN) {
80: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
81: }
83: DMMGSolve(dmmg);
85: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
86: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
87: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
88: /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm); */
89: VecAssemblyBegin(DMMGGetx(dmmg));
90: VecAssemblyEnd(DMMGGetx(dmmg));
91: VecView_VTK(DMMGGetRHS(dmmg), "rhs.vtk", bcTypes[user.bcType]);
92: VecView_VTK(DMMGGetx(dmmg), "solution.vtk", bcTypes[user.bcType]);
94: DMMGDestroy(dmmg);
95: PetscFinalize();
97: return 0;
98: }
102: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
103: {
104: DA da = (DA)dmmg->dm;
105: UserContext *user = (UserContext *) dmmg->user;
107: PetscInt i,j,mx,my,xm,ym,xs,ys;
108: PetscScalar Hx,Hy;
109: PetscScalar **array;
112: DAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0);
113: Hx = 1.0 / (PetscReal)(mx-1);
114: Hy = 1.0 / (PetscReal)(my-1);
115: DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
116: DAVecGetArray(da, b, &array);
117: for (j=ys; j<ys+ym; j++){
118: for(i=xs; i<xs+xm; i++){
119: array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
120: }
121: }
122: DAVecRestoreArray(da, b, &array);
123: VecAssemblyBegin(b);
124: VecAssemblyEnd(b);
126: /* force right hand side to be consistent for singular matrix */
127: /* note this is really a hack, normally the model would provide you with a consistent right handside */
128: if (user->bcType == NEUMANN) {
129: MatNullSpace nullspace;
131: KSPGetNullSpace(dmmg->,"Norm of error %A, Iterations %D\n",norm,its);
290: /*
291: Free work space. All PETSc objects should be destroyed when they
292: are no longer needed.
293: */
294: KSPDestroy(ksp);
295: VecDestroy(u);
296: VecDestroy(x);
297: VecDestroy(b);
298: MatDestroy(C);
300: /*
301: Indicate to PETSc profiling that we're concluding the second stage
302: */
303: PetscLogStagePop();
305: PetscFinalize();
306: return 0;
307: }
./usr/share/doc/petsc2.3.2-doc/src/ksp/ksp/examples/tutorials/ex29.c.html 0000644 0000000 0000000 00000102421 10525170325 024560 0 ustar root root Actual source code: ex29.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Processors: n
5: T*/
7: /*
8: Added at the request of Marc Garbey.
10: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
12: div \rho grad u = f, 0 < x,y < 1,
14: with forcing function
16: f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}
18: with Dirichlet boundary conditions
20: u = f(x,y) for x = 0, x = 1, y = 0, y = 1
22: or pure Neumman boundary conditions
24: This uses multigrid to solve the linear system
25: */
27: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
29: #include petscda.h
30: #include petscksp.h
31: #include petscmg.h
32: #include petscdmmg.h
38: typedef enum {DIRICHLET, NEUMANN} BCType;
40: typedef struct {
41: PetscScalar rho;
42: PetscScalar nu;
43: BCType bcType;
44: } UserContext;
48: int main(int argc,char **argv)
49: {
50: DMMG *dmmg;
51: DA da;
52: UserContext user;
53: PetscReal norm;
54: const char *bcTypes[2] = {"dirichlet","neumann"};
56: PetscInt l,bc;
58: PetscInitialize(&argc,&argv,(char *)0,help);
60: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
61: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
62: DMMGSetDM(dmmg,(DM)da);
63: DADestroy(da);
64: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
65: DMMGSetUser(dmmg,l,&user);
66: }
68: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMMG");
69: user.rho = 1.0;
70: PetscOptionsScalar("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, PETSC_NULL);
71: user.nu = 0.1;
72: PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, PETSC_NULL);
73: bc = (PetscInt)DIRICHLET;
74: PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
75: user.bcType = (BCType)bc;
76: PetscOptionsEnd();
78: DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
79: if (user.bcType == NEUMANN) {
80: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
81: }
83: DMMGSolve(dmmg);
85: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
86: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
87: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
88: /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm); */
89: VecAssemblyBegin(DMMGGetx(dmmg));
90: VecAssemblyEnd(DMMGGetx(dmmg));
91: VecView_VTK(DMMGGetRHS(dmmg), "rhs.vtk", bcTypes[user.bcType]);
92: VecView_VTK(DMMGGetx(dmmg), "solution.vtk", bcTypes[user.bcType]);
94: DMMGDestroy(dmmg);
95: PetscFinalize();
97: return 0;
98: }
102: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
103: {
104: DA da = (DA)dmmg->dm;
105: UserContext *user = (UserContext *) dmmg->user;
107: PetscInt i,j,mx,my,xm,ym,xs,ys;
108: PetscScalar Hx,Hy;
109: PetscScalar **array;
112: DAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0);
113: Hx = 1.0 / (PetscReal)(mx-1);
114: Hy = 1.0 / (PetscReal)(my-1);
115: DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
116: DAVecGetArray(da, b, &array);
117: for (j=ys; j<ys+ym; j++){
118: for(i=xs; i<xs+xm; i++){
119: array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
120: }
121: }
122: DAVecRestoreArray(da, b, &array);
123: VecAssemblyBegin(b);
124: VecAssemblyEnd(b);
126: /* force right hand side to be consistent for singular matrix */
127: /* note this is really a hack, normally the model would provide you with a consistent right handside */
128: if (user->bcType == NEUMANN) {
129: MatNullSpace nullspace;
131: KSPGetNullSpace(dmmg->,"Norm of error %A, Iterations %D\n",norm,its);
290: /*
291: Free work space. All PETSc objects should be destroyed when they
292: are no longer needed.
293: */
294: KSPDestroy(ksp);
295: VecDestroy(u);
296: VecDestroy(x);
297: VecDestroy(b);
298: MatDestroy(C);
300: /*
301: Indicate to PETSc profiling that we're concluding the second stage
302: */
303: PetscLogStagePop();
305: PetscFinalize();
306: return 0;
307: }
./usr/share/doc/petsc2.3.2-doc/src/ksp/ksp/examples/tutorials/ex29.c.html 0000644 0000000 0000000 00000102421 10525170325 024560 0 ustar root root Actual source code: ex29.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Processors: n
5: T*/
7: /*
8: Added at the request of Marc Garbey.
10: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
12: div \rho grad u = f, 0 < x,y < 1,
14: with forcing function
16: f = e^{-(1 - x)^2/\nu} e^{