Actual source code: ex21.c
1: static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
5: #include <petscdmredundant.h>
6: #include <petscdmcomposite.h>
7: #include <petscpf.h>
8: #include <petscsnes.h>
10: /*
12: w - design variables (what we change to get an optimal solution)
13: u - state variables (i.e. the PDE solution)
14: lambda - the Lagrange multipliers
16: U = (w u lambda)
18: fu, fw, flambda contain the gradient of L(w,u,lambda)
20: FU = (fw fu flambda)
22: In this example the PDE is
23: Uxx = 2,
24: u(0) = w(0), thus this is the free parameter
25: u(1) = 0
26: the function we wish to minimize is
27: \integral u^{2}
29: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
31: Use the usual centered finite differences.
33: Note we treat the problem as non-linear though it happens to be linear
35: See ex22.c for the same code, but that interlaces the u and the lambda
37: */
39: typedef struct {
40: DM red1, da1, da2;
41: DM packer;
42: PetscViewer u_viewer, lambda_viewer;
43: PetscViewer fu_viewer, flambda_viewer;
44: } UserCtx;
46: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
47: extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
49: int main(int argc, char **argv)
50: {
51: PetscInt its;
52: Vec U, FU;
53: SNES snes;
54: UserCtx user;
56: PetscFunctionBeginUser;
57: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
59: /* Create a global vector that includes a single redundant array and two da arrays */
60: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.packer));
61: PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 0, 1, &user.red1));
62: PetscCall(DMCompositeAddDM(user.packer, user.red1));
63: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da1));
64: PetscCall(DMSetFromOptions(user.da1));
65: PetscCall(DMSetUp(user.da1));
66: PetscCall(DMCompositeAddDM(user.packer, user.da1));
67: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da2));
68: PetscCall(DMSetFromOptions(user.da2));
69: PetscCall(DMSetUp(user.da2));
70: PetscCall(DMDASetFieldName(user.da1, 0, "u"));
71: PetscCall(DMDASetFieldName(user.da2, 0, "lambda"));
72: PetscCall(DMCompositeAddDM(user.packer, user.da2));
73: PetscCall(DMCreateGlobalVector(user.packer, &U));
74: PetscCall(VecDuplicate(U, &FU));
76: /* create graphics windows */
77: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "u - state variables", -1, -1, -1, -1, &user.u_viewer));
78: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "lambda - Lagrange multipliers", -1, -1, -1, -1, &user.lambda_viewer));
79: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "fu - derivative w.r.t. state variables", -1, -1, -1, -1, &user.fu_viewer));
80: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "flambda - derivative w.r.t. Lagrange multipliers", -1, -1, -1, -1, &user.flambda_viewer));
82: /* create nonlinear solver */
83: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
84: PetscCall(SNESSetFunction(snes, FU, FormFunction, &user));
85: PetscCall(SNESSetFromOptions(snes));
86: PetscCall(SNESMonitorSet(snes, Monitor, &user, 0));
87: PetscCall(SNESSolve(snes, NULL, U));
88: PetscCall(SNESGetIterationNumber(snes, &its));
89: PetscCall(SNESDestroy(&snes));
91: PetscCall(DMDestroy(&user.red1));
92: PetscCall(DMDestroy(&user.da1));
93: PetscCall(DMDestroy(&user.da2));
94: PetscCall(DMDestroy(&user.packer));
95: PetscCall(VecDestroy(&U));
96: PetscCall(VecDestroy(&FU));
97: PetscCall(PetscViewerDestroy(&user.u_viewer));
98: PetscCall(PetscViewerDestroy(&user.lambda_viewer));
99: PetscCall(PetscViewerDestroy(&user.fu_viewer));
100: PetscCall(PetscViewerDestroy(&user.flambda_viewer));
101: PetscCall(PetscFinalize());
102: return 0;
103: }
105: /*
106: Evaluates FU = Gradient(L(w,u,lambda))
108: */
109: PetscErrorCode FormFunction(SNES snes, Vec U, Vec FU, void *dummy)
110: {
111: UserCtx *user = (UserCtx *)dummy;
112: PetscInt xs, xm, i, N;
113: PetscScalar *u, *lambda, *w, *fu, *fw, *flambda, d, h;
114: Vec vw, vu, vlambda, vfw, vfu, vflambda;
116: PetscFunctionBeginUser;
117: PetscCall(DMCompositeGetLocalVectors(user->packer, &vw, &vu, &vlambda));
118: PetscCall(DMCompositeGetLocalVectors(user->packer, &vfw, &vfu, &vflambda));
119: PetscCall(DMCompositeScatter(user->packer, U, vw, vu, vlambda));
121: PetscCall(DMDAGetCorners(user->da1, &xs, NULL, NULL, &xm, NULL, NULL));
122: PetscCall(DMDAGetInfo(user->da1, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
123: PetscCall(VecGetArray(vw, &w));
124: PetscCall(VecGetArray(vfw, &fw));
125: PetscCall(DMDAVecGetArray(user->da1, vu, &u));
126: PetscCall(DMDAVecGetArray(user->da1, vfu, &fu));
127: PetscCall(DMDAVecGetArray(user->da1, vlambda, &lambda));
128: PetscCall(DMDAVecGetArray(user->da1, vflambda, &flambda));
129: d = (N - 1.0);
130: h = 1.0 / d;
132: /* derivative of L() w.r.t. w */
133: if (xs == 0) { /* only first processor computes this */
134: fw[0] = -2. * d * lambda[0];
135: }
137: /* derivative of L() w.r.t. u */
138: for (i = xs; i < xs + xm; i++) {
139: if (i == 0) flambda[0] = h * u[0] + 2. * d * lambda[0] - d * lambda[1];
140: else if (i == 1) flambda[1] = 2. * h * u[1] + 2. * d * lambda[1] - d * lambda[2];
141: else if (i == N - 1) flambda[N - 1] = h * u[N - 1] + 2. * d * lambda[N - 1] - d * lambda[N - 2];
142: else if (i == N - 2) flambda[N - 2] = 2. * h * u[N - 2] + 2. * d * lambda[N - 2] - d * lambda[N - 3];
143: else flambda[i] = 2. * h * u[i] - d * (lambda[i + 1] - 2.0 * lambda[i] + lambda[i - 1]);
144: }
146: /* derivative of L() w.r.t. lambda */
147: for (i = xs; i < xs + xm; i++) {
148: if (i == 0) fu[0] = 2.0 * d * (u[0] - w[0]);
149: else if (i == N - 1) fu[N - 1] = 2.0 * d * u[N - 1];
150: else fu[i] = -(d * (u[i + 1] - 2.0 * u[i] + u[i - 1]) - 2.0 * h);
151: }
153: PetscCall(VecRestoreArray(vw, &w));
154: PetscCall(VecRestoreArray(vfw, &fw));
155: PetscCall(DMDAVecRestoreArray(user->da1, vu, &u));
156: PetscCall(DMDAVecRestoreArray(user->da1, vfu, &fu));
157: PetscCall(DMDAVecRestoreArray(user->da1, vlambda, &lambda));
158: PetscCall(DMDAVecRestoreArray(user->da1, vflambda, &flambda));
160: PetscCall(DMCompositeGather(user->packer, INSERT_VALUES, FU, vfw, vfu, vflambda));
161: PetscCall(DMCompositeRestoreLocalVectors(user->packer, &vw, &vu, &vlambda));
162: PetscCall(DMCompositeRestoreLocalVectors(user->packer, &vfw, &vfu, &vflambda));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal rnorm, void *dummy)
167: {
168: UserCtx *user = (UserCtx *)dummy;
169: Vec w, u, lambda, U, F;
171: PetscFunctionBeginUser;
172: PetscCall(SNESGetSolution(snes, &U));
173: PetscCall(DMCompositeGetAccess(user->packer, U, &w, &u, &lambda));
174: PetscCall(VecView(u, user->u_viewer));
175: PetscCall(VecView(lambda, user->lambda_viewer));
176: PetscCall(DMCompositeRestoreAccess(user->packer, U, &w, &u, &lambda));
178: PetscCall(SNESGetFunction(snes, &F, 0, 0));
179: PetscCall(DMCompositeGetAccess(user->packer, F, &w, &u, &lambda));
180: PetscCall(VecView(u, user->fu_viewer));
181: PetscCall(VecView(lambda, user->flambda_viewer));
182: PetscCall(DMCompositeRestoreAccess(user->packer, F, &w, &u, &lambda));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*TEST
188: test:
189: nsize: 4
190: args: -snes_linesearch_monitor -snes_mf -pc_type none -snes_monitor_short -nox -ksp_monitor_short -snes_converged_reason
191: requires: !single
193: TEST*/