Actual source code: ex21.c
petsc-3.7.3 2016-08-01
2: static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
6: #include <petscdmredundant.h>
7: #include <petscdmcomposite.h>
8: #include <petscpf.h>
9: #include <petscsnes.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 lambda)
19: fu, fw, flambda contain the gradient of L(w,u,lambda)
21: FU = (fw fu flambda)
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 ex22.c for the same code, but that interlaces the u and the lambda
38: */
40: typedef struct {
41: DM red1,da1,da2;
42: DM packer;
43: PetscViewer u_viewer,lambda_viewer;
44: PetscViewer fu_viewer,flambda_viewer;
45: } UserCtx;
47: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
48: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
53: int main(int argc,char **argv)
54: {
56: PetscInt its;
57: Vec U,FU;
58: SNES snes;
59: UserCtx user;
61: PetscInitialize(&argc,&argv,(char*)0,help);
63: /* Create a global vector that includes a single redundant array and two da arrays */
64: DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);
65: DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);
66: DMCompositeAddDM(user.packer,user.red1);
67: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-5,1,1,NULL,&user.da1);
68: DMCompositeAddDM(user.packer,user.da1);
69: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-5,1,1,NULL,&user.da2);
70: DMCompositeAddDM(user.packer,user.da2);
71: DMCreateGlobalVector(user.packer,&U);
72: VecDuplicate(U,&FU);
74: /* create graphics windows */
75: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
76: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);
77: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);
78: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);
81: /* create nonlinear solver */
82: SNESCreate(PETSC_COMM_WORLD,&snes);
83: SNESSetFunction(snes,FU,FormFunction,&user);
84: SNESSetFromOptions(snes);
85: SNESMonitorSet(snes,Monitor,&user,0);
86: SNESSolve(snes,NULL,U);
87: SNESGetIterationNumber(snes,&its);
88: SNESDestroy(&snes);
90: DMDestroy(&user.red1);
91: DMDestroy(&user.da1);
92: DMDestroy(&user.da2);
93: DMDestroy(&user.packer);
94: VecDestroy(&U);
95: VecDestroy(&FU);
96: PetscViewerDestroy(&user.u_viewer);
97: PetscViewerDestroy(&user.lambda_viewer);
98: PetscViewerDestroy(&user.fu_viewer);
99: PetscViewerDestroy(&user.flambda_viewer);
100: PetscFinalize();
101: return 0;
102: }
106: /*
107: Evaluates FU = Gradiant(L(w,u,lambda))
109: */
110: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy)
111: {
112: UserCtx *user = (UserCtx*)dummy;
114: PetscInt xs,xm,i,N;
115: PetscScalar *u,*lambda,*w,*fu,*fw,*flambda,d,h;
116: Vec vw,vu,vlambda,vfw,vfu,vflambda;
119: DMCompositeGetLocalVectors(user->packer,&vw,&vu,&vlambda);
120: DMCompositeGetLocalVectors(user->packer,&vfw,&vfu,&vflambda);
121: DMCompositeScatter(user->packer,U,vw,vu,vlambda);
123: DMDAGetCorners(user->da1,&xs,NULL,NULL,&xm,NULL,NULL);
124: DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0);
125: VecGetArray(vw,&w);
126: VecGetArray(vfw,&fw);
127: DMDAVecGetArray(user->da1,vu,&u);
128: DMDAVecGetArray(user->da1,vfu,&fu);
129: DMDAVecGetArray(user->da1,vlambda,&lambda);
130: DMDAVecGetArray(user->da1,vflambda,&flambda);
131: d = (N-1.0);
132: h = 1.0/d;
134: /* derivative of L() w.r.t. w */
135: if (xs == 0) { /* only first processor computes this */
136: fw[0] = -2.*d*lambda[0];
137: }
139: /* derivative of L() w.r.t. u */
140: for (i=xs; i<xs+xm; i++) {
141: if (i == 0) flambda[0] = h*u[0] + 2.*d*lambda[0] - d*lambda[1];
142: else if (i == 1) flambda[1] = 2.*h*u[1] + 2.*d*lambda[1] - d*lambda[2];
143: else if (i == N-1) flambda[N-1] = h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2];
144: else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3];
145: else flambda[i] = 2.*h*u[i] - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]);
146: }
148: /* derivative of L() w.r.t. lambda */
149: for (i=xs; i<xs+xm; i++) {
150: if (i == 0) fu[0] = 2.0*d*(u[0] - w[0]);
151: else if (i == N-1) fu[N-1] = 2.0*d*u[N-1];
152: else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);
153: }
155: VecRestoreArray(vw,&w);
156: VecRestoreArray(vfw,&fw);
157: DMDAVecRestoreArray(user->da1,vu,&u);
158: DMDAVecRestoreArray(user->da1,vfu,&fu);
159: DMDAVecRestoreArray(user->da1,vlambda,&lambda);
160: DMDAVecRestoreArray(user->da1,vflambda,&flambda);
162: DMCompositeGather(user->packer,FU,INSERT_VALUES,vfw,vfu,vflambda);
163: DMCompositeRestoreLocalVectors(user->packer,&vw,&vu,&vlambda);
164: DMCompositeRestoreLocalVectors(user->packer,&vfw,&vfu,&vflambda);
165: return(0);
166: }
170: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
171: {
172: UserCtx *user = (UserCtx*)dummy;
174: Vec w,u,lambda,U,F;
177: SNESGetSolution(snes,&U);
178: DMCompositeGetAccess(user->packer,U,&w,&u,&lambda);
179: VecView(u,user->u_viewer);
180: VecView(lambda,user->lambda_viewer);
181: DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda);
183: SNESGetFunction(snes,&F,0,0);
184: DMCompositeGetAccess(user->packer,F,&w,&u,&lambda);
185: VecView(u,user->fu_viewer);
186: VecView(lambda,user->flambda_viewer);
187: DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda);
188: return(0);
189: }