Actual source code: ex21.c

  2: static char help[] = "Solves PDE optimization problem.\n\n";

  4: #include <petscdmda.h>
  5: #include <petscdmcomposite.h>
  6: #include <petscpf.h>
  7: #include <petscsnes.h>

  9: /*

 11:        w - design variables (what we change to get an optimal solution)
 12:        u - state variables (i.e. the PDE solution)
 13:        lambda - the Lagrange multipliers

 15:             U = (w u lambda)

 17:        fu, fw, flambda contain the gradient of L(w,u,lambda)

 19:             FU = (fw fu flambda)

 21:        In this example the PDE is 
 22:                              Uxx = 2, 
 23:                             u(0) = w(0), thus this is the free parameter
 24:                             u(1) = 0
 25:        the function we wish to minimize is 
 26:                             \integral u^{2}

 28:        The exact solution for u is given by u(x) = x*x - 1.25*x + .25

 30:        Use the usual centered finite differences.

 32:        Note we treat the problem as non-linear though it happens to be linear

 34:        See ex22.c for the same code, but that interlaces the u and the lambda

 36: */

 38: typedef struct {
 39:   DM           da1,da2;
 40:   PetscInt     nredundant;
 41:   DM           packer;
 42:   PetscViewer  u_viewer,lambda_viewer;
 43:   PetscViewer  fu_viewer,flambda_viewer;
 44: } UserCtx;



 52: int main(int argc,char **argv)
 53: {
 55:   PetscInt       its;
 56:   Vec            U,FU;
 57:   SNES           snes;
 58:   UserCtx        user;

 60:   PetscInitialize(&argc,&argv,(char*)0,help);

 62:   /* Create a global vector that includes a single redundant array and two da arrays */
 63:   DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);
 64:   user.nredundant = 1;
 65:   DMCompositeAddArray(user.packer,0,user.nredundant);
 66:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-5,1,1,PETSC_NULL,&user.da1);
 67:   DMCompositeAddDM(user.packer,(DM)user.da1);
 68:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-5,1,1,PETSC_NULL,&user.da2);
 69:   DMCompositeAddDM(user.packer,(DM)user.da2);
 70:   DMCreateGlobalVector(user.packer,&U);
 71:   VecDuplicate(U,&FU);

 73:   /* create graphics windows */
 74:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
 75:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);
 76:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);
 77:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);


 80:   /* create nonlinear solver */
 81:   SNESCreate(PETSC_COMM_WORLD,&snes);
 82:   SNESSetFunction(snes,FU,FormFunction,&user);
 83:   SNESSetFromOptions(snes);
 84:   SNESMonitorSet(snes,Monitor,&user,0);
 85:   SNESSolve(snes,PETSC_NULL,U);
 86:   SNESGetIterationNumber(snes,&its);
 87:   SNESDestroy(&snes);

 89:   DMDestroy(&user.da1);
 90:   DMDestroy(&user.da2);
 91:   DMDestroy(&user.packer);
 92:   VecDestroy(&U);
 93:   VecDestroy(&FU);
 94:   PetscViewerDestroy(&user.u_viewer);
 95:   PetscViewerDestroy(&user.lambda_viewer);
 96:   PetscViewerDestroy(&user.fu_viewer);
 97:   PetscViewerDestroy(&user.flambda_viewer);
 98:   PetscFinalize();
 99:   return 0;
100: }
101: 
104: /*
105:       Evaluates FU = Gradiant(L(w,u,lambda))

107: */
108: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
109: {
110:   UserCtx        *user = (UserCtx*)dummy;
112:   PetscInt       xs,xm,i,N;
113:   PetscScalar    *u,*lambda,*w,*fu,*fw,*flambda,d,h;
114:   Vec            vu,vlambda,vfu,vflambda;

117:   DMCompositeGetLocalVectors(user->packer,&w,&vu,&vlambda);
118:   DMCompositeGetLocalVectors(user->packer,&fw,&vfu,&vflambda);
119:   DMCompositeScatter(user->packer,U,w,vu,vlambda);

121:   DMDAGetCorners(user->da1,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
122:   DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0);
123:   DMDAVecGetArray(user->da1,vu,&u);
124:   DMDAVecGetArray(user->da1,vfu,&fu);
125:   DMDAVecGetArray(user->da1,vlambda,&lambda);
126:   DMDAVecGetArray(user->da1,vflambda,&flambda);
127:   d    = (N-1.0);
128:   h    = 1.0/d;

130:   /* derivative of L() w.r.t. w */
131:   if (xs == 0) { /* only first processor computes this */
132:     fw[0] = -2.*d*lambda[0];
133:   }

135:   /* derivative of L() w.r.t. u */
136:   for (i=xs; i<xs+xm; i++) {
137:     if      (i == 0)   flambda[0]   =    h*u[0]   + 2.*d*lambda[0]   - d*lambda[1];
138:     else if (i == 1)   flambda[1]   = 2.*h*u[1]   + 2.*d*lambda[1]   - d*lambda[2];
139:     else if (i == N-1) flambda[N-1] =    h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2];
140:     else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3];
141:     else               flambda[i]   = 2.*h*u[i]   - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]);
142:   }

144:   /* derivative of L() w.r.t. lambda */
145:   for (i=xs; i<xs+xm; i++) {
146:     if      (i == 0)   fu[0]   = 2.0*d*(u[0] - w[0]);
147:     else if (i == N-1) fu[N-1] = 2.0*d*u[N-1];
148:     else               fu[i]   = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);
149:   }

151:   DMDAVecRestoreArray(user->da1,vu,&u);
152:   DMDAVecRestoreArray(user->da1,vfu,&fu);
153:   DMDAVecRestoreArray(user->da1,vlambda,&lambda);
154:   DMDAVecRestoreArray(user->da1,vflambda,&flambda);

156:   DMCompositeGather(user->packer,FU,INSERT_VALUES,fw,vfu,vflambda);
157:   DMCompositeRestoreLocalVectors(user->packer,&w,&vu,&vlambda);
158:   DMCompositeRestoreLocalVectors(user->packer,&fw,&vfu,&vflambda);
159:   return(0);
160: }

164: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
165: {
166:   UserCtx        *user = (UserCtx*)dummy;
168:   PetscScalar    *w;
169:   Vec            u,lambda,U,F;

172:   SNESGetSolution(snes,&U);
173:   DMCompositeGetAccess(user->packer,U,&w,&u,&lambda);
174:   VecView(u,user->u_viewer);
175:   VecView(lambda,user->lambda_viewer);
176:   DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda);

178:   SNESGetFunction(snes,&F,0,0);
179:   DMCompositeGetAccess(user->packer,F,&w,&u,&lambda);
180:   VecView(u,user->fu_viewer);
181:   VecView(lambda,user->flambda_viewer);
182:   DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda);
183:   return(0);
184: }