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