Actual source code: ex22.c
2: static char help[] = "Solves PDE optimization problem.\n\n";
4: #include <petscdmda.h>
5: #include <petscdmredundant.h>
6: #include <petscdmcomposite.h>
7: #include <petscpf.h>
8: #include <petscsnes.h>
9: #include <private/dmimpl.h>
11: /*
13: w - design variables (what we change to get an optimal solution)
14: u - state variables (i.e. the PDE solution)
15: lambda - the Lagrange multipliers
17: U = (w [u_0 lambda_0 u_1 lambda_1 .....])
19: fu, fw, flambda contain the gradient of L(w,u,lambda)
21: FU = (fw [fu_0 flambda_0 .....])
23: In this example the PDE is
24: Uxx = 2,
25: u(0) = w(0), thus this is the free parameter
26: u(1) = 0
27: the function we wish to minimize is
28: \integral u^{2}
30: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
32: Use the usual centered finite differences.
34: Note we treat the problem as non-linear though it happens to be linear
36: See ex21.c for the same code, but that does NOT interlaces the u and the lambda
38: The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
39: */
41: typedef struct {
42: PetscViewer u_lambda_viewer;
43: PetscViewer fu_lambda_viewer;
44: } UserCtx;
46: extern PetscErrorCode ComputeFunction(DM,Vec,Vec);
47: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
49: /*
50: Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the
51: smoother on all levels. This is because (1) in the matrix free case no matrix entries are
52: available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal
53: entry for the control variable is zero which means default SOR will not work.
55: */
56: char common_options[] = "-ksp_type fgmres\
57: -snes_grid_sequence 5 \
58: -pc_type mg\
59: -mg_levels_pc_type none \
60: -mg_coarse_pc_type none \
61: -pc_mg_type full \
62: -mg_coarse_ksp_type gmres \
63: -mg_levels_ksp_type gmres \
64: -mg_coarse_ksp_max_it 6 \
65: -mg_levels_ksp_max_it 3";
67: char matrix_free_options[] = "-mat_mffd_compute_normu no \
68: -mat_mffd_type wp";
70: extern PetscErrorCode DMCreateMatrix_MF(DM,const MatType,Mat*);
71: extern PetscErrorCode DMComputeJacobian_MF(DM,Vec,Mat,Mat,MatStructure *);
75: int main(int argc,char **argv)
76: {
78: UserCtx user;
79: DM red,da;
80: SNES snes;
81: DM packer;
82: PetscBool use_monitor = PETSC_FALSE;
84: PetscInitialize(&argc,&argv,PETSC_NULL,help);
85: PetscOptionsSetFromOptions();
87: /* Hardwire several options; can be changed at command line */
88: PetscOptionsInsertString(common_options);
89: PetscOptionsInsertString(matrix_free_options);
90: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
91: PetscOptionsGetBool(PETSC_NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);
93: /* Create a global vector that includes a single redundant array and two da arrays */
94: DMCompositeCreate(PETSC_COMM_WORLD,&packer);
95: DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);
96: DMCompositeAddDM(packer,red);
97: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-5,2,1,PETSC_NULL,&da);
98: DMCompositeAddDM(packer,(DM)da);
99: DMSetApplicationContext(packer,&user);
100: DMSetFunction(packer,ComputeFunction);
101: DMSetJacobian(packer,DMComputeJacobian_MF);
102: packer->ops->creatematrix = DMCreateMatrix_MF;
104: /* create nonlinear multi-level solver */
105: SNESCreate(PETSC_COMM_WORLD,&snes);
106: SNESSetDM(snes,(DM)packer);
107: SNESSetFromOptions(snes);
109: if (use_monitor) {
110: /* create graphics windows */
111: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);
112: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);
113: SNESMonitorSet(snes,Monitor,0,0);
114: }
116: SNESSolve(snes,PETSC_NULL,PETSC_NULL);
117: SNESDestroy(&snes);
119: DMDestroy(&red);
120: DMDestroy(&da);
121: DMDestroy(&packer);
122: if (use_monitor) {
123: PetscViewerDestroy(&user.u_lambda_viewer);
124: PetscViewerDestroy(&user.fu_lambda_viewer);
125: }
126: PetscFinalize();
127: return 0;
128: }
130: typedef struct {
131: PetscScalar u;
132: PetscScalar lambda;
133: } ULambda;
137: /*
138: Evaluates FU = Gradiant(L(w,u,lambda))
140: This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
141: DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).
143: */
144: PetscErrorCode ComputeFunction(DM packer,Vec U,Vec FU)
145: {
147: PetscInt xs,xm,i,N;
148: ULambda *u_lambda,*fu_lambda;
149: PetscScalar d,h,*w,*fw;
150: Vec vw,vfw,vu_lambda,vfu_lambda;
151: DM red,da;
154: DMCompositeGetEntries(packer,&red,&da);
155: DMCompositeGetLocalVectors(packer,&vw,&vu_lambda);
156: DMCompositeScatter(packer,U,vw,vu_lambda);
157: DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda);
159: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
160: DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
161: VecGetArray(vw,&w);
162: VecGetArray(vfw,&fw);
163: DMDAVecGetArray(da,vu_lambda,&u_lambda);
164: DMDAVecGetArray(da,vfu_lambda,&fu_lambda);
165: d = N-1.0;
166: h = 1.0/d;
168: /* derivative of L() w.r.t. w */
169: if (xs == 0) { /* only first processor computes this */
170: fw[0] = -2.0*d*u_lambda[0].lambda;
171: }
173: /* derivative of L() w.r.t. u */
174: for (i=xs; i<xs+xm; i++) {
175: if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda;
176: 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;
177: 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;
178: 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;
179: 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);
180: }
182: /* derivative of L() w.r.t. lambda */
183: for (i=xs; i<xs+xm; i++) {
184: if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]);
185: else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
186: 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);
187: }
189: VecRestoreArray(vw,&w);
190: VecRestoreArray(vfw,&fw);
191: DMDAVecRestoreArray(da,vu_lambda,&u_lambda);
192: DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);
193: DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda);
194: DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda);
195: PetscLogFlops(13.0*N);
196: return(0);
197: }
201: /*
202: Computes the exact solution
203: */
204: PetscErrorCode u_solution(void *dummy,PetscInt n,const PetscScalar *x,PetscScalar *u)
205: {
206: PetscInt i;
208: for (i=0; i<n; i++) {
209: u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
210: }
211: return(0);
212: }
216: PetscErrorCode ExactSolution(DM packer,Vec U)
217: {
218: PF pf;
219: Vec x,u_global;
220: PetscScalar *w;
221: DM da;
223: PetscInt m;
226: DMCompositeGetEntries(packer,&m,&da);
228: PFCreate(PETSC_COMM_WORLD,1,2,&pf);
229: PFSetType(pf,PFQUICK,(void*)u_solution);
230: DMDAGetCoordinates(da,&x);
231: if (!x) {
232: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
233: DMDAGetCoordinates(da,&x);
234: }
235: DMCompositeGetAccess(packer,U,&w,&u_global,0);
236: if (w) w[0] = .25;
237: PFApplyVec(pf,x,u_global);
238: PFDestroy(&pf);
239: DMCompositeRestoreAccess(packer,U,&w,&u_global,0);
240: return(0);
241: }
245: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
246: {
247: UserCtx *user;
249: PetscInt m,N;
250: PetscScalar *w,*dw;
251: Vec u_lambda,U,F,Uexact;
252: DM packer;
253: PetscReal norm;
254: DM da;
257: SNESGetDM(snes,&packer);
258: DMGetApplicationContext(packer,&user);
259: SNESGetSolution(snes,&U);
260: DMCompositeGetAccess(packer,U,&w,&u_lambda);
261: VecView(u_lambda,user->u_lambda_viewer);
262: DMCompositeRestoreAccess(packer,U,&w,&u_lambda);
264: SNESGetFunction(snes,&F,0,0);
265: DMCompositeGetAccess(packer,F,&w,&u_lambda);
266: /* VecView(u_lambda,user->fu_lambda_viewer); */
267: DMCompositeRestoreAccess(packer,U,&w,&u_lambda);
269: DMCompositeGetEntries(packer,&m,&da);
270: DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
271: VecDuplicate(U,&Uexact);
272: ExactSolution(packer,Uexact);
273: VecAXPY(Uexact,-1.0,U);
274: DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);
275: VecStrideNorm(u_lambda,0,NORM_2,&norm);
276: norm = norm/sqrt(N-1.);
277: if (dw) PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G Error at x = 0 %G\n",norm,PetscRealPart(dw[0]));
278: VecView(u_lambda,user->fu_lambda_viewer);
279: DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);
280: VecDestroy(&Uexact);
281: return(0);
282: }
286: PetscErrorCode DMCreateMatrix_MF(DM packer,const MatType stype,Mat *A )
287: {
289: DM da,red;
290: PetscInt m,M,dof;
293: DMCompositeGetEntries(packer,&red,&da);
294: DMRedundantGetSize(red,PETSC_NULL,&m);
295: DMDAGetInfo(da,0,&M,0,0,0,0,0,&dof,0,0,0,0,0);
296: MatCreateMFFD(PETSC_COMM_WORLD,m+M*dof,m+M*dof,PETSC_DETERMINE,PETSC_DETERMINE,A);
297: MatMFFDSetFunction(*A,(PetscErrorCode (*)(void*, Vec,Vec))DMComputeFunction,packer);
298: MatSetUp(*A);
299: return(0);
300: }
304: PetscErrorCode DMComputeJacobian_MF(DM packer,Vec x,Mat A,Mat B,MatStructure *str )
305: {
309: MatMFFDSetBase(A,x,PETSC_NULL);
310: return(0);
311: }