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