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