1 static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmda.h>
5 #include <petscdmredundant.h>
6 #include <petscdmcomposite.h>
7 #include <petscpf.h>
8 #include <petscsnes.h>
9
10 /*
11
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
15
16 U = (w u lambda)
17
18 fu, fw, flambda contain the gradient of L(w,u,lambda)
19
20 FU = (fw fu flambda)
21
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}
28
29 The exact solution for u is given by u(x) = x*x - 1.25*x + .25
30
31 Use the usual centered finite differences.
32
33 Note we treat the problem as non-linear though it happens to be linear
34
35 See ex22.c for the same code, but that interlaces the u and the lambda
36
37 */
38
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;
45
46 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
47 extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
48
main(int argc,char ** argv)49 int main(int argc, char **argv)
50 {
51 PetscInt its;
52 Vec U, FU;
53 SNES snes;
54 UserCtx user;
55
56 PetscFunctionBeginUser;
57 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
58
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));
75
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));
81
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));
90
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 }
104
105 /*
106 Evaluates FU = Gradient(L(w,u,lambda))
107
108 */
FormFunction(SNES snes,Vec U,Vec FU,void * dummy)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;
115
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));
120
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;
131
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 }
136
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 }
145
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 }
152
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));
159
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 }
165
Monitor(SNES snes,PetscInt its,PetscReal rnorm,void * dummy)166 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal rnorm, void *dummy)
167 {
168 UserCtx *user = (UserCtx *)dummy;
169 Vec w, u, lambda, U, F;
170
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));
177
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 }
185
186 /*TEST
187
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
192
193 TEST*/
194