Actual source code: ex21.c
2: static char help[] = "Solves PDE optimization problem.\n\n";
4: #include petscda.h
5: #include petscpf.h
6: #include petscsnes.h
8: /*
10: w - design variables (what we change to get an optimal solution)
11: u - state variables (i.e. the PDE solution)
12: lambda - the Lagrange multipliers
14: U = (w u lambda)
16: fu, fw, flambda contain the gradient of L(w,u,lambda)
18: FU = (fw fu flambda)
20: In this example the PDE is
21: Uxx = 2,
22: u(0) = w(0), thus this is the free parameter
23: u(1) = 0
24: the function we wish to minimize is
25: \integral u^{2}
27: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
29: Use the usual centered finite differences.
31: Note we treat the problem as non-linear though it happens to be linear
33: See ex22.c for the same code, but that interlaces the u and the lambda
35: */
37: typedef struct {
38: DA da1,da2;
39: PetscInt nredundant;
40: DMComposite packer;
41: PetscViewer u_viewer,lambda_viewer;
42: PetscViewer fu_viewer,flambda_viewer;
43: } UserCtx;
51: int main(int argc,char **argv)
52: {
54: PetscInt its;
55: Vec U,FU;
56: SNES snes;
57: UserCtx user;
59: PetscInitialize(&argc,&argv,(char*)0,help);
61: /* Create a global vector that includes a single redundant array and two da arrays */
62: DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);
63: user.nredundant = 1;
64: DMCompositeAddArray(user.packer,0,user.nredundant);
65: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&user.da1);
66: DMCompositeAddDM(user.packer,(DM)user.da1);
67: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&user.da2);
68: DMCompositeAddDM(user.packer,(DM)user.da2);
69: DMCompositeCreateGlobalVector(user.packer,&U);
70: VecDuplicate(U,&FU);
72: /* create graphics windows */
73: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
74: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);
75: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);
76: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);
79: /* create nonlinear solver */
80: SNESCreate(PETSC_COMM_WORLD,&snes);
81: SNESSetFunction(snes,FU,FormFunction,&user);
82: SNESSetFromOptions(snes);
83: SNESMonitorSet(snes,Monitor,&user,0);
84: SNESSolve(snes,PETSC_NULL,U);
85: SNESGetIterationNumber(snes,&its);
86: SNESDestroy(snes);
88: DADestroy(user.da1);
89: DADestroy(user.da2);
90: DMCompositeDestroy(user.packer);
91: VecDestroy(U);
92: VecDestroy(FU);
93: PetscViewerDestroy(user.u_viewer);
94: PetscViewerDestroy(user.lambda_viewer);
95: PetscViewerDestroy(user.fu_viewer);
96: PetscViewerDestroy(user.flambda_viewer);
97: PetscFinalize();
98: return 0;
99: }
100:
101: /*
102: Evaluates FU = Gradiant(L(w,u,lambda))
104: */
105: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
106: {
107: UserCtx *user = (UserCtx*)dummy;
109: PetscInt xs,xm,i,N;
110: PetscScalar *u,*lambda,*w,*fu,*fw,*flambda,d,h;
111: Vec vu,vlambda,vfu,vflambda;
114: DMCompositeGetLocalVectors(user->packer,&w,&vu,&vlambda);
115: DMCompositeGetLocalVectors(user->packer,&fw,&vfu,&vflambda);
116: DMCompositeScatter(user->packer,U,w,vu,vlambda);
118: DAGetCorners(user->da1,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
119: DAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0);
120: DAVecGetArray(user->da1,vu,&u);
121: DAVecGetArray(user->da1,vfu,&fu);
122: DAVecGetArray(user->da1,vlambda,&lambda);
123: DAVecGetArray(user->da1,vflambda,&flambda);
124: d = (N-1.0);
125: h = 1.0/d;
127: /* derivative of L() w.r.t. w */
128: if (xs == 0) { /* only first processor computes this */
129: fw[0] = -2.*d*lambda[0];
130: }
132: /* derivative of L() w.r.t. u */
133: for (i=xs; i<xs+xm; i++) {
134: if (i == 0) flambda[0] = h*u[0] + 2.*d*lambda[0] - d*lambda[1];
135: else if (i == 1) flambda[1] = 2.*h*u[1] + 2.*d*lambda[1] - d*lambda[2];
136: else if (i == N-1) flambda[N-1] = h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2];
137: else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3];
138: else flambda[i] = 2.*h*u[i] - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]);
139: }
141: /* derivative of L() w.r.t. lambda */
142: for (i=xs; i<xs+xm; i++) {
143: if (i == 0) fu[0] = 2.0*d*(u[0] - w[0]);
144: else if (i == N-1) fu[N-1] = 2.0*d*u[N-1];
145: else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);
146: }
148: DAVecRestoreArray(user->da1,vu,&u);
149: DAVecRestoreArray(user->da1,vfu,&fu);
150: DAVecRestoreArray(user->da1,vlambda,&lambda);
151: DAVecRestoreArray(user->da1,vflambda,&flambda);
153: DMCompositeGather(user->packer,FU,fw,vfu,vflambda);
154: DMCompositeRestoreLocalVectors(user->packer,&w,&vu,&vlambda);
155: DMCompositeRestoreLocalVectors(user->packer,&fw,&vfu,&vflambda);
156: return(0);
157: }
159: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
160: {
161: UserCtx *user = (UserCtx*)dummy;
163: PetscScalar *w;
164: Vec u,lambda,U,F;
167: SNESGetSolution(snes,&U);
168: DMCompositeGetAccess(user->packer,U,&w,&u,&lambda);
169: VecView(u,user->u_viewer);
170: VecView(lambda,user->lambda_viewer);
171: DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda);
173: SNESGetFunction(snes,&F,0,0);
174: DMCompositeGetAccess(user->packer,F,&w,&u,&lambda);
175: VecView(u,user->fu_viewer);
176: VecView(lambda,user->flambda_viewer);
177: DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda);
178: return(0);
179: }