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