Actual source code: ex21.c
petsc-3.10.5 2019-03-28
2: static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n";
4: /*T
5: requires: !single
6: T*/
8: #include <petscdm.h>
9: #include <petscdmda.h>
10: #include <petscdmredundant.h>
11: #include <petscdmcomposite.h>
12: #include <petscpf.h>
13: #include <petscsnes.h>
15: /*
17: w - design variables (what we change to get an optimal solution)
18: u - state variables (i.e. the PDE solution)
19: lambda - the Lagrange multipliers
21: U = (w u lambda)
23: fu, fw, flambda contain the gradient of L(w,u,lambda)
25: FU = (fw fu flambda)
27: In this example the PDE is
28: Uxx = 2,
29: u(0) = w(0), thus this is the free parameter
30: u(1) = 0
31: the function we wish to minimize is
32: \integral u^{2}
34: The exact solution for u is given by u(x) = x*x - 1.25*x + .25
36: Use the usual centered finite differences.
38: Note we treat the problem as non-linear though it happens to be linear
40: See ex22.c for the same code, but that interlaces the u and the lambda
42: */
44: typedef struct {
45: DM red1,da1,da2;
46: DM packer;
47: PetscViewer u_viewer,lambda_viewer;
48: PetscViewer fu_viewer,flambda_viewer;
49: } UserCtx;
51: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
52: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
55: int main(int argc,char **argv)
56: {
58: PetscInt its;
59: Vec U,FU;
60: SNES snes;
61: UserCtx user;
63: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
65: /* Create a global vector that includes a single redundant array and two da arrays */
66: DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);
67: DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);
68: DMCompositeAddDM(user.packer,user.red1);
69: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,1,1,NULL,&user.da1);
70: DMSetFromOptions(user.da1);
71: DMSetUp(user.da1);
72: DMCompositeAddDM(user.packer,user.da1);
73: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,1,1,NULL,&user.da2);
74: DMSetFromOptions(user.da2);
75: DMSetUp(user.da2);
76: DMDASetFieldName(user.da1,0,"u");
77: DMDASetFieldName(user.da2,0,"lambda");
78: DMCompositeAddDM(user.packer,user.da2);
79: DMCreateGlobalVector(user.packer,&U);
80: VecDuplicate(U,&FU);
82: /* create graphics windows */
83: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
84: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);
85: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);
86: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);
89: /* create nonlinear solver */
90: SNESCreate(PETSC_COMM_WORLD,&snes);
91: SNESSetFunction(snes,FU,FormFunction,&user);
92: SNESSetFromOptions(snes);
93: SNESMonitorSet(snes,Monitor,&user,0);
94: SNESSolve(snes,NULL,U);
95: SNESGetIterationNumber(snes,&its);
96: SNESDestroy(&snes);
98: DMDestroy(&user.red1);
99: DMDestroy(&user.da1);
100: DMDestroy(&user.da2);
101: DMDestroy(&user.packer);
102: VecDestroy(&U);
103: VecDestroy(&FU);
104: PetscViewerDestroy(&user.u_viewer);
105: PetscViewerDestroy(&user.lambda_viewer);
106: PetscViewerDestroy(&user.fu_viewer);
107: PetscViewerDestroy(&user.flambda_viewer);
108: PetscFinalize();
109: return ierr;
110: }
112: /*
113: Evaluates FU = Gradiant(L(w,u,lambda))
115: */
116: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy)
117: {
118: UserCtx *user = (UserCtx*)dummy;
120: PetscInt xs,xm,i,N;
121: PetscScalar *u,*lambda,*w,*fu,*fw,*flambda,d,h;
122: Vec vw,vu,vlambda,vfw,vfu,vflambda;
125: DMCompositeGetLocalVectors(user->packer,&vw,&vu,&vlambda);
126: DMCompositeGetLocalVectors(user->packer,&vfw,&vfu,&vflambda);
127: DMCompositeScatter(user->packer,U,vw,vu,vlambda);
129: DMDAGetCorners(user->da1,&xs,NULL,NULL,&xm,NULL,NULL);
130: DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0);
131: VecGetArray(vw,&w);
132: VecGetArray(vfw,&fw);
133: DMDAVecGetArray(user->da1,vu,&u);
134: DMDAVecGetArray(user->da1,vfu,&fu);
135: DMDAVecGetArray(user->da1,vlambda,&lambda);
136: DMDAVecGetArray(user->da1,vflambda,&flambda);
137: d = (N-1.0);
138: h = 1.0/d;
140: /* derivative of L() w.r.t. w */
141: if (xs == 0) { /* only first processor computes this */
142: fw[0] = -2.*d*lambda[0];
143: }
145: /* derivative of L() w.r.t. u */
146: for (i=xs; i<xs+xm; i++) {
147: if (i == 0) flambda[0] = h*u[0] + 2.*d*lambda[0] - d*lambda[1];
148: else if (i == 1) flambda[1] = 2.*h*u[1] + 2.*d*lambda[1] - d*lambda[2];
149: else if (i == N-1) flambda[N-1] = h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2];
150: else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3];
151: else flambda[i] = 2.*h*u[i] - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]);
152: }
154: /* derivative of L() w.r.t. lambda */
155: for (i=xs; i<xs+xm; i++) {
156: if (i == 0) fu[0] = 2.0*d*(u[0] - w[0]);
157: else if (i == N-1) fu[N-1] = 2.0*d*u[N-1];
158: else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);
159: }
161: VecRestoreArray(vw,&w);
162: VecRestoreArray(vfw,&fw);
163: DMDAVecRestoreArray(user->da1,vu,&u);
164: DMDAVecRestoreArray(user->da1,vfu,&fu);
165: DMDAVecRestoreArray(user->da1,vlambda,&lambda);
166: DMDAVecRestoreArray(user->da1,vflambda,&flambda);
168: DMCompositeGather(user->packer,INSERT_VALUES,FU,vfw,vfu,vflambda);
169: DMCompositeRestoreLocalVectors(user->packer,&vw,&vu,&vlambda);
170: DMCompositeRestoreLocalVectors(user->packer,&vfw,&vfu,&vflambda);
171: return(0);
172: }
174: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
175: {
176: UserCtx *user = (UserCtx*)dummy;
178: Vec w,u,lambda,U,F;
181: SNESGetSolution(snes,&U);
182: DMCompositeGetAccess(user->packer,U,&w,&u,&lambda);
183: VecView(u,user->u_viewer);
184: VecView(lambda,user->lambda_viewer);
185: DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda);
187: SNESGetFunction(snes,&F,0,0);
188: DMCompositeGetAccess(user->packer,F,&w,&u,&lambda);
189: VecView(u,user->fu_viewer);
190: VecView(lambda,user->flambda_viewer);
191: DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda);
192: return(0);
193: }
195: /*TEST
197: test:
198: nsize: 4
199: args: -snes_linesearch_monitor -snes_mf -snes_monitor_short -nox -ksp_monitor_short -snes_converged_reason
200: requires: !single
202: TEST*/