Actual source code: ex22.c
2: static char help[] = "Solves PDE optimization problem.\n\n";
4: #include petscda.h
5: #include petscpf.h
6: #include petscsnes.h
7: #include petscdmmg.h
9: /*
11: w - design variables (what we change to get an optimal solution)
12: u - state variables (i.e. the PDE solution)
13: lambda - the Lagrange multipliers
15: U = (w [u_0 lambda_0 u_1 lambda_1 .....])
17: fu, fw, flambda contain the gradient of L(w,u,lambda)
19: FU = (fw [fu_0 flambda_0 .....])
21: In this example the PDE is
22: Uxx = 2,
23: u(0) = w(0), thus this is the free parameter
24: u(1) = 0
25: the function we wish to minimize is
26: \integral u^{2}
28: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
30: Use the usual centered finite differences.
32: Note we treat the problem as non-linear though it happens to be linear
34: See ex21.c for the same code, but that does NOT interlaces the u and the lambda
36: The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
37: */
39: typedef struct {
40: PetscViewer u_lambda_viewer;
41: PetscViewer fu_lambda_viewer;
42: } UserCtx;
47: /*
48: Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the
49: smoother on all levels. This is because (1) in the matrix free case no matrix entries are
50: available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal
51: entry for the control variable is zero which means default SOR will not work.
53: */
54: char common_options[] = "-dmmg_grid_sequence \
55: -dmmg_nlevels 5 \
56: -mg_levels_pc_type none \
57: -mg_coarse_pc_type none \
58: -pc_mg_type full \
59: -mg_coarse_ksp_type gmres \
60: -mg_levels_ksp_type gmres \
61: -mg_coarse_ksp_max_it 6 \
62: -mg_levels_ksp_max_it 3";
64: char matrix_free_options[] = "-snes_mf_compute_normu no \
65: -snes_mf_type wp \
66: -dmmg_jacobian_mf_fd";
68: /*
69: Currently only global coloring is supported with VecPack
70: */
71: char matrix_based_options[] = "-dmmg_iscoloring_type global";
78: int main(int argc,char **argv)
79: {
81: UserCtx user;
82: DA da;
83: DMMG *dmmg;
84: VecPack packer;
85: PetscTruth use_matrix_based = PETSC_FALSE,use_monitor = PETSC_FALSE;
86: PetscInt i;
88: PetscInitialize(&argc,&argv,PETSC_NULL,help);
89: PetscOptionsSetFromOptions();
91: /* Hardwire several options; can be changed at command line */
92: PetscOptionsInsertString(common_options);
93: PetscOptionsGetTruth(PETSC_NULL,"-use_matrix_based",&use_matrix_based,PETSC_IGNORE);
94: if (use_matrix_based) {
95: PetscOptionsInsertString(matrix_based_options);
96: } else {
97: PetscOptionsInsertString(matrix_free_options);
98: }
99: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
100: PetscOptionsGetTruth(PETSC_NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);
102: /* Create a global vector that includes a single redundant array and two da arrays */
103: VecPackCreate(PETSC_COMM_WORLD,&packer);
104: VecPackAddArray(packer,0,1);
105: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,2,1,PETSC_NULL,&da);
106: VecPackAddDA(packer,da);
109: /* create nonlinear multi-level solver */
110: DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
111: DMMGSetDM(dmmg,(DM)packer);
112: DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
114: if (use_monitor) {
115: /* create graphics windows */
116: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);
117: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);
118: for (i=0; i<DMMGGetLevels(dmmg); i++) {
119: SNESMonitorSet(dmmg[i]->snes,Monitor,dmmg[i],0);
120: }
121: }
123: DMMGSolve(dmmg);
124: DMMGDestroy(dmmg);
126: DADestroy(da);
127: VecPackDestroy(packer);
128: if (use_monitor) {
129: PetscViewerDestroy(user.u_lambda_viewer);
130: PetscViewerDestroy(user.fu_lambda_viewer);
131: }
133: PetscFinalize();
134: return 0;
135: }
137: typedef struct {
138: PetscScalar u;
139: PetscScalar lambda;
140: } ULambda;
141:
142: /*
143: Evaluates FU = Gradiant(L(w,u,lambda))
145: This local function acts on the ghosted version of U (accessed via VecPackGetLocalVectors() and
146: VecPackScatter()) BUT the global, nonghosted version of FU (via VecPackGetAccess()).
148: */
149: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
150: {
151: DMMG dmmg = (DMMG)dummy;
153: PetscInt xs,xm,i,N,nredundant;
154: ULambda *u_lambda,*fu_lambda;
155: PetscScalar d,h,*w,*fw;
156: Vec vu_lambda,vfu_lambda;
157: DA da;
158: VecPack packer = (VecPack)dmmg->dm;
161: VecPackGetEntries(packer,&nredundant,&da);
162: VecPackGetLocalVectors(packer,&w,&vu_lambda);
163: VecPackScatter(packer,U,w,vu_lambda);
164: VecPackGetAccess(packer,FU,&fw,&vfu_lambda);
166: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
167: DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
168: DAVecGetArray(da,vu_lambda,&u_lambda);
169: DAVecGetArray(da,vfu_lambda,&fu_lambda);
170: d = N-1.0;
171: h = 1.0/d;
173: /* derivative of L() w.r.t. w */
174: if (xs == 0) { /* only first processor computes this */
175: fw[0] = -2.0*d*u_lambda[0].lambda;
176: }
178: /* derivative of L() w.r.t. u */
179: for (i=xs; i<xs+xm; i++) {
180: if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda;
181: else if (i == 1) fu_lambda[1].lambda = 2.*h*u_lambda[1].u + 2.*d*u_lambda[1].lambda - d*u_lambda[2].lambda;
182: else if (i == N-1) fu_lambda[N-1].lambda = h*u_lambda[N-1].u + 2.*d*u_lambda[N-1].lambda - d*u_lambda[N-2].lambda;
183: else if (i == N-2) fu_lambda[N-2].lambda = 2.*h*u_lambda[N-2].u + 2.*d*u_lambda[N-2].lambda - d*u_lambda[N-3].lambda;
184: else fu_lambda[i].lambda = 2.*h*u_lambda[i].u - d*(u_lambda[i+1].lambda - 2.0*u_lambda[i].lambda + u_lambda[i-1].lambda);
185: }
187: /* derivative of L() w.r.t. lambda */
188: for (i=xs; i<xs+xm; i++) {
189: if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]);
190: else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
191: else fu_lambda[i].u = -(d*(u_lambda[i+1].u - 2.0*u_lambda[i].u + u_lambda[i-1].u) - 2.0*h);
192: }
194: DAVecRestoreArray(da,vu_lambda,&u_lambda);
195: DAVecRestoreArray(da,vfu_lambda,&fu_lambda);
196: VecPackRestoreLocalVectors(packer,&w,&vu_lambda);
197: VecPackRestoreAccess(packer,FU,&fw,&vfu_lambda);
198: PetscLogFlops(13*N);
199: return(0);
200: }
202: /*
203: Computes the exact solution
204: */
205: PetscErrorCode u_solution(void *dummy,PetscInt n,PetscScalar *x,PetscScalar *u)
206: {
207: PetscInt i;
209: for (i=0; i<n; i++) {
210: u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
211: }
212: return(0);
213: }
215: PetscErrorCode ExactSolution(VecPack packer,Vec U)
216: {
217: PF pf;
218: Vec x,u_global;
219: PetscScalar *w;
220: DA da;
222: PetscInt m;
225: VecPackGetEntries(packer,&m,&da);
227: PFCreate(PETSC_COMM_WORLD,1,2,&pf);
228: PFSetType(pf,PFQUICK,(void*)u_solution);
229: DAGetCoordinates(da,&x);
230: if (!x) {
231: DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
232: DAGetCoordinates(da,&x);
233: }
234: VecPackGetAccess(packer,U,&w,&u_global,0);
235: if (w) w[0] = .25;
236: PFApplyVec(pf,x,u_global);
237: PFDestroy(pf);
238: VecPackRestoreAccess(packer,U,&w,&u_global,0);
239: return(0);
240: }
243: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
244: {
245: DMMG dmmg = (DMMG)dummy;
246: UserCtx *user = (UserCtx*)dmmg->user;
248: PetscInt m,N;
249: PetscScalar *w,*dw;
250: Vec u_lambda,U,F,Uexact;
251: VecPack packer = (VecPack)dmmg->dm;
252: PetscReal norm;
253: DA da;
256: SNESGetSolution(snes,&U);
257: VecPackGetAccess(packer,U,&w,&u_lambda);
258: VecView(u_lambda,user->u_lambda_viewer);
259: VecPackRestoreAccess(packer,U,&w,&u_lambda);
261: SNESGetFunction(snes,&F,0,0);
262: VecPackGetAccess(packer,F,&w,&u_lambda);
263: /* VecView(u_lambda,user->fu_lambda_viewer); */
264: VecPackRestoreAccess(packer,U,&w,&u_lambda);
266: VecPackGetEntries(packer,&m,&da);
267: DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
268: VecDuplicate(U,&Uexact);
269: ExactSolution(packer,Uexact);
270: VecAXPY(Uexact,-1.0,U);
271: VecPackGetAccess(packer,Uexact,&dw,&u_lambda);
272: VecStrideNorm(u_lambda,0,NORM_2,&norm);
273: norm = norm/sqrt(N-1.);
274: if (dw) PetscPrintf(dmmg->comm,"Norm of error %G Error at x = 0 %G\n",norm,PetscRealPart(dw[0]));
275: VecView(u_lambda,user->fu_lambda_viewer);
276: VecPackRestoreAccess(packer,Uexact,&dw,&u_lambda);
277: VecDestroy(Uexact);
278: return(0);
279: }