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 { 76 UserCtx user; 77 DM red, da; 78 SNES snes; 79 DM packer; 80 PetscBool use_monitor = PETSC_FALSE; 81 82 PetscFunctionBeginUser; 83 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 84 85 /* Hardwire several options; can be changed at command line */ 86 PetscCall(PetscOptionsInsertString(NULL, common_options)); 87 PetscCall(PetscOptionsInsertString(NULL, matrix_free_options)); 88 PetscCall(PetscOptionsInsert(NULL, &argc, &argv, NULL)); 89 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_monitor", &use_monitor, PETSC_IGNORE)); 90 91 /* Create a global vector that includes a single redundant array and two da arrays */ 92 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &packer)); 93 PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 0, 1, &red)); 94 PetscCall(DMSetOptionsPrefix(red, "red_")); 95 PetscCall(DMCompositeAddDM(packer, red)); 96 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 5, 2, 1, NULL, &da)); 97 PetscCall(DMSetOptionsPrefix(red, "da_")); 98 PetscCall(DMSetFromOptions(da)); 99 PetscCall(DMSetUp(da)); 100 PetscCall(DMDASetFieldName(da, 0, "u")); 101 PetscCall(DMDASetFieldName(da, 1, "lambda")); 102 PetscCall(DMCompositeAddDM(packer, (DM)da)); 103 PetscCall(DMSetApplicationContext(packer, &user)); 104 105 packer->ops->creatematrix = DMCreateMatrix_MF; 106 107 /* create nonlinear multi-level solver */ 108 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 109 PetscCall(SNESSetDM(snes, packer)); 110 PetscCall(SNESSetFunction(snes, NULL, ComputeFunction, NULL)); 111 PetscCall(SNESSetJacobian(snes, NULL, NULL, ComputeJacobian_MF, NULL)); 112 113 PetscCall(SNESSetFromOptions(snes)); 114 115 if (use_monitor) { 116 /* create graphics windows */ 117 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "u_lambda - state variables and Lagrange multipliers", -1, -1, -1, -1, &user.u_lambda_viewer)); 118 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)); 119 PetscCall(SNESMonitorSet(snes, Monitor, 0, 0)); 120 } 121 122 PetscCall(SNESSolve(snes, NULL, NULL)); 123 PetscCall(SNESDestroy(&snes)); 124 125 PetscCall(DMDestroy(&red)); 126 PetscCall(DMDestroy(&da)); 127 PetscCall(DMDestroy(&packer)); 128 if (use_monitor) { 129 PetscCall(PetscViewerDestroy(&user.u_lambda_viewer)); 130 PetscCall(PetscViewerDestroy(&user.fu_lambda_viewer)); 131 } 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 typedef struct { 137 PetscScalar u; 138 PetscScalar lambda; 139 } ULambda; 140 141 /* 142 Evaluates FU = Gradient(L(w,u,lambda)) 143 144 This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and 145 DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()). 146 147 */ 148 PetscErrorCode ComputeFunction(SNES snes, Vec U, Vec FU, void *ctx) 149 { 150 PetscInt xs, xm, i, N; 151 ULambda *u_lambda, *fu_lambda; 152 PetscScalar d, h, *w, *fw; 153 Vec vw, vfw, vu_lambda, vfu_lambda; 154 DM packer, red, da; 155 156 PetscFunctionBeginUser; 157 PetscCall(VecGetDM(U, &packer)); 158 PetscCall(DMCompositeGetEntries(packer, &red, &da)); 159 PetscCall(DMCompositeGetLocalVectors(packer, &vw, &vu_lambda)); 160 PetscCall(DMCompositeScatter(packer, U, vw, vu_lambda)); 161 PetscCall(DMCompositeGetAccess(packer, FU, &vfw, &vfu_lambda)); 162 163 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 164 PetscCall(DMDAGetInfo(da, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 165 PetscCall(VecGetArray(vw, &w)); 166 PetscCall(VecGetArray(vfw, &fw)); 167 PetscCall(DMDAVecGetArray(da, vu_lambda, &u_lambda)); 168 PetscCall(DMDAVecGetArray(da, vfu_lambda, &fu_lambda)); 169 d = N - 1.0; 170 h = 1.0 / d; 171 172 /* derivative of L() w.r.t. w */ 173 if (xs == 0) { /* only first processor computes this */ 174 fw[0] = -2.0 * d * u_lambda[0].lambda; 175 } 176 177 /* derivative of L() w.r.t. u */ 178 for (i = xs; i < xs + xm; i++) { 179 if (i == 0) fu_lambda[0].lambda = h * u_lambda[0].u + 2. * d * u_lambda[0].lambda - d * u_lambda[1].lambda; 180 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; 181 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; 182 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; 183 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); 184 } 185 186 /* derivative of L() w.r.t. lambda */ 187 for (i = xs; i < xs + xm; i++) { 188 if (i == 0) fu_lambda[0].u = 2.0 * d * (u_lambda[0].u - w[0]); 189 else if (i == N - 1) fu_lambda[N - 1].u = 2.0 * d * u_lambda[N - 1].u; 190 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); 191 } 192 193 PetscCall(VecRestoreArray(vw, &w)); 194 PetscCall(VecRestoreArray(vfw, &fw)); 195 PetscCall(DMDAVecRestoreArray(da, vu_lambda, &u_lambda)); 196 PetscCall(DMDAVecRestoreArray(da, vfu_lambda, &fu_lambda)); 197 PetscCall(DMCompositeRestoreLocalVectors(packer, &vw, &vu_lambda)); 198 PetscCall(DMCompositeRestoreAccess(packer, FU, &vfw, &vfu_lambda)); 199 PetscCall(PetscLogFlops(13.0 * N)); 200 PetscFunctionReturn(0); 201 } 202 203 /* 204 Computes the exact solution 205 */ 206 PetscErrorCode u_solution(void *dummy, PetscInt n, const PetscScalar *x, PetscScalar *u) 207 { 208 PetscInt i; 209 210 PetscFunctionBeginUser; 211 for (i = 0; i < n; i++) u[2 * i] = x[i] * x[i] - 1.25 * x[i] + .25; 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode ExactSolution(DM packer, Vec U) 216 { 217 PF pf; 218 Vec x, u_global; 219 PetscScalar *w; 220 DM da; 221 PetscInt m; 222 223 PetscFunctionBeginUser; 224 PetscCall(DMCompositeGetEntries(packer, &m, &da)); 225 226 PetscCall(PFCreate(PETSC_COMM_WORLD, 1, 2, &pf)); 227 /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */ 228 PetscCall(PFSetType(pf, PFQUICK, (void *)(PETSC_UINTPTR_T)u_solution)); 229 PetscCall(DMGetCoordinates(da, &x)); 230 if (!x) { 231 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 232 PetscCall(DMGetCoordinates(da, &x)); 233 } 234 PetscCall(DMCompositeGetAccess(packer, U, &w, &u_global, 0)); 235 if (w) w[0] = .25; 236 PetscCall(PFApplyVec(pf, x, u_global)); 237 PetscCall(PFDestroy(&pf)); 238 PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_global, 0)); 239 PetscFunctionReturn(0); 240 } 241 242 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal rnorm, void *dummy) 243 { 244 UserCtx *user; 245 PetscInt m, N; 246 PetscScalar *w, *dw; 247 Vec u_lambda, U, F, Uexact; 248 DM packer; 249 PetscReal norm; 250 DM da; 251 252 PetscFunctionBeginUser; 253 PetscCall(SNESGetDM(snes, &packer)); 254 PetscCall(DMGetApplicationContext(packer, &user)); 255 PetscCall(SNESGetSolution(snes, &U)); 256 PetscCall(DMCompositeGetAccess(packer, U, &w, &u_lambda)); 257 PetscCall(VecView(u_lambda, user->u_lambda_viewer)); 258 PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_lambda)); 259 260 PetscCall(SNESGetFunction(snes, &F, 0, 0)); 261 PetscCall(DMCompositeGetAccess(packer, F, &w, &u_lambda)); 262 /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */ 263 PetscCall(DMCompositeRestoreAccess(packer, U, &w, &u_lambda)); 264 265 PetscCall(DMCompositeGetEntries(packer, &m, &da)); 266 PetscCall(DMDAGetInfo(da, 0, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 267 PetscCall(VecDuplicate(U, &Uexact)); 268 PetscCall(ExactSolution(packer, Uexact)); 269 PetscCall(VecAXPY(Uexact, -1.0, U)); 270 PetscCall(DMCompositeGetAccess(packer, Uexact, &dw, &u_lambda)); 271 PetscCall(VecStrideNorm(u_lambda, 0, NORM_2, &norm)); 272 norm = norm / PetscSqrtReal((PetscReal)N - 1.); 273 if (dw) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Error at x = 0 %g\n", (double)norm, (double)PetscRealPart(dw[0]))); 274 PetscCall(VecView(u_lambda, user->fu_lambda_viewer)); 275 PetscCall(DMCompositeRestoreAccess(packer, Uexact, &dw, &u_lambda)); 276 PetscCall(VecDestroy(&Uexact)); 277 PetscFunctionReturn(0); 278 } 279 280 PetscErrorCode DMCreateMatrix_MF(DM packer, Mat *A) 281 { 282 Vec t; 283 PetscInt m; 284 285 PetscFunctionBeginUser; 286 PetscCall(DMGetGlobalVector(packer, &t)); 287 PetscCall(VecGetLocalSize(t, &m)); 288 PetscCall(DMRestoreGlobalVector(packer, &t)); 289 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, A)); 290 PetscCall(MatSetUp(*A)); 291 PetscCall(MatSetDM(*A, packer)); 292 PetscFunctionReturn(0); 293 } 294 295 PetscErrorCode ComputeJacobian_MF(SNES snes, Vec x, Mat A, Mat B, void *ctx) 296 { 297 PetscFunctionBeginUser; 298 PetscCall(MatMFFDSetFunction(A, (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunction, snes)); 299 PetscCall(MatMFFDSetBase(A, x, NULL)); 300 PetscFunctionReturn(0); 301 } 302 303 /*TEST 304 305 test: 306 nsize: 2 307 args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view 308 requires: !single 309 310 TEST*/ 311