1
2 static const char help[] = "Solves PDE optimization problem using full-space method, interlaces state and adjoint variables.\n\n";
3
4
5 #include <petscdm.h>
6 #include <petscdmda.h>
7 #include <petscdmredundant.h>
8 #include <petscdmcomposite.h>
9 #include <petscpf.h>
10 #include <petscsnes.h>
11 #include <petsc/private/dmimpl.h>
12
13 /*
14
15 w - design variables (what we change to get an optimal solution)
16 u - state variables (i.e. the PDE solution)
17 lambda - the Lagrange multipliers
18
19 U = (w [u_0 lambda_0 u_1 lambda_1 .....])
20
21 fu, fw, flambda contain the gradient of L(w,u,lambda)
22
23 FU = (fw [fu_0 flambda_0 .....])
24
25 In this example the PDE is
26 Uxx = 2,
27 u(0) = w(0), thus this is the free parameter
28 u(1) = 0
29 the function we wish to minimize is
30 \integral u^{2}
31
32 The exact solution for u is given by u(x) = x*x - 1.25*x + .25
33
34 Use the usual centered finite differences.
35
36 Note we treat the problem as non-linear though it happens to be linear
37
38 See ex21.c for the same code, but that does NOT interlaces the u and the lambda
39
40 The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
41 */
42
43 typedef struct {
44 PetscViewer u_lambda_viewer;
45 PetscViewer fu_lambda_viewer;
46 } UserCtx;
47
48 extern PetscErrorCode ComputeFunction(SNES,Vec,Vec,void*);
49 extern PetscErrorCode ComputeJacobian_MF(SNES,Vec,Mat,Mat,void*);
50 extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
51
52 /*
53 Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the
54 smoother on all levels. This is because (1) in the matrix free case no matrix entries are
55 available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal
56 entry for the control variable is zero which means default SOR will not work.
57
58 */
59 char common_options[] = "-ksp_type fgmres\
60 -snes_grid_sequence 2 \
61 -pc_type mg\
62 -mg_levels_pc_type none \
63 -mg_coarse_pc_type none \
64 -pc_mg_type full \
65 -mg_coarse_ksp_type gmres \
66 -mg_levels_ksp_type gmres \
67 -mg_coarse_ksp_max_it 6 \
68 -mg_levels_ksp_max_it 3";
69
70 char matrix_free_options[] = "-mat_mffd_compute_normu no \
71 -mat_mffd_type wp";
72
73 extern PetscErrorCode DMCreateMatrix_MF(DM,Mat*);
74
main(int argc,char ** argv)75 int main(int argc,char **argv)
76 {
77 PetscErrorCode ierr;
78 UserCtx user;
79 DM red,da;
80 SNES snes;
81 DM packer;
82 PetscBool use_monitor = PETSC_FALSE;
83
84 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
85
86 /* Hardwire several options; can be changed at command line */
87 ierr = PetscOptionsInsertString(NULL,common_options);CHKERRQ(ierr);
88 ierr = PetscOptionsInsertString(NULL,matrix_free_options);CHKERRQ(ierr);
89 ierr = PetscOptionsInsert(NULL,&argc,&argv,NULL);CHKERRQ(ierr);
90 ierr = PetscOptionsGetBool(NULL,NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);CHKERRQ(ierr);
91
92 /* Create a global vector that includes a single redundant array and two da arrays */
93 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&packer);CHKERRQ(ierr);
94 ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);CHKERRQ(ierr);
95 ierr = DMSetOptionsPrefix(red,"red_");CHKERRQ(ierr);
96 ierr = DMCompositeAddDM(packer,red);CHKERRQ(ierr);
97 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,2,1,NULL,&da);CHKERRQ(ierr);
98 ierr = DMSetOptionsPrefix(red,"da_");CHKERRQ(ierr);
99 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
100 ierr = DMSetUp(da);CHKERRQ(ierr);
101 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
102 ierr = DMDASetFieldName(da,1,"lambda");CHKERRQ(ierr);
103 ierr = DMCompositeAddDM(packer,(DM)da);CHKERRQ(ierr);
104 ierr = DMSetApplicationContext(packer,&user);CHKERRQ(ierr);
105
106 packer->ops->creatematrix = DMCreateMatrix_MF;
107
108 /* create nonlinear multi-level solver */
109 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
110 ierr = SNESSetDM(snes,packer);CHKERRQ(ierr);
111 ierr = SNESSetFunction(snes,NULL,ComputeFunction,NULL);CHKERRQ(ierr);
112 ierr = SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL);CHKERRQ(ierr);
113
114 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
115
116 if (use_monitor) {
117 /* create graphics windows */
118 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);CHKERRQ(ierr);
119 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);CHKERRQ(ierr);
120 ierr = SNESMonitorSet(snes,Monitor,0,0);CHKERRQ(ierr);
121 }
122
123 ierr = SNESSolve(snes,NULL,NULL);CHKERRQ(ierr);
124 ierr = SNESDestroy(&snes);CHKERRQ(ierr);
125
126 ierr = DMDestroy(&red);CHKERRQ(ierr);
127 ierr = DMDestroy(&da);CHKERRQ(ierr);
128 ierr = DMDestroy(&packer);CHKERRQ(ierr);
129 if (use_monitor) {
130 ierr = PetscViewerDestroy(&user.u_lambda_viewer);CHKERRQ(ierr);
131 ierr = PetscViewerDestroy(&user.fu_lambda_viewer);CHKERRQ(ierr);
132 }
133 ierr = PetscFinalize();
134 return ierr;
135 }
136
137 typedef struct {
138 PetscScalar u;
139 PetscScalar lambda;
140 } ULambda;
141
142 /*
143 Evaluates FU = Gradiant(L(w,u,lambda))
144
145 This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and
146 DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).
147
148 */
ComputeFunction(SNES snes,Vec U,Vec FU,void * ctx)149 PetscErrorCode ComputeFunction(SNES snes,Vec U,Vec FU,void *ctx)
150 {
151 PetscErrorCode ierr;
152 PetscInt xs,xm,i,N;
153 ULambda *u_lambda,*fu_lambda;
154 PetscScalar d,h,*w,*fw;
155 Vec vw,vfw,vu_lambda,vfu_lambda;
156 DM packer,red,da;
157
158 PetscFunctionBeginUser;
159 ierr = VecGetDM(U, &packer);CHKERRQ(ierr);
160 ierr = DMCompositeGetEntries(packer,&red,&da);CHKERRQ(ierr);
161 ierr = DMCompositeGetLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr);
162 ierr = DMCompositeScatter(packer,U,vw,vu_lambda);CHKERRQ(ierr);
163 ierr = DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr);
164
165 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
166 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
167 ierr = VecGetArray(vw,&w);CHKERRQ(ierr);
168 ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr);
169 ierr = DMDAVecGetArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr);
170 ierr = DMDAVecGetArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr);
171 d = N-1.0;
172 h = 1.0/d;
173
174 /* derivative of L() w.r.t. w */
175 if (xs == 0) { /* only first processor computes this */
176 fw[0] = -2.0*d*u_lambda[0].lambda;
177 }
178
179 /* derivative of L() w.r.t. u */
180 for (i=xs; i<xs+xm; i++) {
181 if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda;
182 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;
183 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;
184 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;
185 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);
186 }
187
188 /* derivative of L() w.r.t. lambda */
189 for (i=xs; i<xs+xm; i++) {
190 if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]);
191 else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;
192 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);
193 }
194
195 ierr = VecRestoreArray(vw,&w);CHKERRQ(ierr);
196 ierr = VecRestoreArray(vfw,&fw);CHKERRQ(ierr);
197 ierr = DMDAVecRestoreArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr);
198 ierr = DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr);
199 ierr = DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr);
200 ierr = DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr);
201 ierr = PetscLogFlops(13.0*N);CHKERRQ(ierr);
202 PetscFunctionReturn(0);
203 }
204
205 /*
206 Computes the exact solution
207 */
u_solution(void * dummy,PetscInt n,const PetscScalar * x,PetscScalar * u)208 PetscErrorCode u_solution(void *dummy,PetscInt n,const PetscScalar *x,PetscScalar *u)
209 {
210 PetscInt i;
211
212 PetscFunctionBeginUser;
213 for (i=0; i<n; i++) u[2*i] = x[i]*x[i] - 1.25*x[i] + .25;
214 PetscFunctionReturn(0);
215 }
216
ExactSolution(DM packer,Vec U)217 PetscErrorCode ExactSolution(DM packer,Vec U)
218 {
219 PF pf;
220 Vec x,u_global;
221 PetscScalar *w;
222 DM da;
223 PetscErrorCode ierr;
224 PetscInt m;
225
226 PetscFunctionBeginUser;
227 ierr = DMCompositeGetEntries(packer,&m,&da);CHKERRQ(ierr);
228
229 ierr = PFCreate(PETSC_COMM_WORLD,1,2,&pf);CHKERRQ(ierr);
230 /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */
231 ierr = PFSetType(pf,PFQUICK,(void*)(PETSC_UINTPTR_T)u_solution);CHKERRQ(ierr);
232 ierr = DMGetCoordinates(da,&x);CHKERRQ(ierr);
233 if (!x) {
234 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
235 ierr = DMGetCoordinates(da,&x);CHKERRQ(ierr);
236 }
237 ierr = DMCompositeGetAccess(packer,U,&w,&u_global,0);CHKERRQ(ierr);
238 if (w) w[0] = .25;
239 ierr = PFApplyVec(pf,x,u_global);CHKERRQ(ierr);
240 ierr = PFDestroy(&pf);CHKERRQ(ierr);
241 ierr = DMCompositeRestoreAccess(packer,U,&w,&u_global,0);CHKERRQ(ierr);
242 PetscFunctionReturn(0);
243 }
244
Monitor(SNES snes,PetscInt its,PetscReal rnorm,void * dummy)245 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
246 {
247 UserCtx *user;
248 PetscErrorCode ierr;
249 PetscInt m,N;
250 PetscScalar *w,*dw;
251 Vec u_lambda,U,F,Uexact;
252 DM packer;
253 PetscReal norm;
254 DM da;
255
256 PetscFunctionBeginUser;
257 ierr = SNESGetDM(snes,&packer);CHKERRQ(ierr);
258 ierr = DMGetApplicationContext(packer,&user);CHKERRQ(ierr);
259 ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr);
260 ierr = DMCompositeGetAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr);
261 ierr = VecView(u_lambda,user->u_lambda_viewer);CHKERRQ(ierr);
262 ierr = DMCompositeRestoreAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr);
263
264 ierr = SNESGetFunction(snes,&F,0,0);CHKERRQ(ierr);
265 ierr = DMCompositeGetAccess(packer,F,&w,&u_lambda);CHKERRQ(ierr);
266 /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */
267 ierr = DMCompositeRestoreAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr);
268
269 ierr = DMCompositeGetEntries(packer,&m,&da);CHKERRQ(ierr);
270 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
271 ierr = VecDuplicate(U,&Uexact);CHKERRQ(ierr);
272 ierr = ExactSolution(packer,Uexact);CHKERRQ(ierr);
273 ierr = VecAXPY(Uexact,-1.0,U);CHKERRQ(ierr);
274 ierr = DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);CHKERRQ(ierr);
275 ierr = VecStrideNorm(u_lambda,0,NORM_2,&norm);CHKERRQ(ierr);
276 norm = norm/PetscSqrtReal((PetscReal)N-1.);
277 if (dw) ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Error at x = 0 %g\n",(double)norm,(double)PetscRealPart(dw[0]));CHKERRQ(ierr);
278 ierr = VecView(u_lambda,user->fu_lambda_viewer);CHKERRQ(ierr);
279 ierr = DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);CHKERRQ(ierr);
280 ierr = VecDestroy(&Uexact);CHKERRQ(ierr);
281 PetscFunctionReturn(0);
282 }
283
DMCreateMatrix_MF(DM packer,Mat * A)284 PetscErrorCode DMCreateMatrix_MF(DM packer,Mat *A)
285 {
286 PetscErrorCode ierr;
287 Vec t;
288 PetscInt m;
289
290 PetscFunctionBeginUser;
291 ierr = DMGetGlobalVector(packer,&t);CHKERRQ(ierr);
292 ierr = VecGetLocalSize(t,&m);CHKERRQ(ierr);
293 ierr = DMRestoreGlobalVector(packer,&t);CHKERRQ(ierr);
294 ierr = MatCreateMFFD(PETSC_COMM_WORLD,m,m,PETSC_DETERMINE,PETSC_DETERMINE,A);CHKERRQ(ierr);
295 ierr = MatSetUp(*A);CHKERRQ(ierr);
296 ierr = MatSetDM(*A,packer);CHKERRQ(ierr);
297 PetscFunctionReturn(0);
298 }
299
ComputeJacobian_MF(SNES snes,Vec x,Mat A,Mat B,void * ctx)300 PetscErrorCode ComputeJacobian_MF(SNES snes,Vec x,Mat A,Mat B,void *ctx)
301 {
302 PetscErrorCode ierr;
303
304 PetscFunctionBeginUser;
305 ierr = MatMFFDSetFunction(A,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr);
306 ierr = MatMFFDSetBase(A,x,NULL);CHKERRQ(ierr);
307 PetscFunctionReturn(0);
308 }
309
310
311 /*TEST
312
313 test:
314 nsize: 2
315 args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view
316 requires: !single
317
318 TEST*/
319