Actual source code: ex21.c


  2: static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n";

  4: #include <petscdm.h>
  5: #include <petscdmda.h>
  6: #include <petscdmredundant.h>
  7: #include <petscdmcomposite.h>
  8: #include <petscpf.h>
  9: #include <petscsnes.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 lambda)

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

 21:             FU = (fw fu flambda)

 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 ex22.c for the same code, but that interlaces the u and the lambda

 38: */

 40: typedef struct {
 41:   DM          red1, da1, da2;
 42:   DM          packer;
 43:   PetscViewer u_viewer, lambda_viewer;
 44:   PetscViewer fu_viewer, flambda_viewer;
 45: } UserCtx;

 47: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
 48: extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);

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

 57:   PetscFunctionBeginUser;
 58:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 60:   /* Create a global vector that includes a single redundant array and two da arrays */
 61:   PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.packer));
 62:   PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 0, 1, &user.red1));
 63:   PetscCall(DMCompositeAddDM(user.packer, user.red1));
 64:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da1));
 65:   PetscCall(DMSetFromOptions(user.da1));
 66:   PetscCall(DMSetUp(user.da1));
 67:   PetscCall(DMCompositeAddDM(user.packer, user.da1));
 68:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 1, 1, NULL, &user.da2));
 69:   PetscCall(DMSetFromOptions(user.da2));
 70:   PetscCall(DMSetUp(user.da2));
 71:   PetscCall(DMDASetFieldName(user.da1, 0, "u"));
 72:   PetscCall(DMDASetFieldName(user.da2, 0, "lambda"));
 73:   PetscCall(DMCompositeAddDM(user.packer, user.da2));
 74:   PetscCall(DMCreateGlobalVector(user.packer, &U));
 75:   PetscCall(VecDuplicate(U, &FU));

 77:   /* create graphics windows */
 78:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "u - state variables", -1, -1, -1, -1, &user.u_viewer));
 79:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "lambda - Lagrange multipliers", -1, -1, -1, -1, &user.lambda_viewer));
 80:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "fu - derivative w.r.t. state variables", -1, -1, -1, -1, &user.fu_viewer));
 81:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "flambda - derivative w.r.t. Lagrange multipliers", -1, -1, -1, -1, &user.flambda_viewer));

 83:   /* create nonlinear solver */
 84:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
 85:   PetscCall(SNESSetFunction(snes, FU, FormFunction, &user));
 86:   PetscCall(SNESSetFromOptions(snes));
 87:   PetscCall(SNESMonitorSet(snes, Monitor, &user, 0));
 88:   PetscCall(SNESSolve(snes, NULL, U));
 89:   PetscCall(SNESGetIterationNumber(snes, &its));
 90:   PetscCall(SNESDestroy(&snes));

 92:   PetscCall(DMDestroy(&user.red1));
 93:   PetscCall(DMDestroy(&user.da1));
 94:   PetscCall(DMDestroy(&user.da2));
 95:   PetscCall(DMDestroy(&user.packer));
 96:   PetscCall(VecDestroy(&U));
 97:   PetscCall(VecDestroy(&FU));
 98:   PetscCall(PetscViewerDestroy(&user.u_viewer));
 99:   PetscCall(PetscViewerDestroy(&user.lambda_viewer));
100:   PetscCall(PetscViewerDestroy(&user.fu_viewer));
101:   PetscCall(PetscViewerDestroy(&user.flambda_viewer));
102:   PetscCall(PetscFinalize());
103:   return 0;
104: }

106: /*
107:       Evaluates FU = Gradient(L(w,u,lambda))

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

117:   PetscFunctionBeginUser;
118:   PetscCall(DMCompositeGetLocalVectors(user->packer, &vw, &vu, &vlambda));
119:   PetscCall(DMCompositeGetLocalVectors(user->packer, &vfw, &vfu, &vflambda));
120:   PetscCall(DMCompositeScatter(user->packer, U, vw, vu, vlambda));

122:   PetscCall(DMDAGetCorners(user->da1, &xs, NULL, NULL, &xm, NULL, NULL));
123:   PetscCall(DMDAGetInfo(user->da1, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
124:   PetscCall(VecGetArray(vw, &w));
125:   PetscCall(VecGetArray(vfw, &fw));
126:   PetscCall(DMDAVecGetArray(user->da1, vu, &u));
127:   PetscCall(DMDAVecGetArray(user->da1, vfu, &fu));
128:   PetscCall(DMDAVecGetArray(user->da1, vlambda, &lambda));
129:   PetscCall(DMDAVecGetArray(user->da1, vflambda, &flambda));
130:   d = (N - 1.0);
131:   h = 1.0 / d;

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

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

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

154:   PetscCall(VecRestoreArray(vw, &w));
155:   PetscCall(VecRestoreArray(vfw, &fw));
156:   PetscCall(DMDAVecRestoreArray(user->da1, vu, &u));
157:   PetscCall(DMDAVecRestoreArray(user->da1, vfu, &fu));
158:   PetscCall(DMDAVecRestoreArray(user->da1, vlambda, &lambda));
159:   PetscCall(DMDAVecRestoreArray(user->da1, vflambda, &flambda));

161:   PetscCall(DMCompositeGather(user->packer, INSERT_VALUES, FU, vfw, vfu, vflambda));
162:   PetscCall(DMCompositeRestoreLocalVectors(user->packer, &vw, &vu, &vlambda));
163:   PetscCall(DMCompositeRestoreLocalVectors(user->packer, &vfw, &vfu, &vflambda));
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

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

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

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

187: /*TEST

189:    test:
190:       nsize: 4
191:       args: -snes_linesearch_monitor -snes_mf -pc_type none -snes_monitor_short -nox -ksp_monitor_short -snes_converged_reason
192:       requires: !single

194: TEST*/