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;
57: PetscFunctionBeginUser;
58: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
60: /* Create a global vector that includes a single redundant array and two da arrays */
61: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.packer));
62: PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 0, 1, &user.red1));
63: PetscCall(DMCompositeAddDM(user.packer, user.red1));
64: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da1));
65: PetscCall(DMSetFromOptions(user.da1));
66: PetscCall(DMSetUp(user.da1));
67: PetscCall(DMCompositeAddDM(user.packer, user.da1));
68: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da2));
69: PetscCall(DMSetFromOptions(user.da2));
70: PetscCall(DMSetUp(user.da2));
71: PetscCall(DMDASetFieldName(user.da1, 0, "u"));
72: PetscCall(DMDASetFieldName(user.da2, 0, "lambda"));
73: PetscCall(DMCompositeAddDM(user.packer, user.da2));
74: PetscCall(DMCreateGlobalVector(user.packer, &U));
75: PetscCall(VecDuplicate(U, &FU));
77: /* create graphics windows */
78: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "u - state variables", -1, -1, -1, -1, &user.u_viewer));
79: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "lambda - Lagrange multipliers", -1, -1, -1, -1, &user.lambda_viewer));
80: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "fu - derivative w.r.t. state variables", -1, -1, -1, -1, &user.fu_viewer));
81: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "flambda - derivative w.r.t. Lagrange multipliers", -1, -1, -1, -1, &user.flambda_viewer));
83: /* create nonlinear solver */
84: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
85: PetscCall(SNESSetFunction(snes, FU, FormFunction, &user));
86: PetscCall(SNESSetFromOptions(snes));
87: PetscCall(SNESMonitorSet(snes, Monitor, &user, 0));
88: PetscCall(SNESSolve(snes, NULL, U));
89: PetscCall(SNESGetIterationNumber(snes, &its));
90: PetscCall(SNESDestroy(&snes));
92: PetscCall(DMDestroy(&user.red1));
93: PetscCall(DMDestroy(&user.da1));
94: PetscCall(DMDestroy(&user.da2));
95: PetscCall(DMDestroy(&user.packer));
96: PetscCall(VecDestroy(&U));
97: PetscCall(VecDestroy(&FU));
98: PetscCall(PetscViewerDestroy(&user.u_viewer));
99: PetscCall(PetscViewerDestroy(&user.lambda_viewer));
100: PetscCall(PetscViewerDestroy(&user.fu_viewer));
101: PetscCall(PetscViewerDestroy(&user.flambda_viewer));
102: PetscCall(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;
117: PetscFunctionBeginUser;
118: PetscCall(DMCompositeGetLocalVectors(user->packer, &vw, &vu, &vlambda));
119: PetscCall(DMCompositeGetLocalVectors(user->packer, &vfw, &vfu, &vflambda));
120: PetscCall(DMCompositeScatter(user->packer, U, vw, vu, vlambda));
122: PetscCall(DMDAGetCorners(user->da1, &xs, NULL, NULL, &xm, NULL, NULL));
123: PetscCall(DMDAGetInfo(user->da1, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
124: PetscCall(VecGetArray(vw, &w));
125: PetscCall(VecGetArray(vfw, &fw));
126: PetscCall(DMDAVecGetArray(user->da1, vu, &u));
127: PetscCall(DMDAVecGetArray(user->da1, vfu, &fu));
128: PetscCall(DMDAVecGetArray(user->da1, vlambda, &lambda));
129: PetscCall(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: PetscCall(VecRestoreArray(vw, &w));
155: PetscCall(VecRestoreArray(vfw, &fw));
156: PetscCall(DMDAVecRestoreArray(user->da1, vu, &u));
157: PetscCall(DMDAVecRestoreArray(user->da1, vfu, &fu));
158: PetscCall(DMDAVecRestoreArray(user->da1, vlambda, &lambda));
159: PetscCall(DMDAVecRestoreArray(user->da1, vflambda, &flambda));
161: PetscCall(DMCompositeGather(user->packer, INSERT_VALUES, FU, vfw, vfu, vflambda));
162: PetscCall(DMCompositeRestoreLocalVectors(user->packer, &vw, &vu, &vlambda));
163: PetscCall(DMCompositeRestoreLocalVectors(user->packer, &vfw, &vfu, &vflambda));
164: PetscFunctionReturn(PETSC_SUCCESS);
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;
172: PetscFunctionBeginUser;
173: PetscCall(SNESGetSolution(snes, &U));
174: PetscCall(DMCompositeGetAccess(user->packer, U, &w, &u, &lambda));
175: PetscCall(VecView(u, user->u_viewer));
176: PetscCall(VecView(lambda, user->lambda_viewer));
177: PetscCall(DMCompositeRestoreAccess(user->packer, U, &w, &u, &lambda));
179: PetscCall(SNESGetFunction(snes, &F, 0, 0));
180: PetscCall(DMCompositeGetAccess(user->packer, F, &w, &u, &lambda));
181: PetscCall(VecView(u, user->fu_viewer));
182: PetscCall(VecView(lambda, user->flambda_viewer));
183: PetscCall(DMCompositeRestoreAccess(user->packer, F, &w, &u, &lambda));
184: PetscFunctionReturn(PETSC_SUCCESS);
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*/