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