Actual source code: ex23.c
2: static char help[] = "Solves PDE problem from ex22.c\n\n";
4: #include <petscdmda.h>
5: #include <petscpf.h>
6: #include <petscsnes.h>
7: #include <petscdmmg.h>
9: /*
11: In this example the PDE is
12: Uxx + U^2 = 2,
13: u(0) = .25
14: u(1) = 0
16: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
18: Use the usual centered finite differences.
21: See ex22.c for a design optimization problem
23: */
25: typedef struct {
26: PetscViewer u_viewer;
27: PetscViewer fu_viewer;
28: } UserCtx;
35: int main(int argc,char **argv)
36: {
38: UserCtx user;
39: DM da;
40: DMMG *dmmg;
42: PetscInitialize(&argc,&argv,PETSC_NULL,help);
44: /* Hardwire several options; can be changed at command line */
45: PetscOptionsSetValue("-dmmg_grid_sequence",PETSC_NULL);
46: PetscOptionsSetValue("-ksp_type","fgmres");
47: PetscOptionsSetValue("-ksp_max_it","5");
48: PetscOptionsSetValue("-pc_mg_type","full");
49: PetscOptionsSetValue("-mg_coarse_ksp_type","gmres");
50: PetscOptionsSetValue("-mg_levels_ksp_type","gmres");
51: PetscOptionsSetValue("-mg_coarse_ksp_max_it","3");
52: PetscOptionsSetValue("-mg_levels_ksp_max_it","3");
53: /* PetscOptionsSetValue("-snes_mf_type","wp"); */
54: /* PetscOptionsSetValue("-snes_mf_compute_normu","no"); */
55: PetscOptionsSetValue("-snes_ls","basic");
56: /* PetscOptionsSetValue("-dmmg_jacobian_mf_fd",0); */
57: /* PetscOptionsSetValue("-snes_ls","basicnonorms"); */
58: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
59:
60: /* Create a global vector from a da arrays */
61: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-5,1,1,PETSC_NULL,&da);
63: /* create graphics windows */
64: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
65: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - discretization of function",-1,-1,-1,-1,&user.fu_viewer);
67: /* create nonlinear multi-level solver */
68: DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
69: DMMGSetDM(dmmg,(DM)da);
70: DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
71: DMMGSetFromOptions(dmmg);
72: DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,0);
73: DMMGSolve(dmmg);
74: DMMGDestroy(dmmg);
76: DMDestroy(&da);
77: PetscViewerDestroy(&user.u_viewer);
78: PetscViewerDestroy(&user.fu_viewer);
80: PetscFinalize();
81: return 0;
82: }
86: /*
87: This local function acts on the ghosted version of U (accessed via DMGetLocalVector())
88: BUT the global, nonghosted version of FU
90: */
91: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
92: {
93: DMMG dmmg = (DMMG)dummy;
95: PetscInt xs,xm,i,N;
96: PetscScalar *u,*fu,d,h;
97: Vec vu;
98: DM da = dmmg->dm;
101: DMGetLocalVector(da,&vu);
102: DMGlobalToLocalBegin(da,U,INSERT_VALUES,vu);
103: DMGlobalToLocalEnd(da,U,INSERT_VALUES,vu);
105: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
106: DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
107: DMDAVecGetArray(da,vu,&u);
108: DMDAVecGetArray(da,FU,&fu);
109: d = N-1.0;
110: h = 1.0/d;
112: for (i=xs; i<xs+xm; i++) {
113: if (i == 0) fu[0] = 2.0*d*(u[0] - .25) + h*u[0]*u[0];
114: else if (i == N-1) fu[N-1] = 2.0*d*u[N-1] + h*u[N-1]*u[N-1];
115: else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];
116: }
118: DMDAVecRestoreArray(da,vu,&u);
119: DMDAVecRestoreArray(da,FU,&fu);
120: DMRestoreLocalVector(da,&vu);
121: PetscLogFlops(9.0*N);
122: return(0);
123: }
127: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *pt,PetscScalar *u,PetscScalar *fu,void* dummy)
128: {
129: PetscInt i,N = info->mx;
130: PetscScalar d,h;
133: d = N-1.0;
134: h = 1.0/d;
136: i = pt->i;
137: if (i == 0) *fu = 2.0*d*(u[0] - .25) + h*u[0]*u[0];
138: else if (i == N-1) *fu = 2.0*d*u[N-1] + h*u[N-1]*u[N-1];
139: else *fu = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];
141: return(0);
142: }