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, 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;
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: PetscFunctionBeginUser;
208: for (PetscInt i = 0; i < n; i++) u[2 * i] = x[i] * x[i] - 1.25 * x[i] + .25;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: PetscErrorCode ExactSolution(DM packer, Vec U)
213: {
214: PF pf;
215: Vec x, u_global;
216: PetscScalar *w;
217: DM da;
218: PetscInt m;
220: PetscFunctionBeginUser;
221: PetscCall(DMCompositeGetEntries(packer, &m, &da));
223: PetscCall(PFCreate(PETSC_COMM_WORLD, 1, 2, &pf));
224: /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */
225: PetscCall(PFSetType(pf, PFQUICK, (void *)(PETSC_UINTPTR_T)u_solution));
226: PetscCall(DMGetCoordinates(da, &x));
227: if (!x) {
228: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
229: PetscCall(DMGetCoordinates(da, &x));
230: }
231: PetscCall(DMCompositeGetAccess(packer, U, &w, &u_global, 0));
232: if (w) w[0] = .25;
233: PetscCall(PFApplyVec(pf, x, u_global));
234: PetscCall(PFDestroy(&pf));
235: PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_global, 0));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal rnorm, void *dummy)
240: {
241: UserCtx *user;
242: PetscInt m, N;
243: PetscScalar *w, *dw;
244: Vec u_lambda, U, F, Uexact;
245: DM packer;
246: PetscReal norm;
247: DM da;
249: PetscFunctionBeginUser;
250: PetscCall(SNESGetDM(snes, &packer));
251: PetscCall(DMGetApplicationContext(packer, &user));
252: PetscCall(SNESGetSolution(snes, &U));
253: PetscCall(DMCompositeGetAccess(packer, U, &w, &u_lambda));
254: PetscCall(VecView(u_lambda, user->u_lambda_viewer));
255: PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_lambda));
257: PetscCall(SNESGetFunction(snes, &F, 0, 0));
258: PetscCall(DMCompositeGetAccess(packer, F, &w, &u_lambda));
259: /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */
260: PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_lambda));
262: PetscCall(DMCompositeGetEntries(packer, &m, &da));
263: PetscCall(DMDAGetInfo(da, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
264: PetscCall(VecDuplicate(U, &Uexact));
265: PetscCall(ExactSolution(packer, Uexact));
266: PetscCall(VecAXPY(Uexact, -1.0, U));
267: PetscCall(DMCompositeGetAccess(packer, Uexact, &dw, &u_lambda));
268: PetscCall(VecStrideNorm(u_lambda, 0, NORM_2, &norm));
269: norm = norm / PetscSqrtReal((PetscReal)N - 1.);
270: if (dw) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Error at x = 0 %g\n", (double)norm, (double)PetscRealPart(dw[0])));
271: PetscCall(VecView(u_lambda, user->fu_lambda_viewer));
272: PetscCall(DMCompositeRestoreAccess(packer, Uexact, &dw, &u_lambda));
273: PetscCall(VecDestroy(&Uexact));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode DMCreateMatrix_MF(DM packer, Mat *A)
278: {
279: Vec t;
280: PetscInt m;
282: PetscFunctionBeginUser;
283: PetscCall(DMGetGlobalVector(packer, &t));
284: PetscCall(VecGetLocalSize(t, &m));
285: PetscCall(DMRestoreGlobalVector(packer, &t));
286: PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, A));
287: PetscCall(MatSetUp(*A));
288: PetscCall(MatSetDM(*A, packer));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: PetscErrorCode ComputeJacobian_MF(SNES snes, Vec x, Mat A, Mat B, PetscCtx ctx)
293: {
294: PetscFunctionBeginUser;
295: PetscCall(MatMFFDSetFunction(A, (PetscErrorCode (*)(void *, Vec, Vec))SNESComputeFunction, snes));
296: PetscCall(MatMFFDSetBase(A, x, NULL));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*TEST
302: test:
303: nsize: 2
304: args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view
305: requires: !single
307: TEST*/