1 2 static const char help[] = "Solves PDE optimization problem using full-space method, interlaces state and adjoint variables.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 #include <petscdmredundant.h> 7 #include <petscdmcomposite.h> 8 #include <petscpf.h> 9 #include <petscsnes.h> 10 #include <petsc/private/dmimpl.h> 11 12 /* 13 14 w - design variables (what we change to get an optimal solution) 15 u - state variables (i.e. the PDE solution) 16 lambda - the Lagrange multipliers 17 18 U = (w [u_0 lambda_0 u_1 lambda_1 .....]) 19 20 fu, fw, flambda contain the gradient of L(w,u,lambda) 21 22 FU = (fw [fu_0 flambda_0 .....]) 23 24 In this example the PDE is 25 Uxx = 2, 26 u(0) = w(0), thus this is the free parameter 27 u(1) = 0 28 the function we wish to minimize is 29 \integral u^{2} 30 31 The exact solution for u is given by u(x) = x*x - 1.25*x + .25 32 33 Use the usual centered finite differences. 34 35 Note we treat the problem as non-linear though it happens to be linear 36 37 See ex21.c for the same code, but that does NOT interlaces the u and the lambda 38 39 The vectors u_lambda and fu_lambda contain the u and the lambda interlaced 40 */ 41 42 typedef struct { 43 PetscViewer u_lambda_viewer; 44 PetscViewer fu_lambda_viewer; 45 } UserCtx; 46 47 extern PetscErrorCode ComputeFunction(SNES, Vec, Vec, void *); 48 extern PetscErrorCode ComputeJacobian_MF(SNES, Vec, Mat, Mat, void *); 49 extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *); 50 51 /* 52 Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the 53 smoother on all levels. This is because (1) in the matrix free case no matrix entries are 54 available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal 55 entry for the control variable is zero which means default SOR will not work. 56 57 */ 58 char common_options[] = "-ksp_type fgmres\ 59 -snes_grid_sequence 2 \ 60 -pc_type mg\ 61 -mg_levels_pc_type none \ 62 -mg_coarse_pc_type none \ 63 -pc_mg_type full \ 64 -mg_coarse_ksp_type gmres \ 65 -mg_levels_ksp_type gmres \ 66 -mg_coarse_ksp_max_it 6 \ 67 -mg_levels_ksp_max_it 3"; 68 69 char matrix_free_options[] = "-mat_mffd_compute_normu no \ 70 -mat_mffd_type wp"; 71 72 extern PetscErrorCode DMCreateMatrix_MF(DM, Mat *); 73 74 int main(int argc, char **argv) { 75 UserCtx user; 76 DM red, da; 77 SNES snes; 78 DM packer; 79 PetscBool use_monitor = PETSC_FALSE; 80 81 PetscFunctionBeginUser; 82 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 83 84 /* Hardwire several options; can be changed at command line */ 85 PetscCall(PetscOptionsInsertString(NULL, common_options)); 86 PetscCall(PetscOptionsInsertString(NULL, matrix_free_options)); 87 PetscCall(PetscOptionsInsert(NULL, &argc, &argv, NULL)); 88 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_monitor", &use_monitor, PETSC_IGNORE)); 89 90 /* Create a global vector that includes a single redundant array and two da arrays */ 91 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &packer)); 92 PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 0, 1, &red)); 93 PetscCall(DMSetOptionsPrefix(red, "red_")); 94 PetscCall(DMCompositeAddDM(packer, red)); 95 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 2, 1, NULL, &da)); 96 PetscCall(DMSetOptionsPrefix(red, "da_")); 97 PetscCall(DMSetFromOptions(da)); 98 PetscCall(DMSetUp(da)); 99 PetscCall(DMDASetFieldName(da, 0, "u")); 100 PetscCall(DMDASetFieldName(da, 1, "lambda")); 101 PetscCall(DMCompositeAddDM(packer, (DM)da)); 102 PetscCall(DMSetApplicationContext(packer, &user)); 103 104 packer->ops->creatematrix = DMCreateMatrix_MF; 105 106 /* create nonlinear multi-level solver */ 107 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 108 PetscCall(SNESSetDM(snes, packer)); 109 PetscCall(SNESSetFunction(snes, NULL, ComputeFunction, NULL)); 110 PetscCall(SNESSetJacobian(snes, NULL, NULL, ComputeJacobian_MF, NULL)); 111 112 PetscCall(SNESSetFromOptions(snes)); 113 114 if (use_monitor) { 115 /* create graphics windows */ 116 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "u_lambda - state variables and Lagrange multipliers", -1, -1, -1, -1, &user.u_lambda_viewer)); 117 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "fu_lambda - derivate w.r.t. state variables and Lagrange multipliers", -1, -1, -1, -1, &user.fu_lambda_viewer)); 118 PetscCall(SNESMonitorSet(snes, Monitor, 0, 0)); 119 } 120 121 PetscCall(SNESSolve(snes, NULL, NULL)); 122 PetscCall(SNESDestroy(&snes)); 123 124 PetscCall(DMDestroy(&red)); 125 PetscCall(DMDestroy(&da)); 126 PetscCall(DMDestroy(&packer)); 127 if (use_monitor) { 128 PetscCall(PetscViewerDestroy(&user.u_lambda_viewer)); 129 PetscCall(PetscViewerDestroy(&user.fu_lambda_viewer)); 130 } 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 typedef struct { 136 PetscScalar u; 137 PetscScalar lambda; 138 } ULambda; 139 140 /* 141 Evaluates FU = Gradiant(L(w,u,lambda)) 142 143 This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and 144 DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()). 145 146 */ 147 PetscErrorCode ComputeFunction(SNES snes, Vec U, Vec FU, void *ctx) { 148 PetscInt xs, xm, i, N; 149 ULambda *u_lambda, *fu_lambda; 150 PetscScalar d, h, *w, *fw; 151 Vec vw, vfw, vu_lambda, vfu_lambda; 152 DM packer, red, da; 153 154 PetscFunctionBeginUser; 155 PetscCall(VecGetDM(U, &packer)); 156 PetscCall(DMCompositeGetEntries(packer, &red, &da)); 157 PetscCall(DMCompositeGetLocalVectors(packer, &vw, &vu_lambda)); 158 PetscCall(DMCompositeScatter(packer, U, vw, vu_lambda)); 159 PetscCall(DMCompositeGetAccess(packer, FU, &vfw, &vfu_lambda)); 160 161 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 162 PetscCall(DMDAGetInfo(da, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 163 PetscCall(VecGetArray(vw, &w)); 164 PetscCall(VecGetArray(vfw, &fw)); 165 PetscCall(DMDAVecGetArray(da, vu_lambda, &u_lambda)); 166 PetscCall(DMDAVecGetArray(da, vfu_lambda, &fu_lambda)); 167 d = N - 1.0; 168 h = 1.0 / d; 169 170 /* derivative of L() w.r.t. w */ 171 if (xs == 0) { /* only first processor computes this */ 172 fw[0] = -2.0 * d * u_lambda[0].lambda; 173 } 174 175 /* derivative of L() w.r.t. u */ 176 for (i = xs; i < xs + xm; i++) { 177 if (i == 0) fu_lambda[0].lambda = h * u_lambda[0].u + 2. * d * u_lambda[0].lambda - d * u_lambda[1].lambda; 178 else if (i == 1) fu_lambda[1].lambda = 2. * h * u_lambda[1].u + 2. * d * u_lambda[1].lambda - d * u_lambda[2].lambda; 179 else if (i == N - 1) fu_lambda[N - 1].lambda = h * u_lambda[N - 1].u + 2. * d * u_lambda[N - 1].lambda - d * u_lambda[N - 2].lambda; 180 else if (i == N - 2) fu_lambda[N - 2].lambda = 2. * h * u_lambda[N - 2].u + 2. * d * u_lambda[N - 2].lambda - d * u_lambda[N - 3].lambda; 181 else fu_lambda[i].lambda = 2. * h * u_lambda[i].u - d * (u_lambda[i + 1].lambda - 2.0 * u_lambda[i].lambda + u_lambda[i - 1].lambda); 182 } 183 184 /* derivative of L() w.r.t. lambda */ 185 for (i = xs; i < xs + xm; i++) { 186 if (i == 0) fu_lambda[0].u = 2.0 * d * (u_lambda[0].u - w[0]); 187 else if (i == N - 1) fu_lambda[N - 1].u = 2.0 * d * u_lambda[N - 1].u; 188 else fu_lambda[i].u = -(d * (u_lambda[i + 1].u - 2.0 * u_lambda[i].u + u_lambda[i - 1].u) - 2.0 * h); 189 } 190 191 PetscCall(VecRestoreArray(vw, &w)); 192 PetscCall(VecRestoreArray(vfw, &fw)); 193 PetscCall(DMDAVecRestoreArray(da, vu_lambda, &u_lambda)); 194 PetscCall(DMDAVecRestoreArray(da, vfu_lambda, &fu_lambda)); 195 PetscCall(DMCompositeRestoreLocalVectors(packer, &vw, &vu_lambda)); 196 PetscCall(DMCompositeRestoreAccess(packer, FU, &vfw, &vfu_lambda)); 197 PetscCall(PetscLogFlops(13.0 * N)); 198 PetscFunctionReturn(0); 199 } 200 201 /* 202 Computes the exact solution 203 */ 204 PetscErrorCode u_solution(void *dummy, PetscInt n, const PetscScalar *x, PetscScalar *u) { 205 PetscInt i; 206 207 PetscFunctionBeginUser; 208 for (i = 0; i < n; i++) u[2 * i] = x[i] * x[i] - 1.25 * x[i] + .25; 209 PetscFunctionReturn(0); 210 } 211 212 PetscErrorCode ExactSolution(DM packer, Vec U) { 213 PF pf; 214 Vec x, u_global; 215 PetscScalar *w; 216 DM da; 217 PetscInt m; 218 219 PetscFunctionBeginUser; 220 PetscCall(DMCompositeGetEntries(packer, &m, &da)); 221 222 PetscCall(PFCreate(PETSC_COMM_WORLD, 1, 2, &pf)); 223 /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */ 224 PetscCall(PFSetType(pf, PFQUICK, (void *)(PETSC_UINTPTR_T)u_solution)); 225 PetscCall(DMGetCoordinates(da, &x)); 226 if (!x) { 227 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 228 PetscCall(DMGetCoordinates(da, &x)); 229 } 230 PetscCall(DMCompositeGetAccess(packer, U, &w, &u_global, 0)); 231 if (w) w[0] = .25; 232 PetscCall(PFApplyVec(pf, x, u_global)); 233 PetscCall(PFDestroy(&pf)); 234 PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_global, 0)); 235 PetscFunctionReturn(0); 236 } 237 238 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal rnorm, void *dummy) { 239 UserCtx *user; 240 PetscInt m, N; 241 PetscScalar *w, *dw; 242 Vec u_lambda, U, F, Uexact; 243 DM packer; 244 PetscReal norm; 245 DM da; 246 247 PetscFunctionBeginUser; 248 PetscCall(SNESGetDM(snes, &packer)); 249 PetscCall(DMGetApplicationContext(packer, &user)); 250 PetscCall(SNESGetSolution(snes, &U)); 251 PetscCall(DMCompositeGetAccess(packer, U, &w, &u_lambda)); 252 PetscCall(VecView(u_lambda, user->u_lambda_viewer)); 253 PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_lambda)); 254 255 PetscCall(SNESGetFunction(snes, &F, 0, 0)); 256 PetscCall(DMCompositeGetAccess(packer, F, &w, &u_lambda)); 257 /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */ 258 PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_lambda)); 259 260 PetscCall(DMCompositeGetEntries(packer, &m, &da)); 261 PetscCall(DMDAGetInfo(da, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 262 PetscCall(VecDuplicate(U, &Uexact)); 263 PetscCall(ExactSolution(packer, Uexact)); 264 PetscCall(VecAXPY(Uexact, -1.0, U)); 265 PetscCall(DMCompositeGetAccess(packer, Uexact, &dw, &u_lambda)); 266 PetscCall(VecStrideNorm(u_lambda, 0, NORM_2, &norm)); 267 norm = norm / PetscSqrtReal((PetscReal)N - 1.); 268 if (dw) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Error at x = 0 %g\n", (double)norm, (double)PetscRealPart(dw[0]))); 269 PetscCall(VecView(u_lambda, user->fu_lambda_viewer)); 270 PetscCall(DMCompositeRestoreAccess(packer, Uexact, &dw, &u_lambda)); 271 PetscCall(VecDestroy(&Uexact)); 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode DMCreateMatrix_MF(DM packer, Mat *A) { 276 Vec t; 277 PetscInt m; 278 279 PetscFunctionBeginUser; 280 PetscCall(DMGetGlobalVector(packer, &t)); 281 PetscCall(VecGetLocalSize(t, &m)); 282 PetscCall(DMRestoreGlobalVector(packer, &t)); 283 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, A)); 284 PetscCall(MatSetUp(*A)); 285 PetscCall(MatSetDM(*A, packer)); 286 PetscFunctionReturn(0); 287 } 288 289 PetscErrorCode ComputeJacobian_MF(SNES snes, Vec x, Mat A, Mat B, void *ctx) { 290 PetscFunctionBeginUser; 291 PetscCall(MatMFFDSetFunction(A, (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunction, snes)); 292 PetscCall(MatMFFDSetBase(A, x, NULL)); 293 PetscFunctionReturn(0); 294 } 295 296 /*TEST 297 298 test: 299 nsize: 2 300 args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view 301 requires: !single 302 303 TEST*/ 304