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