Actual source code: ex21.c

petsc-3.10.2 2018-10-09
Report Typos and Errors

  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*/