Actual source code: ex22.c
1: static const char help[] = "Solves PDE optimization problem using full-space method, interlaces state and adjoint variables.\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>
9: #include <petsc/private/dmimpl.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_0 lambda_0 u_1 lambda_1 .....])
19: fu, fw, flambda contain the gradient of L(w,u,lambda)
21: FU = (fw [fu_0 flambda_0 .....])
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 ex21.c for the same code, but that does NOT interlaces the u and the lambda
38: The vectors u_lambda and fu_lambda contain the u and the lambda interlaced
39: */
41: typedef struct {
42: PetscViewer u_lambda_viewer;
43: PetscViewer fu_lambda_viewer;
44: } UserCtx;
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 *);
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.
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";
68: char matrix_free_options[] = "-mat_mffd_compute_normu no \
69: -mat_mffd_type wp";
71: extern PetscErrorCode DMCreateMatrix_MF(DM, Mat *);
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;
81: PetscFunctionBeginUser;
82: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
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));
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));
104: packer->ops->creatematrix = DMCreateMatrix_MF;
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));
112: PetscCall(SNESSetFromOptions(snes));
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: }
121: PetscCall(SNESSolve(snes, NULL, NULL));
122: PetscCall(SNESDestroy(&snes));
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: }
135: typedef struct {
136: PetscScalar u;
137: PetscScalar lambda;
138: } ULambda;
140: /*
141: Evaluates FU = Gradient(L(w,u,lambda))
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()).
146: */
147: PetscErrorCode ComputeFunction(SNES snes, Vec U, Vec FU, void *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;
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));
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;
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: }
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: }
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: }
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: }
202: /*
203: Computes the exact solution
204: */
205: PetscErrorCode u_solution(void *dummy, PetscInt n, const PetscScalar *x, PetscScalar *u)
206: {
207: PetscInt i;
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: }
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;
222: PetscFunctionBeginUser;
223: PetscCall(DMCompositeGetEntries(packer, &m, &da));
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: }
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;
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));
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));
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: }
279: PetscErrorCode DMCreateMatrix_MF(DM packer, Mat *A)
280: {
281: Vec t;
282: PetscInt m;
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: }
294: PetscErrorCode ComputeJacobian_MF(SNES snes, Vec x, Mat A, Mat B, void *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: }
302: /*TEST
304: test:
305: nsize: 2
306: args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view
307: requires: !single
309: TEST*/