1 #include <petsctao.h> 2 3 typedef struct { 4 PetscInt n; /* Number of variables */ 5 PetscInt m; /* Number of constraints */ 6 PetscInt mx; /* grid points in each direction */ 7 PetscInt nt; /* Number of time steps */ 8 PetscInt ndata; /* Number of data points per sample */ 9 IS s_is; 10 IS d_is; 11 VecScatter state_scatter; 12 VecScatter design_scatter; 13 VecScatter *uxi_scatter, *uyi_scatter, *ux_scatter, *uy_scatter, *ui_scatter; 14 VecScatter *yi_scatter; 15 16 Mat Js, Jd, JsBlockPrec, JsInv, JsBlock; 17 PetscBool jformed, c_formed; 18 19 PetscReal alpha; /* Regularization parameter */ 20 PetscReal gamma; 21 PetscReal ht; /* Time step */ 22 PetscReal T; /* Final time */ 23 Mat Q, QT; 24 Mat L, LT; 25 Mat Div, Divwork, Divxy[2]; 26 Mat Grad, Gradxy[2]; 27 Mat M; 28 Mat *C, *Cwork; 29 /* Mat Hs,Hd,Hsd; */ 30 Vec q; 31 Vec ur; /* reference */ 32 33 Vec d; 34 Vec dwork; 35 36 Vec y; /* state variables */ 37 Vec ywork; 38 Vec ytrue; 39 Vec *yi, *yiwork, *ziwork; 40 Vec *uxi, *uyi, *uxiwork, *uyiwork, *ui, *uiwork; 41 42 Vec u; /* design variables */ 43 Vec uwork, vwork; 44 Vec utrue; 45 46 Vec js_diag; 47 48 Vec c; /* constraint vector */ 49 Vec cwork; 50 51 Vec lwork; 52 53 KSP solver; 54 PC prec; 55 PetscInt block_index; 56 57 PetscInt ksp_its; 58 PetscInt ksp_its_initial; 59 } AppCtx; 60 61 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *); 62 PetscErrorCode FormGradient(Tao, Vec, Vec, void *); 63 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *); 64 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void *); 65 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void *); 66 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *); 67 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *); 68 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 69 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 70 PetscErrorCode HyperbolicInitialize(AppCtx *user); 71 PetscErrorCode HyperbolicDestroy(AppCtx *user); 72 PetscErrorCode HyperbolicMonitor(Tao, void *); 73 74 PetscErrorCode StateMatMult(Mat, Vec, Vec); 75 PetscErrorCode StateMatBlockMult(Mat, Vec, Vec); 76 PetscErrorCode StateMatBlockMultTranspose(Mat, Vec, Vec); 77 PetscErrorCode StateMatMultTranspose(Mat, Vec, Vec); 78 PetscErrorCode StateMatGetDiagonal(Mat, Vec); 79 PetscErrorCode StateMatDuplicate(Mat, MatDuplicateOption, Mat *); 80 PetscErrorCode StateMatInvMult(Mat, Vec, Vec); 81 PetscErrorCode StateMatInvTransposeMult(Mat, Vec, Vec); 82 PetscErrorCode StateMatBlockPrecMult(PC, Vec, Vec); 83 84 PetscErrorCode DesignMatMult(Mat, Vec, Vec); 85 PetscErrorCode DesignMatMultTranspose(Mat, Vec, Vec); 86 87 PetscErrorCode Scatter_yi(Vec, Vec *, VecScatter *, PetscInt); /* y to y1,y2,...,y_nt */ 88 PetscErrorCode Gather_yi(Vec, Vec *, VecScatter *, PetscInt); 89 PetscErrorCode Scatter_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt); /* u to ux_1,uy_1,ux_2,uy_2,...,u */ 90 PetscErrorCode Gather_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt); 91 92 static char help[] = ""; 93 94 int main(int argc, char **argv) 95 { 96 Vec x, x0; 97 Tao tao; 98 AppCtx user; 99 IS is_allstate, is_alldesign; 100 PetscInt lo, hi, hi2, lo2, ksp_old; 101 PetscInt ntests = 1; 102 PetscInt i; 103 PetscLogStage stages[1]; 104 105 PetscFunctionBeginUser; 106 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 107 user.mx = 32; 108 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "hyperbolic example", NULL); 109 PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL)); 110 user.nt = 16; 111 PetscCall(PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL)); 112 user.ndata = 64; 113 PetscCall(PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.ndata, NULL)); 114 user.alpha = 10.0; 115 PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL)); 116 user.T = 1.0 / 32.0; 117 PetscCall(PetscOptionsReal("-Tfinal", "Final time", "", user.T, &user.T, NULL)); 118 PetscCall(PetscOptionsInt("-ntests", "Number of times to repeat TaoSolve", "", ntests, &ntests, NULL)); 119 PetscOptionsEnd(); 120 121 user.m = user.mx * user.mx * user.nt; /* number of constraints */ 122 user.n = user.mx * user.mx * 3 * user.nt; /* number of variables */ 123 user.ht = user.T / user.nt; /* Time step */ 124 user.gamma = user.T * user.ht / (user.mx * user.mx); 125 126 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u)); 127 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y)); 128 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c)); 129 PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m)); 130 PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m)); 131 PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m)); 132 PetscCall(VecSetFromOptions(user.u)); 133 PetscCall(VecSetFromOptions(user.y)); 134 PetscCall(VecSetFromOptions(user.c)); 135 136 /* Create scatters for reduced spaces. 137 If the state vector y and design vector u are partitioned as 138 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors), 139 then the solution vector x is organized as 140 [y_1; u_1; y_2; u_2; ...; y_np; u_np]. 141 The index sets user.s_is and user.d_is correspond to the indices of the 142 state and design variables owned by the current processor. 143 */ 144 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 145 146 PetscCall(VecGetOwnershipRange(user.y, &lo, &hi)); 147 PetscCall(VecGetOwnershipRange(user.u, &lo2, &hi2)); 148 149 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_allstate)); 150 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + lo2, 1, &user.s_is)); 151 152 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, lo2, 1, &is_alldesign)); 153 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, hi + lo2, 1, &user.d_is)); 154 155 PetscCall(VecSetSizes(x, hi - lo + hi2 - lo2, user.n)); 156 PetscCall(VecSetFromOptions(x)); 157 158 PetscCall(VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter)); 159 PetscCall(VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter)); 160 PetscCall(ISDestroy(&is_alldesign)); 161 PetscCall(ISDestroy(&is_allstate)); 162 163 /* Create TAO solver and set desired solution method */ 164 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 165 PetscCall(TaoSetType(tao, TAOLCL)); 166 167 /* Set up initial vectors and matrices */ 168 PetscCall(HyperbolicInitialize(&user)); 169 170 PetscCall(Gather(x, user.y, user.state_scatter, user.u, user.design_scatter)); 171 PetscCall(VecDuplicate(x, &x0)); 172 PetscCall(VecCopy(x, x0)); 173 174 /* Set solution vector with an initial guess */ 175 PetscCall(TaoSetSolution(tao, x)); 176 PetscCall(TaoSetObjective(tao, FormFunction, &user)); 177 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user)); 178 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user)); 179 PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, &user)); 180 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user)); 181 PetscCall(TaoSetFromOptions(tao)); 182 PetscCall(TaoSetStateDesignIS(tao, user.s_is, user.d_is)); 183 184 /* SOLVE THE APPLICATION */ 185 PetscCall(PetscLogStageRegister("Trials", &stages[0])); 186 PetscCall(PetscLogStagePush(stages[0])); 187 user.ksp_its_initial = user.ksp_its; 188 ksp_old = user.ksp_its; 189 for (i = 0; i < ntests; i++) { 190 PetscCall(TaoSolve(tao)); 191 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp_old)); 192 PetscCall(VecCopy(x0, x)); 193 PetscCall(TaoSetSolution(tao, x)); 194 } 195 PetscCall(PetscLogStagePop()); 196 PetscCall(PetscBarrier((PetscObject)x)); 197 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations within initialization: ")); 198 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial)); 199 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Total KSP iterations over %" PetscInt_FMT " trial(s): ", ntests)); 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its)); 201 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations per trial: ")); 202 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initial) / ntests)); 203 204 PetscCall(TaoDestroy(&tao)); 205 PetscCall(VecDestroy(&x)); 206 PetscCall(VecDestroy(&x0)); 207 PetscCall(HyperbolicDestroy(&user)); 208 PetscCall(PetscFinalize()); 209 return 0; 210 } 211 /* ------------------------------------------------------------------- */ 212 /* 213 dwork = Qy - d 214 lwork = L*(u-ur).^2 215 f = 1/2 * (dwork.dork + alpha*y.lwork) 216 */ 217 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr) 218 { 219 PetscReal d1 = 0, d2 = 0; 220 AppCtx *user = (AppCtx *)ptr; 221 222 PetscFunctionBegin; 223 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 224 PetscCall(MatMult(user->Q, user->y, user->dwork)); 225 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 226 PetscCall(VecDot(user->dwork, user->dwork, &d1)); 227 228 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); 229 PetscCall(VecPointwiseMult(user->uwork, user->uwork, user->uwork)); 230 PetscCall(MatMult(user->L, user->uwork, user->lwork)); 231 PetscCall(VecDot(user->y, user->lwork, &d2)); 232 *f = 0.5 * (d1 + user->alpha * d2); 233 PetscFunctionReturn(PETSC_SUCCESS); 234 } 235 236 /* ------------------------------------------------------------------- */ 237 /* 238 state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2 239 design: g_d = alpha*(L'y).*(u-ur) 240 */ 241 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr) 242 { 243 AppCtx *user = (AppCtx *)ptr; 244 245 PetscFunctionBegin; 246 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 247 PetscCall(MatMult(user->Q, user->y, user->dwork)); 248 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 249 250 PetscCall(MatMult(user->QT, user->dwork, user->ywork)); 251 252 PetscCall(MatMult(user->LT, user->y, user->uwork)); 253 PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u)); 254 PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork)); 255 PetscCall(VecScale(user->uwork, user->alpha)); 256 257 PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork)); 258 PetscCall(MatMult(user->L, user->vwork, user->lwork)); 259 PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork)); 260 261 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 262 PetscFunctionReturn(PETSC_SUCCESS); 263 } 264 265 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) 266 { 267 PetscReal d1, d2; 268 AppCtx *user = (AppCtx *)ptr; 269 270 PetscFunctionBegin; 271 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 272 PetscCall(MatMult(user->Q, user->y, user->dwork)); 273 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 274 275 PetscCall(MatMult(user->QT, user->dwork, user->ywork)); 276 277 PetscCall(VecDot(user->dwork, user->dwork, &d1)); 278 279 PetscCall(MatMult(user->LT, user->y, user->uwork)); 280 PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u)); 281 PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork)); 282 PetscCall(VecScale(user->uwork, user->alpha)); 283 284 PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork)); 285 PetscCall(MatMult(user->L, user->vwork, user->lwork)); 286 PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork)); 287 288 PetscCall(VecDot(user->y, user->lwork, &d2)); 289 290 *f = 0.5 * (d1 + user->alpha * d2); 291 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 295 /* ------------------------------------------------------------------- */ 296 /* A 297 MatShell object 298 */ 299 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr) 300 { 301 PetscInt i; 302 AppCtx *user = (AppCtx *)ptr; 303 304 PetscFunctionBegin; 305 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 306 PetscCall(Scatter_yi(user->u, user->ui, user->ui_scatter, user->nt)); 307 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 308 for (i = 0; i < user->nt; i++) { 309 PetscCall(MatCopy(user->Divxy[0], user->C[i], SUBSET_NONZERO_PATTERN)); 310 PetscCall(MatCopy(user->Divxy[1], user->Cwork[i], SAME_NONZERO_PATTERN)); 311 312 PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i])); 313 PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i])); 314 PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], SUBSET_NONZERO_PATTERN)); 315 PetscCall(MatScale(user->C[i], user->ht)); 316 PetscCall(MatShift(user->C[i], 1.0)); 317 } 318 PetscFunctionReturn(PETSC_SUCCESS); 319 } 320 321 /* ------------------------------------------------------------------- */ 322 /* B */ 323 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr) 324 { 325 AppCtx *user = (AppCtx *)ptr; 326 327 PetscFunctionBegin; 328 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 329 PetscFunctionReturn(PETSC_SUCCESS); 330 } 331 332 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y) 333 { 334 PetscInt i; 335 AppCtx *user; 336 337 PetscFunctionBegin; 338 PetscCall(MatShellGetContext(J_shell, &user)); 339 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 340 user->block_index = 0; 341 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 342 343 for (i = 1; i < user->nt; i++) { 344 user->block_index = i; 345 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 346 PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1])); 347 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1])); 348 } 349 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 350 PetscFunctionReturn(PETSC_SUCCESS); 351 } 352 353 PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y) 354 { 355 PetscInt i; 356 AppCtx *user; 357 358 PetscFunctionBegin; 359 PetscCall(MatShellGetContext(J_shell, &user)); 360 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 361 362 for (i = 0; i < user->nt - 1; i++) { 363 user->block_index = i; 364 PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i])); 365 PetscCall(MatMult(user->M, user->yi[i + 1], user->ziwork[i + 1])); 366 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i + 1])); 367 } 368 369 i = user->nt - 1; 370 user->block_index = i; 371 PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i])); 372 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 373 PetscFunctionReturn(PETSC_SUCCESS); 374 } 375 376 PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y) 377 { 378 PetscInt i; 379 AppCtx *user; 380 381 PetscFunctionBegin; 382 PetscCall(MatShellGetContext(J_shell, &user)); 383 i = user->block_index; 384 PetscCall(VecPointwiseMult(user->uxiwork[i], X, user->uxi[i])); 385 PetscCall(VecPointwiseMult(user->uyiwork[i], X, user->uyi[i])); 386 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 387 PetscCall(MatMult(user->Div, user->uiwork[i], Y)); 388 PetscCall(VecAYPX(Y, user->ht, X)); 389 PetscFunctionReturn(PETSC_SUCCESS); 390 } 391 392 PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y) 393 { 394 PetscInt i; 395 AppCtx *user; 396 397 PetscFunctionBegin; 398 PetscCall(MatShellGetContext(J_shell, &user)); 399 i = user->block_index; 400 PetscCall(MatMult(user->Grad, X, user->uiwork[i])); 401 PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 402 PetscCall(VecPointwiseMult(user->uxiwork[i], user->uxi[i], user->uxiwork[i])); 403 PetscCall(VecPointwiseMult(user->uyiwork[i], user->uyi[i], user->uyiwork[i])); 404 PetscCall(VecWAXPY(Y, 1.0, user->uxiwork[i], user->uyiwork[i])); 405 PetscCall(VecAYPX(Y, user->ht, X)); 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408 409 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y) 410 { 411 PetscInt i; 412 AppCtx *user; 413 414 PetscFunctionBegin; 415 PetscCall(MatShellGetContext(J_shell, &user)); 416 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt)); 417 PetscCall(Scatter_uxi_uyi(X, user->uxiwork, user->uxi_scatter, user->uyiwork, user->uyi_scatter, user->nt)); 418 for (i = 0; i < user->nt; i++) { 419 PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i])); 420 PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i])); 421 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 422 PetscCall(MatMult(user->Div, user->uiwork[i], user->ziwork[i])); 423 PetscCall(VecScale(user->ziwork[i], user->ht)); 424 } 425 PetscCall(Gather_yi(Y, user->ziwork, user->yi_scatter, user->nt)); 426 PetscFunctionReturn(PETSC_SUCCESS); 427 } 428 429 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 430 { 431 PetscInt i; 432 AppCtx *user; 433 434 PetscFunctionBegin; 435 PetscCall(MatShellGetContext(J_shell, &user)); 436 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt)); 437 PetscCall(Scatter_yi(X, user->yiwork, user->yi_scatter, user->nt)); 438 for (i = 0; i < user->nt; i++) { 439 PetscCall(MatMult(user->Grad, user->yiwork[i], user->uiwork[i])); 440 PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 441 PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i])); 442 PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i])); 443 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i])); 444 PetscCall(VecScale(user->uiwork[i], user->ht)); 445 } 446 PetscCall(Gather_yi(Y, user->uiwork, user->ui_scatter, user->nt)); 447 PetscFunctionReturn(PETSC_SUCCESS); 448 } 449 450 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y) 451 { 452 PetscInt i; 453 AppCtx *user; 454 455 PetscFunctionBegin; 456 PetscCall(PCShellGetContext(PC_shell, &user)); 457 i = user->block_index; 458 PetscCheck(user->c_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed"); 459 PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y)); 460 PetscFunctionReturn(PETSC_SUCCESS); 461 } 462 463 PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y) 464 { 465 PetscInt i; 466 AppCtx *user; 467 468 PetscFunctionBegin; 469 PetscCall(PCShellGetContext(PC_shell, &user)); 470 471 i = user->block_index; 472 PetscCheck(user->c_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed"); 473 PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y)); 474 PetscFunctionReturn(PETSC_SUCCESS); 475 } 476 477 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y) 478 { 479 AppCtx *user; 480 PetscInt its, i; 481 482 PetscFunctionBegin; 483 PetscCall(MatShellGetContext(J_shell, &user)); 484 485 if (Y == user->ytrue) { 486 /* First solve is done using true solution to set up problem */ 487 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, PETSC_CURRENT, PETSC_CURRENT)); 488 } else { 489 PetscCall(KSPSetTolerances(user->solver, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT)); 490 } 491 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 492 PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 493 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 494 495 user->block_index = 0; 496 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0])); 497 498 PetscCall(KSPGetIterationNumber(user->solver, &its)); 499 user->ksp_its = user->ksp_its + its; 500 for (i = 1; i < user->nt; i++) { 501 PetscCall(MatMult(user->M, user->yiwork[i - 1], user->ziwork[i - 1])); 502 PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i - 1])); 503 user->block_index = i; 504 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 505 506 PetscCall(KSPGetIterationNumber(user->solver, &its)); 507 user->ksp_its = user->ksp_its + its; 508 } 509 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 510 PetscFunctionReturn(PETSC_SUCCESS); 511 } 512 513 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y) 514 { 515 AppCtx *user; 516 PetscInt its, i; 517 518 PetscFunctionBegin; 519 PetscCall(MatShellGetContext(J_shell, &user)); 520 521 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt)); 522 PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 523 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 524 525 i = user->nt - 1; 526 user->block_index = i; 527 PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i])); 528 529 PetscCall(KSPGetIterationNumber(user->solver, &its)); 530 user->ksp_its = user->ksp_its + its; 531 532 for (i = user->nt - 2; i >= 0; i--) { 533 PetscCall(MatMult(user->M, user->yiwork[i + 1], user->ziwork[i + 1])); 534 PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i + 1])); 535 user->block_index = i; 536 PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i])); 537 538 PetscCall(KSPGetIterationNumber(user->solver, &its)); 539 user->ksp_its = user->ksp_its + its; 540 } 541 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt)); 542 PetscFunctionReturn(PETSC_SUCCESS); 543 } 544 545 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell) 546 { 547 AppCtx *user; 548 549 PetscFunctionBegin; 550 PetscCall(MatShellGetContext(J_shell, &user)); 551 552 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell)); 553 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult)); 554 PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate)); 555 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTranspose)); 556 PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagonal)); 557 PetscFunctionReturn(PETSC_SUCCESS); 558 } 559 560 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X) 561 { 562 AppCtx *user; 563 564 PetscFunctionBegin; 565 PetscCall(MatShellGetContext(J_shell, &user)); 566 PetscCall(VecCopy(user->js_diag, X)); 567 PetscFunctionReturn(PETSC_SUCCESS); 568 } 569 570 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr) 571 { 572 /* con = Ay - q, A = [C(u1) 0 0 ... 0; 573 -M C(u2) 0 ... 0; 574 0 -M C(u3) ... 0; 575 ... ; 576 0 ... -M C(u_nt)] 577 C(u) = eye + ht*Div*[diag(u1); diag(u2)] */ 578 PetscInt i; 579 AppCtx *user = (AppCtx *)ptr; 580 581 PetscFunctionBegin; 582 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 583 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt)); 584 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 585 586 user->block_index = 0; 587 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 588 589 for (i = 1; i < user->nt; i++) { 590 user->block_index = i; 591 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 592 PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1])); 593 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1])); 594 } 595 596 PetscCall(Gather_yi(C, user->yiwork, user->yi_scatter, user->nt)); 597 PetscCall(VecAXPY(C, -1.0, user->q)); 598 PetscFunctionReturn(PETSC_SUCCESS); 599 } 600 601 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 602 { 603 PetscFunctionBegin; 604 PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 605 PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 606 PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 607 PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 608 PetscFunctionReturn(PETSC_SUCCESS); 609 } 610 611 PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt) 612 { 613 PetscInt i; 614 615 PetscFunctionBegin; 616 for (i = 0; i < nt; i++) { 617 PetscCall(VecScatterBegin(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD)); 618 PetscCall(VecScatterEnd(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD)); 619 PetscCall(VecScatterBegin(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD)); 620 PetscCall(VecScatterEnd(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD)); 621 } 622 PetscFunctionReturn(PETSC_SUCCESS); 623 } 624 625 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 626 { 627 PetscFunctionBegin; 628 PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 629 PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 630 PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 631 PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 632 PetscFunctionReturn(PETSC_SUCCESS); 633 } 634 635 PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt) 636 { 637 PetscInt i; 638 639 PetscFunctionBegin; 640 for (i = 0; i < nt; i++) { 641 PetscCall(VecScatterBegin(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 642 PetscCall(VecScatterEnd(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 643 PetscCall(VecScatterBegin(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 644 PetscCall(VecScatterEnd(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE)); 645 } 646 PetscFunctionReturn(PETSC_SUCCESS); 647 } 648 649 PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 650 { 651 PetscInt i; 652 653 PetscFunctionBegin; 654 for (i = 0; i < nt; i++) { 655 PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 656 PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 657 } 658 PetscFunctionReturn(PETSC_SUCCESS); 659 } 660 661 PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 662 { 663 PetscInt i; 664 665 PetscFunctionBegin; 666 for (i = 0; i < nt; i++) { 667 PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 668 PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 669 } 670 PetscFunctionReturn(PETSC_SUCCESS); 671 } 672 673 PetscErrorCode HyperbolicInitialize(AppCtx *user) 674 { 675 PetscInt n, i, j, linear_index, istart, iend, iblock, lo, hi; 676 Vec XX, YY, XXwork, YYwork, yi, uxi, ui, bc; 677 PetscReal h, sum; 678 PetscScalar hinv, neg_hinv, quarter = 0.25, one = 1.0, half_hinv, neg_half_hinv; 679 PetscScalar vx, vy, zero = 0.0; 680 IS is_from_y, is_to_yi, is_from_u, is_to_uxi, is_to_uyi; 681 682 PetscFunctionBegin; 683 user->jformed = PETSC_FALSE; 684 user->c_formed = PETSC_FALSE; 685 686 user->ksp_its = 0; 687 user->ksp_its_initial = 0; 688 689 n = user->mx * user->mx; 690 691 h = 1.0 / user->mx; 692 hinv = user->mx; 693 neg_hinv = -hinv; 694 half_hinv = hinv / 2.0; 695 neg_half_hinv = neg_hinv / 2.0; 696 697 /* Generate Grad matrix */ 698 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad)); 699 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, 2 * n, n)); 700 PetscCall(MatSetFromOptions(user->Grad)); 701 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 3, NULL, 3, NULL)); 702 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 3, NULL)); 703 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend)); 704 705 for (i = istart; i < iend; i++) { 706 if (i < n) { 707 iblock = i / user->mx; 708 j = iblock * user->mx + ((i + user->mx - 1) % user->mx); 709 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 710 j = iblock * user->mx + ((i + 1) % user->mx); 711 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 712 } 713 if (i >= n) { 714 j = (i - user->mx) % n; 715 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 716 j = (j + 2 * user->mx) % n; 717 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 718 } 719 } 720 721 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY)); 722 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY)); 723 724 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[0])); 725 PetscCall(MatSetSizes(user->Gradxy[0], PETSC_DECIDE, PETSC_DECIDE, n, n)); 726 PetscCall(MatSetFromOptions(user->Gradxy[0])); 727 PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[0], 3, NULL, 3, NULL)); 728 PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[0], 3, NULL)); 729 PetscCall(MatGetOwnershipRange(user->Gradxy[0], &istart, &iend)); 730 731 for (i = istart; i < iend; i++) { 732 iblock = i / user->mx; 733 j = iblock * user->mx + ((i + user->mx - 1) % user->mx); 734 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 735 j = iblock * user->mx + ((i + 1) % user->mx); 736 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 737 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &i, &zero, INSERT_VALUES)); 738 } 739 PetscCall(MatAssemblyBegin(user->Gradxy[0], MAT_FINAL_ASSEMBLY)); 740 PetscCall(MatAssemblyEnd(user->Gradxy[0], MAT_FINAL_ASSEMBLY)); 741 742 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[1])); 743 PetscCall(MatSetSizes(user->Gradxy[1], PETSC_DECIDE, PETSC_DECIDE, n, n)); 744 PetscCall(MatSetFromOptions(user->Gradxy[1])); 745 PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[1], 3, NULL, 3, NULL)); 746 PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[1], 3, NULL)); 747 PetscCall(MatGetOwnershipRange(user->Gradxy[1], &istart, &iend)); 748 749 for (i = istart; i < iend; i++) { 750 j = (i + n - user->mx) % n; 751 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &half_hinv, INSERT_VALUES)); 752 j = (j + 2 * user->mx) % n; 753 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES)); 754 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &i, &zero, INSERT_VALUES)); 755 } 756 PetscCall(MatAssemblyBegin(user->Gradxy[1], MAT_FINAL_ASSEMBLY)); 757 PetscCall(MatAssemblyEnd(user->Gradxy[1], MAT_FINAL_ASSEMBLY)); 758 759 /* Generate Div matrix */ 760 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div)); 761 PetscCall(MatTranspose(user->Gradxy[0], MAT_INITIAL_MATRIX, &user->Divxy[0])); 762 PetscCall(MatTranspose(user->Gradxy[1], MAT_INITIAL_MATRIX, &user->Divxy[1])); 763 764 /* Off-diagonal averaging matrix */ 765 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->M)); 766 PetscCall(MatSetSizes(user->M, PETSC_DECIDE, PETSC_DECIDE, n, n)); 767 PetscCall(MatSetFromOptions(user->M)); 768 PetscCall(MatMPIAIJSetPreallocation(user->M, 4, NULL, 4, NULL)); 769 PetscCall(MatSeqAIJSetPreallocation(user->M, 4, NULL)); 770 PetscCall(MatGetOwnershipRange(user->M, &istart, &iend)); 771 772 for (i = istart; i < iend; i++) { 773 /* kron(Id,Av) */ 774 iblock = i / user->mx; 775 j = iblock * user->mx + ((i + user->mx - 1) % user->mx); 776 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 777 j = iblock * user->mx + ((i + 1) % user->mx); 778 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 779 780 /* kron(Av,Id) */ 781 j = (i + user->mx) % n; 782 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 783 j = (i + n - user->mx) % n; 784 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES)); 785 } 786 PetscCall(MatAssemblyBegin(user->M, MAT_FINAL_ASSEMBLY)); 787 PetscCall(MatAssemblyEnd(user->M, MAT_FINAL_ASSEMBLY)); 788 789 /* Generate 2D grid */ 790 PetscCall(VecCreate(PETSC_COMM_WORLD, &XX)); 791 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q)); 792 PetscCall(VecSetSizes(XX, PETSC_DECIDE, n)); 793 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt)); 794 PetscCall(VecSetFromOptions(XX)); 795 PetscCall(VecSetFromOptions(user->q)); 796 797 PetscCall(VecDuplicate(XX, &YY)); 798 PetscCall(VecDuplicate(XX, &XXwork)); 799 PetscCall(VecDuplicate(XX, &YYwork)); 800 PetscCall(VecDuplicate(XX, &user->d)); 801 PetscCall(VecDuplicate(XX, &user->dwork)); 802 803 PetscCall(VecGetOwnershipRange(XX, &istart, &iend)); 804 for (linear_index = istart; linear_index < iend; linear_index++) { 805 i = linear_index % user->mx; 806 j = (linear_index - i) / user->mx; 807 vx = h * (i + 0.5); 808 vy = h * (j + 0.5); 809 PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES)); 810 PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES)); 811 } 812 813 PetscCall(VecAssemblyBegin(XX)); 814 PetscCall(VecAssemblyEnd(XX)); 815 PetscCall(VecAssemblyBegin(YY)); 816 PetscCall(VecAssemblyEnd(YY)); 817 818 /* Compute final density function yT 819 yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2)) 820 yT = yT / (h^2*sum(yT)) */ 821 PetscCall(VecCopy(XX, XXwork)); 822 PetscCall(VecCopy(YY, YYwork)); 823 824 PetscCall(VecShift(XXwork, -0.25)); 825 PetscCall(VecShift(YYwork, -0.25)); 826 827 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 828 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 829 830 PetscCall(VecCopy(XXwork, user->dwork)); 831 PetscCall(VecAXPY(user->dwork, 1.0, YYwork)); 832 PetscCall(VecScale(user->dwork, -30.0)); 833 PetscCall(VecExp(user->dwork)); 834 PetscCall(VecCopy(user->dwork, user->d)); 835 836 PetscCall(VecCopy(XX, XXwork)); 837 PetscCall(VecCopy(YY, YYwork)); 838 839 PetscCall(VecShift(XXwork, -0.75)); 840 PetscCall(VecShift(YYwork, -0.75)); 841 842 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 843 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 844 845 PetscCall(VecCopy(XXwork, user->dwork)); 846 PetscCall(VecAXPY(user->dwork, 1.0, YYwork)); 847 PetscCall(VecScale(user->dwork, -30.0)); 848 PetscCall(VecExp(user->dwork)); 849 850 PetscCall(VecAXPY(user->d, 1.0, user->dwork)); 851 PetscCall(VecShift(user->d, 1.0)); 852 PetscCall(VecSum(user->d, &sum)); 853 PetscCall(VecScale(user->d, 1.0 / (h * h * sum))); 854 855 /* Initial conditions of forward problem */ 856 PetscCall(VecDuplicate(XX, &bc)); 857 PetscCall(VecCopy(XX, XXwork)); 858 PetscCall(VecCopy(YY, YYwork)); 859 860 PetscCall(VecShift(XXwork, -0.5)); 861 PetscCall(VecShift(YYwork, -0.5)); 862 863 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 864 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 865 866 PetscCall(VecWAXPY(bc, 1.0, XXwork, YYwork)); 867 PetscCall(VecScale(bc, -50.0)); 868 PetscCall(VecExp(bc)); 869 PetscCall(VecShift(bc, 1.0)); 870 PetscCall(VecSum(bc, &sum)); 871 PetscCall(VecScale(bc, 1.0 / (h * h * sum))); 872 873 /* Create scatter from y to y_1,y_2,...,y_nt */ 874 /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */ 875 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->yi_scatter)); 876 PetscCall(VecCreate(PETSC_COMM_WORLD, &yi)); 877 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx)); 878 PetscCall(VecSetFromOptions(yi)); 879 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi)); 880 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork)); 881 PetscCall(VecDuplicateVecs(yi, user->nt, &user->ziwork)); 882 for (i = 0; i < user->nt; i++) { 883 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi)); 884 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi)); 885 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + i * user->mx * user->mx, 1, &is_from_y)); 886 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i])); 887 PetscCall(ISDestroy(&is_to_yi)); 888 PetscCall(ISDestroy(&is_from_y)); 889 } 890 891 /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */ 892 /* TODO: reorder for better parallelism */ 893 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uxi_scatter)); 894 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uyi_scatter)); 895 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->ux_scatter)); 896 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uy_scatter)); 897 PetscCall(PetscMalloc1(2 * user->nt * user->mx * user->mx, &user->ui_scatter)); 898 PetscCall(VecCreate(PETSC_COMM_WORLD, &uxi)); 899 PetscCall(VecCreate(PETSC_COMM_WORLD, &ui)); 900 PetscCall(VecSetSizes(uxi, PETSC_DECIDE, user->mx * user->mx)); 901 PetscCall(VecSetSizes(ui, PETSC_DECIDE, 2 * user->mx * user->mx)); 902 PetscCall(VecSetFromOptions(uxi)); 903 PetscCall(VecSetFromOptions(ui)); 904 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxi)); 905 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyi)); 906 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxiwork)); 907 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyiwork)); 908 PetscCall(VecDuplicateVecs(ui, user->nt, &user->ui)); 909 PetscCall(VecDuplicateVecs(ui, user->nt, &user->uiwork)); 910 for (i = 0; i < user->nt; i++) { 911 PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi)); 912 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi)); 913 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u)); 914 PetscCall(VecScatterCreate(user->u, is_from_u, user->uxi[i], is_to_uxi, &user->uxi_scatter[i])); 915 916 PetscCall(ISDestroy(&is_to_uxi)); 917 PetscCall(ISDestroy(&is_from_u)); 918 919 PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi)); 920 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi)); 921 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + (2 * i + 1) * user->mx * user->mx, 1, &is_from_u)); 922 PetscCall(VecScatterCreate(user->u, is_from_u, user->uyi[i], is_to_uyi, &user->uyi_scatter[i])); 923 924 PetscCall(ISDestroy(&is_to_uyi)); 925 PetscCall(ISDestroy(&is_from_u)); 926 927 PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi)); 928 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi)); 929 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_from_u)); 930 PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uxi[i], is_to_uxi, &user->ux_scatter[i])); 931 932 PetscCall(ISDestroy(&is_to_uxi)); 933 PetscCall(ISDestroy(&is_from_u)); 934 935 PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi)); 936 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi)); 937 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + user->mx * user->mx, 1, &is_from_u)); 938 PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uyi[i], is_to_uyi, &user->uy_scatter[i])); 939 940 PetscCall(ISDestroy(&is_to_uyi)); 941 PetscCall(ISDestroy(&is_from_u)); 942 943 PetscCall(VecGetOwnershipRange(user->ui[i], &lo, &hi)); 944 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi)); 945 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u)); 946 PetscCall(VecScatterCreate(user->u, is_from_u, user->ui[i], is_to_uxi, &user->ui_scatter[i])); 947 948 PetscCall(ISDestroy(&is_to_uxi)); 949 PetscCall(ISDestroy(&is_from_u)); 950 } 951 952 /* RHS of forward problem */ 953 PetscCall(MatMult(user->M, bc, user->yiwork[0])); 954 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0)); 955 PetscCall(Gather_yi(user->q, user->yiwork, user->yi_scatter, user->nt)); 956 957 /* Compute true velocity field utrue */ 958 PetscCall(VecDuplicate(user->u, &user->utrue)); 959 for (i = 0; i < user->nt; i++) { 960 PetscCall(VecCopy(YY, user->uxi[i])); 961 PetscCall(VecScale(user->uxi[i], 150.0 * i * user->ht)); 962 PetscCall(VecCopy(XX, user->uyi[i])); 963 PetscCall(VecShift(user->uyi[i], -10.0)); 964 PetscCall(VecScale(user->uyi[i], 15.0 * i * user->ht)); 965 } 966 PetscCall(Gather_uxi_uyi(user->utrue, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 967 968 /* Initial guess and reference model */ 969 PetscCall(VecDuplicate(user->utrue, &user->ur)); 970 for (i = 0; i < user->nt; i++) { 971 PetscCall(VecCopy(XX, user->uxi[i])); 972 PetscCall(VecShift(user->uxi[i], i * user->ht)); 973 PetscCall(VecCopy(YY, user->uyi[i])); 974 PetscCall(VecShift(user->uyi[i], -i * user->ht)); 975 } 976 PetscCall(Gather_uxi_uyi(user->ur, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 977 978 /* Generate regularization matrix L */ 979 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->LT)); 980 PetscCall(MatSetSizes(user->LT, PETSC_DECIDE, PETSC_DECIDE, 2 * n * user->nt, n * user->nt)); 981 PetscCall(MatSetFromOptions(user->LT)); 982 PetscCall(MatMPIAIJSetPreallocation(user->LT, 1, NULL, 1, NULL)); 983 PetscCall(MatSeqAIJSetPreallocation(user->LT, 1, NULL)); 984 PetscCall(MatGetOwnershipRange(user->LT, &istart, &iend)); 985 986 for (i = istart; i < iend; i++) { 987 iblock = (i + n) / (2 * n); 988 j = i - iblock * n; 989 PetscCall(MatSetValues(user->LT, 1, &i, 1, &j, &user->gamma, INSERT_VALUES)); 990 } 991 992 PetscCall(MatAssemblyBegin(user->LT, MAT_FINAL_ASSEMBLY)); 993 PetscCall(MatAssemblyEnd(user->LT, MAT_FINAL_ASSEMBLY)); 994 995 PetscCall(MatTranspose(user->LT, MAT_INITIAL_MATRIX, &user->L)); 996 997 /* Build work vectors and matrices */ 998 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork)); 999 PetscCall(VecSetType(user->lwork, VECMPI)); 1000 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, user->m)); 1001 PetscCall(VecSetFromOptions(user->lwork)); 1002 1003 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork)); 1004 1005 PetscCall(VecDuplicate(user->y, &user->ywork)); 1006 PetscCall(VecDuplicate(user->u, &user->uwork)); 1007 PetscCall(VecDuplicate(user->u, &user->vwork)); 1008 PetscCall(VecDuplicate(user->u, &user->js_diag)); 1009 PetscCall(VecDuplicate(user->c, &user->cwork)); 1010 1011 /* Create matrix-free shell user->Js for computing A*x */ 1012 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->Js)); 1013 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult)); 1014 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate)); 1015 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTranspose)); 1016 PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagonal)); 1017 1018 /* Diagonal blocks of user->Js */ 1019 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlock)); 1020 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockMult)); 1021 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBlockMultTranspose)); 1022 1023 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U, 1024 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular. 1025 This is an SOR preconditioner for user->JsBlock. */ 1026 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlockPrec)); 1027 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockPrecMult)); 1028 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBlockPrecMultTranspose)); 1029 1030 /* Create a matrix-free shell user->Jd for computing B*x */ 1031 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->n - user->m, user, &user->Jd)); 1032 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (PetscErrorCodeFn *)DesignMatMult)); 1033 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)DesignMatMultTranspose)); 1034 1035 /* User-defined routines for computing user->Js\x and user->Js^T\x*/ 1036 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsInv)); 1037 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (PetscErrorCodeFn *)StateMatInvMult)); 1038 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatInvTransposeMult)); 1039 1040 /* Build matrices for SOR preconditioner */ 1041 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt)); 1042 PetscCall(PetscMalloc1(5 * n, &user->C)); 1043 PetscCall(PetscMalloc1(2 * n, &user->Cwork)); 1044 for (i = 0; i < user->nt; i++) { 1045 PetscCall(MatDuplicate(user->Divxy[0], MAT_COPY_VALUES, &user->C[i])); 1046 PetscCall(MatDuplicate(user->Divxy[1], MAT_COPY_VALUES, &user->Cwork[i])); 1047 1048 PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i])); 1049 PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i])); 1050 PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], DIFFERENT_NONZERO_PATTERN)); 1051 PetscCall(MatScale(user->C[i], user->ht)); 1052 PetscCall(MatShift(user->C[i], 1.0)); 1053 } 1054 1055 /* Solver options and tolerances */ 1056 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver)); 1057 PetscCall(KSPSetType(user->solver, KSPGMRES)); 1058 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec)); 1059 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500)); 1060 /* PetscCall(KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500)); */ 1061 PetscCall(KSPGetPC(user->solver, &user->prec)); 1062 PetscCall(PCSetType(user->prec, PCSHELL)); 1063 1064 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult)); 1065 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMultTranspose)); 1066 PetscCall(PCShellSetContext(user->prec, user)); 1067 1068 /* Compute true state function yt given ut */ 1069 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue)); 1070 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt)); 1071 PetscCall(VecSetFromOptions(user->ytrue)); 1072 user->c_formed = PETSC_TRUE; 1073 PetscCall(VecCopy(user->utrue, user->u)); /* Set u=utrue temporarily for StateMatInv */ 1074 PetscCall(VecSet(user->ytrue, 0.0)); /* Initial guess */ 1075 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue)); 1076 PetscCall(VecCopy(user->ur, user->u)); /* Reset u=ur */ 1077 1078 /* Initial guess y0 for state given u0 */ 1079 PetscCall(StateMatInvMult(user->Js, user->q, user->y)); 1080 1081 /* Data discretization */ 1082 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Q)); 1083 PetscCall(MatSetSizes(user->Q, PETSC_DECIDE, PETSC_DECIDE, user->mx * user->mx, user->m)); 1084 PetscCall(MatSetFromOptions(user->Q)); 1085 PetscCall(MatMPIAIJSetPreallocation(user->Q, 0, NULL, 1, NULL)); 1086 PetscCall(MatSeqAIJSetPreallocation(user->Q, 1, NULL)); 1087 1088 PetscCall(MatGetOwnershipRange(user->Q, &istart, &iend)); 1089 1090 for (i = istart; i < iend; i++) { 1091 j = i + user->m - user->mx * user->mx; 1092 PetscCall(MatSetValues(user->Q, 1, &i, 1, &j, &one, INSERT_VALUES)); 1093 } 1094 1095 PetscCall(MatAssemblyBegin(user->Q, MAT_FINAL_ASSEMBLY)); 1096 PetscCall(MatAssemblyEnd(user->Q, MAT_FINAL_ASSEMBLY)); 1097 1098 PetscCall(MatTranspose(user->Q, MAT_INITIAL_MATRIX, &user->QT)); 1099 1100 PetscCall(VecDestroy(&XX)); 1101 PetscCall(VecDestroy(&YY)); 1102 PetscCall(VecDestroy(&XXwork)); 1103 PetscCall(VecDestroy(&YYwork)); 1104 PetscCall(VecDestroy(&yi)); 1105 PetscCall(VecDestroy(&uxi)); 1106 PetscCall(VecDestroy(&ui)); 1107 PetscCall(VecDestroy(&bc)); 1108 1109 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */ 1110 PetscCall(KSPSetFromOptions(user->solver)); 1111 PetscFunctionReturn(PETSC_SUCCESS); 1112 } 1113 1114 PetscErrorCode HyperbolicDestroy(AppCtx *user) 1115 { 1116 PetscInt i; 1117 1118 PetscFunctionBegin; 1119 PetscCall(MatDestroy(&user->Q)); 1120 PetscCall(MatDestroy(&user->QT)); 1121 PetscCall(MatDestroy(&user->Div)); 1122 PetscCall(MatDestroy(&user->Divwork)); 1123 PetscCall(MatDestroy(&user->Grad)); 1124 PetscCall(MatDestroy(&user->L)); 1125 PetscCall(MatDestroy(&user->LT)); 1126 PetscCall(KSPDestroy(&user->solver)); 1127 PetscCall(MatDestroy(&user->Js)); 1128 PetscCall(MatDestroy(&user->Jd)); 1129 PetscCall(MatDestroy(&user->JsBlockPrec)); 1130 PetscCall(MatDestroy(&user->JsInv)); 1131 PetscCall(MatDestroy(&user->JsBlock)); 1132 PetscCall(MatDestroy(&user->Divxy[0])); 1133 PetscCall(MatDestroy(&user->Divxy[1])); 1134 PetscCall(MatDestroy(&user->Gradxy[0])); 1135 PetscCall(MatDestroy(&user->Gradxy[1])); 1136 PetscCall(MatDestroy(&user->M)); 1137 for (i = 0; i < user->nt; i++) { 1138 PetscCall(MatDestroy(&user->C[i])); 1139 PetscCall(MatDestroy(&user->Cwork[i])); 1140 } 1141 PetscCall(PetscFree(user->C)); 1142 PetscCall(PetscFree(user->Cwork)); 1143 PetscCall(VecDestroy(&user->u)); 1144 PetscCall(VecDestroy(&user->uwork)); 1145 PetscCall(VecDestroy(&user->vwork)); 1146 PetscCall(VecDestroy(&user->utrue)); 1147 PetscCall(VecDestroy(&user->y)); 1148 PetscCall(VecDestroy(&user->ywork)); 1149 PetscCall(VecDestroy(&user->ytrue)); 1150 PetscCall(VecDestroyVecs(user->nt, &user->yi)); 1151 PetscCall(VecDestroyVecs(user->nt, &user->yiwork)); 1152 PetscCall(VecDestroyVecs(user->nt, &user->ziwork)); 1153 PetscCall(VecDestroyVecs(user->nt, &user->uxi)); 1154 PetscCall(VecDestroyVecs(user->nt, &user->uyi)); 1155 PetscCall(VecDestroyVecs(user->nt, &user->uxiwork)); 1156 PetscCall(VecDestroyVecs(user->nt, &user->uyiwork)); 1157 PetscCall(VecDestroyVecs(user->nt, &user->ui)); 1158 PetscCall(VecDestroyVecs(user->nt, &user->uiwork)); 1159 PetscCall(VecDestroy(&user->c)); 1160 PetscCall(VecDestroy(&user->cwork)); 1161 PetscCall(VecDestroy(&user->ur)); 1162 PetscCall(VecDestroy(&user->q)); 1163 PetscCall(VecDestroy(&user->d)); 1164 PetscCall(VecDestroy(&user->dwork)); 1165 PetscCall(VecDestroy(&user->lwork)); 1166 PetscCall(VecDestroy(&user->js_diag)); 1167 PetscCall(ISDestroy(&user->s_is)); 1168 PetscCall(ISDestroy(&user->d_is)); 1169 PetscCall(VecScatterDestroy(&user->state_scatter)); 1170 PetscCall(VecScatterDestroy(&user->design_scatter)); 1171 for (i = 0; i < user->nt; i++) { 1172 PetscCall(VecScatterDestroy(&user->uxi_scatter[i])); 1173 PetscCall(VecScatterDestroy(&user->uyi_scatter[i])); 1174 PetscCall(VecScatterDestroy(&user->ux_scatter[i])); 1175 PetscCall(VecScatterDestroy(&user->uy_scatter[i])); 1176 PetscCall(VecScatterDestroy(&user->ui_scatter[i])); 1177 PetscCall(VecScatterDestroy(&user->yi_scatter[i])); 1178 } 1179 PetscCall(PetscFree(user->uxi_scatter)); 1180 PetscCall(PetscFree(user->uyi_scatter)); 1181 PetscCall(PetscFree(user->ux_scatter)); 1182 PetscCall(PetscFree(user->uy_scatter)); 1183 PetscCall(PetscFree(user->ui_scatter)); 1184 PetscCall(PetscFree(user->yi_scatter)); 1185 PetscFunctionReturn(PETSC_SUCCESS); 1186 } 1187 1188 PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr) 1189 { 1190 Vec X; 1191 PetscReal unorm, ynorm; 1192 AppCtx *user = (AppCtx *)ptr; 1193 1194 PetscFunctionBegin; 1195 PetscCall(TaoGetSolution(tao, &X)); 1196 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 1197 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue)); 1198 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); 1199 PetscCall(VecNorm(user->uwork, NORM_2, &unorm)); 1200 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm)); 1201 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm)); 1202 PetscFunctionReturn(PETSC_SUCCESS); 1203 } 1204 1205 /*TEST 1206 1207 build: 1208 requires: !complex 1209 1210 test: 1211 requires: !single 1212 args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5 1213 1214 test: 1215 suffix: guess_pod 1216 requires: !single 1217 args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5 1218 1219 TEST*/ 1220