Actual source code: ex21.c
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 *);
50: int main(int argc, char **argv)
51: {
52: PetscInt its;
53: Vec U, FU;
54: SNES snes;
55: UserCtx user;
58: PetscInitialize(&argc, &argv, (char *)0, help);
60: /* Create a global vector that includes a single redundant array and two da arrays */
61: DMCompositeCreate(PETSC_COMM_WORLD, &user.packer);
62: DMRedundantCreate(PETSC_COMM_WORLD, 0, 1, &user.red1);
63: DMCompositeAddDM(user.packer, user.red1);
64: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da1);
65: DMSetFromOptions(user.da1);
66: DMSetUp(user.da1);
67: DMCompositeAddDM(user.packer, user.da1);
68: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da2);
69: DMSetFromOptions(user.da2);
70: DMSetUp(user.da2);
71: DMDASetFieldName(user.da1, 0, "u");
72: DMDASetFieldName(user.da2, 0, "lambda");
73: DMCompositeAddDM(user.packer, user.da2);
74: DMCreateGlobalVector(user.packer, &U);
75: VecDuplicate(U, &FU);
77: /* create graphics windows */
78: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "u - state variables", -1, -1, -1, -1, &user.u_viewer);
79: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "lambda - Lagrange multipliers", -1, -1, -1, -1, &user.lambda_viewer);
80: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "fu - derivate w.r.t. state variables", -1, -1, -1, -1, &user.fu_viewer);
81: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "flambda - derivate w.r.t. Lagrange multipliers", -1, -1, -1, -1, &user.flambda_viewer);
83: /* create nonlinear solver */
84: SNESCreate(PETSC_COMM_WORLD, &snes);
85: SNESSetFunction(snes, FU, FormFunction, &user);
86: SNESSetFromOptions(snes);
87: SNESMonitorSet(snes, Monitor, &user, 0);
88: SNESSolve(snes, NULL, U);
89: SNESGetIterationNumber(snes, &its);
90: SNESDestroy(&snes);
92: DMDestroy(&user.red1);
93: DMDestroy(&user.da1);
94: DMDestroy(&user.da2);
95: DMDestroy(&user.packer);
96: VecDestroy(&U);
97: VecDestroy(&FU);
98: PetscViewerDestroy(&user.u_viewer);
99: PetscViewerDestroy(&user.lambda_viewer);
100: PetscViewerDestroy(&user.fu_viewer);
101: PetscViewerDestroy(&user.flambda_viewer);
102: PetscFinalize();
103: return 0;
104: }
106: /*
107: Evaluates FU = Gradient(L(w,u,lambda))
109: */
110: PetscErrorCode FormFunction(SNES snes, Vec U, Vec FU, void *dummy)
111: {
112: UserCtx *user = (UserCtx *)dummy;
113: PetscInt xs, xm, i, N;
114: PetscScalar *u, *lambda, *w, *fu, *fw, *flambda, d, h;
115: Vec vw, vu, vlambda, vfw, vfu, vflambda;
118: DMCompositeGetLocalVectors(user->packer, &vw, &vu, &vlambda);
119: DMCompositeGetLocalVectors(user->packer, &vfw, &vfu, &vflambda);
120: DMCompositeScatter(user->packer, U, vw, vu, vlambda);
122: DMDAGetCorners(user->da1, &xs, NULL, NULL, &xm, NULL, NULL);
123: DMDAGetInfo(user->da1, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
124: VecGetArray(vw, &w);
125: VecGetArray(vfw, &fw);
126: DMDAVecGetArray(user->da1, vu, &u);
127: DMDAVecGetArray(user->da1, vfu, &fu);
128: DMDAVecGetArray(user->da1, vlambda, &lambda);
129: DMDAVecGetArray(user->da1, vflambda, &flambda);
130: d = (N - 1.0);
131: h = 1.0 / d;
133: /* derivative of L() w.r.t. w */
134: if (xs == 0) { /* only first processor computes this */
135: fw[0] = -2. * d * lambda[0];
136: }
138: /* derivative of L() w.r.t. u */
139: for (i = xs; i < xs + xm; i++) {
140: if (i == 0) flambda[0] = h * u[0] + 2. * d * lambda[0] - d * lambda[1];
141: else if (i == 1) flambda[1] = 2. * h * u[1] + 2. * d * lambda[1] - d * lambda[2];
142: else if (i == N - 1) flambda[N - 1] = h * u[N - 1] + 2. * d * lambda[N - 1] - d * lambda[N - 2];
143: else if (i == N - 2) flambda[N - 2] = 2. * h * u[N - 2] + 2. * d * lambda[N - 2] - d * lambda[N - 3];
144: else flambda[i] = 2. * h * u[i] - d * (lambda[i + 1] - 2.0 * lambda[i] + lambda[i - 1]);
145: }
147: /* derivative of L() w.r.t. lambda */
148: for (i = xs; i < xs + xm; i++) {
149: if (i == 0) fu[0] = 2.0 * d * (u[0] - w[0]);
150: else if (i == N - 1) fu[N - 1] = 2.0 * d * u[N - 1];
151: else fu[i] = -(d * (u[i + 1] - 2.0 * u[i] + u[i - 1]) - 2.0 * h);
152: }
154: VecRestoreArray(vw, &w);
155: VecRestoreArray(vfw, &fw);
156: DMDAVecRestoreArray(user->da1, vu, &u);
157: DMDAVecRestoreArray(user->da1, vfu, &fu);
158: DMDAVecRestoreArray(user->da1, vlambda, &lambda);
159: DMDAVecRestoreArray(user->da1, vflambda, &flambda);
161: DMCompositeGather(user->packer, INSERT_VALUES, FU, vfw, vfu, vflambda);
162: DMCompositeRestoreLocalVectors(user->packer, &vw, &vu, &vlambda);
163: DMCompositeRestoreLocalVectors(user->packer, &vfw, &vfu, &vflambda);
164: return 0;
165: }
167: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal rnorm, void *dummy)
168: {
169: UserCtx *user = (UserCtx *)dummy;
170: Vec w, u, lambda, U, F;
173: SNESGetSolution(snes, &U);
174: DMCompositeGetAccess(user->packer, U, &w, &u, &lambda);
175: VecView(u, user->u_viewer);
176: VecView(lambda, user->lambda_viewer);
177: DMCompositeRestoreAccess(user->packer, U, &w, &u, &lambda);
179: SNESGetFunction(snes, &F, 0, 0);
180: DMCompositeGetAccess(user->packer, F, &w, &u, &lambda);
181: VecView(u, user->fu_viewer);
182: VecView(lambda, user->flambda_viewer);
183: DMCompositeRestoreAccess(user->packer, F, &w, &u, &lambda);
184: return 0;
185: }
187: /*TEST
189: test:
190: nsize: 4
191: args: -snes_linesearch_monitor -snes_mf -pc_type none -snes_monitor_short -nox -ksp_monitor_short -snes_converged_reason
192: requires: !single
194: TEST*/