1 #include <petsc/private/taoimpl.h> 2 3 typedef struct { 4 PetscInt n; /* Number of variables */ 5 PetscInt m; /* Number of constraints per time step */ 6 PetscInt mx; /* grid points in each direction */ 7 PetscInt nt; /* Number of time steps; as of now, must be divisible by 8 */ 8 PetscInt ndata; /* Number of data points per sample */ 9 PetscInt ns; /* Number of samples */ 10 PetscInt *sample_times; /* Times of samples */ 11 IS s_is; 12 IS d_is; 13 14 VecScatter state_scatter; 15 VecScatter design_scatter; 16 VecScatter *yi_scatter; 17 VecScatter *di_scatter; 18 19 Mat Js, Jd, JsBlockPrec, JsInv, JsBlock; 20 PetscBool jformed, dsg_formed; 21 22 PetscReal alpha; /* Regularization parameter */ 23 PetscReal beta; /* Weight attributed to ||u||^2 in regularization functional */ 24 PetscReal noise; /* Amount of noise to add to data */ 25 PetscReal ht; /* Time step */ 26 27 Mat Qblock, QblockT; 28 Mat L, LT; 29 Mat Div, Divwork; 30 Mat Grad; 31 Mat Av, Avwork, AvT; 32 Mat DSG; 33 Vec q; 34 Vec ur; /* reference */ 35 36 Vec d; 37 Vec dwork; 38 Vec *di; 39 40 Vec y; /* state variables */ 41 Vec ywork; 42 43 Vec ytrue; 44 Vec *yi, *yiwork; 45 46 Vec u; /* design variables */ 47 Vec uwork; 48 49 Vec utrue; 50 Vec js_diag; 51 Vec c; /* constraint vector */ 52 Vec cwork; 53 54 Vec lwork; 55 Vec S; 56 Vec Rwork, Swork, Twork; 57 Vec Av_u; 58 59 KSP solver; 60 PC prec; 61 62 PetscInt ksp_its; 63 PetscInt ksp_its_initial; 64 } AppCtx; 65 66 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *); 67 PetscErrorCode FormGradient(Tao, Vec, Vec, void *); 68 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *); 69 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void *); 70 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void *); 71 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *); 72 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *); 73 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 74 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 75 PetscErrorCode ParabolicInitialize(AppCtx *user); 76 PetscErrorCode ParabolicDestroy(AppCtx *user); 77 PetscErrorCode ParabolicMonitor(Tao, void *); 78 79 PetscErrorCode StateMatMult(Mat, Vec, Vec); 80 PetscErrorCode StateMatBlockMult(Mat, Vec, Vec); 81 PetscErrorCode StateMatMultTranspose(Mat, Vec, Vec); 82 PetscErrorCode StateMatGetDiagonal(Mat, Vec); 83 PetscErrorCode StateMatDuplicate(Mat, MatDuplicateOption, Mat *); 84 PetscErrorCode StateMatInvMult(Mat, Vec, Vec); 85 PetscErrorCode StateMatInvTransposeMult(Mat, Vec, Vec); 86 PetscErrorCode StateMatBlockPrecMult(PC, Vec, Vec); 87 88 PetscErrorCode DesignMatMult(Mat, Vec, Vec); 89 PetscErrorCode DesignMatMultTranspose(Mat, Vec, Vec); 90 91 PetscErrorCode Gather_i(Vec, Vec *, VecScatter *, PetscInt); 92 PetscErrorCode Scatter_i(Vec, Vec *, VecScatter *, PetscInt); 93 94 static char help[] = ""; 95 96 int main(int argc, char **argv) 97 { 98 Vec x, x0; 99 Tao tao; 100 AppCtx user; 101 IS is_allstate, is_alldesign; 102 PetscInt lo, hi, hi2, lo2, ksp_old; 103 PetscInt ntests = 1; 104 PetscInt i; 105 PetscLogStage stages[1]; 106 107 PetscFunctionBeginUser; 108 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 109 user.mx = 8; 110 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "parabolic example", NULL); 111 PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL)); 112 user.nt = 8; 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.ns = 8; 117 PetscCall(PetscOptionsInt("-ns", "Number of samples", "", user.ns, &user.ns, NULL)); 118 user.alpha = 1.0; 119 PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL)); 120 user.beta = 0.01; 121 PetscCall(PetscOptionsReal("-beta", "Weight attributed to ||u||^2 in regularization functional", "", user.beta, &user.beta, NULL)); 122 user.noise = 0.01; 123 PetscCall(PetscOptionsReal("-noise", "Amount of noise to add to data", "", user.noise, &user.noise, NULL)); 124 PetscCall(PetscOptionsInt("-ntests", "Number of times to repeat TaoSolve", "", ntests, &ntests, NULL)); 125 PetscOptionsEnd(); 126 127 user.m = user.mx * user.mx * user.mx; /* number of constraints per time step */ 128 user.n = user.m * (user.nt + 1); /* number of variables */ 129 user.ht = (PetscReal)1 / user.nt; 130 131 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u)); 132 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y)); 133 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c)); 134 PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m * user.nt)); 135 PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m * user.nt)); 136 PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m * user.nt)); 137 PetscCall(VecSetFromOptions(user.u)); 138 PetscCall(VecSetFromOptions(user.y)); 139 PetscCall(VecSetFromOptions(user.c)); 140 141 /* Create scatters for reduced spaces. 142 If the state vector y and design vector u are partitioned as 143 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors), 144 then the solution vector x is organized as 145 [y_1; u_1; y_2; u_2; ...; y_np; u_np]. 146 The index sets user.s_is and user.d_is correspond to the indices of the 147 state and design variables owned by the current processor. 148 */ 149 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 150 151 PetscCall(VecGetOwnershipRange(user.y, &lo, &hi)); 152 PetscCall(VecGetOwnershipRange(user.u, &lo2, &hi2)); 153 154 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_allstate)); 155 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + lo2, 1, &user.s_is)); 156 157 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, lo2, 1, &is_alldesign)); 158 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, hi + lo2, 1, &user.d_is)); 159 160 PetscCall(VecSetSizes(x, hi - lo + hi2 - lo2, user.n)); 161 PetscCall(VecSetFromOptions(x)); 162 163 PetscCall(VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter)); 164 PetscCall(VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter)); 165 PetscCall(ISDestroy(&is_alldesign)); 166 PetscCall(ISDestroy(&is_allstate)); 167 168 /* Create TAO solver and set desired solution method */ 169 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 170 PetscCall(TaoSetType(tao, TAOLCL)); 171 172 /* Set up initial vectors and matrices */ 173 PetscCall(ParabolicInitialize(&user)); 174 175 PetscCall(Gather(x, user.y, user.state_scatter, user.u, user.design_scatter)); 176 PetscCall(VecDuplicate(x, &x0)); 177 PetscCall(VecCopy(x, x0)); 178 179 /* Set solution vector with an initial guess */ 180 PetscCall(TaoSetSolution(tao, x)); 181 PetscCall(TaoSetObjective(tao, FormFunction, &user)); 182 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user)); 183 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user)); 184 185 PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, &user)); 186 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user)); 187 188 PetscCall(TaoSetFromOptions(tao)); 189 PetscCall(TaoSetStateDesignIS(tao, user.s_is, user.d_is)); 190 191 /* SOLVE THE APPLICATION */ 192 PetscCall(PetscLogStageRegister("Trials", &stages[0])); 193 PetscCall(PetscLogStagePush(stages[0])); 194 user.ksp_its_initial = user.ksp_its; 195 ksp_old = user.ksp_its; 196 for (i = 0; i < ntests; i++) { 197 PetscCall(TaoSolve(tao)); 198 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp_old)); 199 PetscCall(VecCopy(x0, x)); 200 PetscCall(TaoSetSolution(tao, x)); 201 } 202 PetscCall(PetscLogStagePop()); 203 PetscCall(PetscBarrier((PetscObject)x)); 204 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations within initialization: ")); 205 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial)); 206 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Total KSP iterations over %" PetscInt_FMT " trial(s): ", ntests)); 207 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its)); 208 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations per trial: ")); 209 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initial) / ntests)); 210 211 PetscCall(TaoDestroy(&tao)); 212 PetscCall(VecDestroy(&x)); 213 PetscCall(VecDestroy(&x0)); 214 PetscCall(ParabolicDestroy(&user)); 215 216 PetscCall(PetscFinalize()); 217 return 0; 218 } 219 /* ------------------------------------------------------------------- */ 220 /* 221 dwork = Qy - d 222 lwork = L*(u-ur) 223 f = 1/2 * (dwork.dork + alpha*lwork.lwork) 224 */ 225 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr) 226 { 227 PetscReal d1 = 0, d2 = 0; 228 PetscInt i, j; 229 AppCtx *user = (AppCtx *)ptr; 230 231 PetscFunctionBegin; 232 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 233 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 234 for (j = 0; j < user->ns; j++) { 235 i = user->sample_times[j]; 236 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j])); 237 } 238 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns)); 239 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 240 PetscCall(VecDot(user->dwork, user->dwork, &d1)); 241 242 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); 243 PetscCall(MatMult(user->L, user->uwork, user->lwork)); 244 PetscCall(VecDot(user->lwork, user->lwork, &d2)); 245 246 *f = 0.5 * (d1 + user->alpha * d2); 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 /* ------------------------------------------------------------------- */ 251 /* 252 state: g_s = Q' *(Qy - d) 253 design: g_d = alpha*L'*L*(u-ur) 254 */ 255 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr) 256 { 257 PetscInt i, j; 258 AppCtx *user = (AppCtx *)ptr; 259 260 PetscFunctionBegin; 261 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 262 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 263 for (j = 0; j < user->ns; j++) { 264 i = user->sample_times[j]; 265 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j])); 266 } 267 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns)); 268 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 269 PetscCall(Scatter_i(user->dwork, user->di, user->di_scatter, user->ns)); 270 PetscCall(VecSet(user->ywork, 0.0)); 271 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); 272 for (j = 0; j < user->ns; j++) { 273 i = user->sample_times[j]; 274 PetscCall(MatMult(user->QblockT, user->di[j], user->yiwork[i])); 275 } 276 PetscCall(Gather_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); 277 278 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); 279 PetscCall(MatMult(user->L, user->uwork, user->lwork)); 280 PetscCall(MatMult(user->LT, user->lwork, user->uwork)); 281 PetscCall(VecScale(user->uwork, user->alpha)); 282 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 286 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) 287 { 288 PetscReal d1, d2; 289 PetscInt i, j; 290 AppCtx *user = (AppCtx *)ptr; 291 292 PetscFunctionBegin; 293 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 294 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 295 for (j = 0; j < user->ns; j++) { 296 i = user->sample_times[j]; 297 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j])); 298 } 299 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns)); 300 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); 301 PetscCall(VecDot(user->dwork, user->dwork, &d1)); 302 PetscCall(Scatter_i(user->dwork, user->di, user->di_scatter, user->ns)); 303 PetscCall(VecSet(user->ywork, 0.0)); 304 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); 305 for (j = 0; j < user->ns; j++) { 306 i = user->sample_times[j]; 307 PetscCall(MatMult(user->QblockT, user->di[j], user->yiwork[i])); 308 } 309 PetscCall(Gather_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); 310 311 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); 312 PetscCall(MatMult(user->L, user->uwork, user->lwork)); 313 PetscCall(VecDot(user->lwork, user->lwork, &d2)); 314 PetscCall(MatMult(user->LT, user->lwork, user->uwork)); 315 PetscCall(VecScale(user->uwork, user->alpha)); 316 *f = 0.5 * (d1 + user->alpha * d2); 317 318 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 319 PetscFunctionReturn(PETSC_SUCCESS); 320 } 321 322 /* ------------------------------------------------------------------- */ 323 /* A 324 MatShell object 325 */ 326 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr) 327 { 328 AppCtx *user = (AppCtx *)ptr; 329 330 PetscFunctionBegin; 331 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 332 PetscCall(VecSet(user->uwork, 0)); 333 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); 334 PetscCall(VecExp(user->uwork)); 335 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); 336 PetscCall(VecCopy(user->Av_u, user->Swork)); 337 PetscCall(VecReciprocal(user->Swork)); 338 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); 339 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Swork)); 340 if (user->dsg_formed) { 341 PetscCall(MatProductNumeric(user->DSG)); 342 } else { 343 PetscCall(MatMatMult(user->Divwork, user->Grad, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &user->DSG)); 344 user->dsg_formed = PETSC_TRUE; 345 } 346 347 /* B = speye(nx^3) + ht*DSG; */ 348 PetscCall(MatScale(user->DSG, user->ht)); 349 PetscCall(MatShift(user->DSG, 1.0)); 350 PetscFunctionReturn(PETSC_SUCCESS); 351 } 352 353 /* ------------------------------------------------------------------- */ 354 /* B */ 355 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr) 356 { 357 AppCtx *user = (AppCtx *)ptr; 358 359 PetscFunctionBegin; 360 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 361 PetscFunctionReturn(PETSC_SUCCESS); 362 } 363 364 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y) 365 { 366 PetscInt i; 367 AppCtx *user; 368 369 PetscFunctionBegin; 370 PetscCall(MatShellGetContext(J_shell, &user)); 371 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); 372 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 373 for (i = 1; i < user->nt; i++) { 374 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 375 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1])); 376 } 377 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 378 PetscFunctionReturn(PETSC_SUCCESS); 379 } 380 381 PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y) 382 { 383 PetscInt i; 384 AppCtx *user; 385 386 PetscFunctionBegin; 387 PetscCall(MatShellGetContext(J_shell, &user)); 388 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); 389 for (i = 0; i < user->nt - 1; i++) { 390 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 391 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i + 1])); 392 } 393 i = user->nt - 1; 394 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 395 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 396 PetscFunctionReturn(PETSC_SUCCESS); 397 } 398 399 PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y) 400 { 401 AppCtx *user; 402 403 PetscFunctionBegin; 404 PetscCall(MatShellGetContext(J_shell, &user)); 405 PetscCall(MatMult(user->Grad, X, user->Swork)); 406 PetscCall(VecPointwiseDivide(user->Swork, user->Swork, user->Av_u)); 407 PetscCall(MatMult(user->Div, user->Swork, Y)); 408 PetscCall(VecAYPX(Y, user->ht, X)); 409 PetscFunctionReturn(PETSC_SUCCESS); 410 } 411 412 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y) 413 { 414 PetscInt i; 415 AppCtx *user; 416 417 PetscFunctionBegin; 418 PetscCall(MatShellGetContext(J_shell, &user)); 419 420 /* sdiag(1./v) */ 421 PetscCall(VecSet(user->uwork, 0)); 422 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); 423 PetscCall(VecExp(user->uwork)); 424 425 /* sdiag(1./((Av*(1./v)).^2)) */ 426 PetscCall(MatMult(user->Av, user->uwork, user->Swork)); 427 PetscCall(VecPointwiseMult(user->Swork, user->Swork, user->Swork)); 428 PetscCall(VecReciprocal(user->Swork)); 429 430 /* (Av * (sdiag(1./v) * b)) */ 431 PetscCall(VecPointwiseMult(user->uwork, user->uwork, X)); 432 PetscCall(MatMult(user->Av, user->uwork, user->Twork)); 433 434 /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */ 435 PetscCall(VecPointwiseMult(user->Swork, user->Twork, user->Swork)); 436 437 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 438 for (i = 0; i < user->nt; i++) { 439 /* (sdiag(Grad*y(:,i)) */ 440 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork)); 441 442 /* ht * Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */ 443 PetscCall(VecPointwiseMult(user->Twork, user->Twork, user->Swork)); 444 PetscCall(MatMult(user->Div, user->Twork, user->yiwork[i])); 445 PetscCall(VecScale(user->yiwork[i], user->ht)); 446 } 447 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 448 449 PetscFunctionReturn(PETSC_SUCCESS); 450 } 451 452 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 453 { 454 PetscInt i; 455 AppCtx *user; 456 457 PetscFunctionBegin; 458 PetscCall(MatShellGetContext(J_shell, &user)); 459 460 /* sdiag(1./((Av*(1./v)).^2)) */ 461 PetscCall(VecSet(user->uwork, 0)); 462 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); 463 PetscCall(VecExp(user->uwork)); 464 PetscCall(MatMult(user->Av, user->uwork, user->Rwork)); 465 PetscCall(VecPointwiseMult(user->Rwork, user->Rwork, user->Rwork)); 466 PetscCall(VecReciprocal(user->Rwork)); 467 468 PetscCall(VecSet(Y, 0.0)); 469 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 470 PetscCall(Scatter_i(X, user->yiwork, user->yi_scatter, user->nt)); 471 for (i = 0; i < user->nt; i++) { 472 /* (Div' * b(:,i)) */ 473 PetscCall(MatMult(user->Grad, user->yiwork[i], user->Swork)); 474 475 /* sdiag(Grad*y(:,i)) */ 476 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork)); 477 478 /* (sdiag(Grad*y(:,i)) * (Div' * b(:,i))) */ 479 PetscCall(VecPointwiseMult(user->Twork, user->Swork, user->Twork)); 480 481 /* (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i)))) */ 482 PetscCall(VecPointwiseMult(user->Twork, user->Rwork, user->Twork)); 483 484 /* (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 485 PetscCall(MatMult(user->AvT, user->Twork, user->yiwork[i])); 486 487 /* sdiag(1./v) * (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 488 PetscCall(VecPointwiseMult(user->yiwork[i], user->uwork, user->yiwork[i])); 489 PetscCall(VecAXPY(Y, user->ht, user->yiwork[i])); 490 } 491 PetscFunctionReturn(PETSC_SUCCESS); 492 } 493 494 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y) 495 { 496 AppCtx *user; 497 498 PetscFunctionBegin; 499 PetscCall(PCShellGetContext(PC_shell, &user)); 500 501 if (user->dsg_formed) { 502 PetscCall(MatSOR(user->DSG, X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y)); 503 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "DSG not formed"); 504 PetscFunctionReturn(PETSC_SUCCESS); 505 } 506 507 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y) 508 { 509 AppCtx *user; 510 PetscInt its, i; 511 512 PetscFunctionBegin; 513 PetscCall(MatShellGetContext(J_shell, &user)); 514 515 if (Y == user->ytrue) { 516 /* First solve is done with true solution to set up problem */ 517 PetscCall(KSPSetTolerances(user->solver, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 518 } else { 519 PetscCall(KSPSetTolerances(user->solver, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 520 } 521 522 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); 523 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0])); 524 PetscCall(KSPGetIterationNumber(user->solver, &its)); 525 user->ksp_its = user->ksp_its + its; 526 527 for (i = 1; i < user->nt; i++) { 528 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i - 1])); 529 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 530 PetscCall(KSPGetIterationNumber(user->solver, &its)); 531 user->ksp_its = user->ksp_its + its; 532 } 533 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 534 PetscFunctionReturn(PETSC_SUCCESS); 535 } 536 537 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y) 538 { 539 AppCtx *user; 540 PetscInt its, i; 541 542 PetscFunctionBegin; 543 PetscCall(MatShellGetContext(J_shell, &user)); 544 545 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); 546 547 i = user->nt - 1; 548 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 549 550 PetscCall(KSPGetIterationNumber(user->solver, &its)); 551 user->ksp_its = user->ksp_its + its; 552 553 for (i = user->nt - 2; i >= 0; i--) { 554 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i + 1])); 555 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 556 557 PetscCall(KSPGetIterationNumber(user->solver, &its)); 558 user->ksp_its = user->ksp_its + its; 559 } 560 561 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 562 PetscFunctionReturn(PETSC_SUCCESS); 563 } 564 565 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell) 566 { 567 AppCtx *user; 568 569 PetscFunctionBegin; 570 PetscCall(MatShellGetContext(J_shell, &user)); 571 572 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell)); 573 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (void (*)(void))StateMatMult)); 574 PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate)); 575 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose)); 576 PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal)); 577 PetscFunctionReturn(PETSC_SUCCESS); 578 } 579 580 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X) 581 { 582 AppCtx *user; 583 584 PetscFunctionBegin; 585 PetscCall(MatShellGetContext(J_shell, &user)); 586 PetscCall(VecCopy(user->js_diag, X)); 587 PetscFunctionReturn(PETSC_SUCCESS); 588 } 589 590 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr) 591 { 592 /* con = Ay - q, A = [B 0 0 ... 0; 593 -I B 0 ... 0; 594 0 -I B ... 0; 595 ... ; 596 0 ... -I B] 597 B = ht * Div * Sigma * Grad + eye */ 598 PetscInt i; 599 AppCtx *user = (AppCtx *)ptr; 600 601 PetscFunctionBegin; 602 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 603 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 604 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 605 for (i = 1; i < user->nt; i++) { 606 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 607 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1])); 608 } 609 PetscCall(Gather_i(C, user->yiwork, user->yi_scatter, user->nt)); 610 PetscCall(VecAXPY(C, -1.0, user->q)); 611 PetscFunctionReturn(PETSC_SUCCESS); 612 } 613 614 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 615 { 616 PetscFunctionBegin; 617 PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 618 PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 619 PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 620 PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 621 PetscFunctionReturn(PETSC_SUCCESS); 622 } 623 624 PetscErrorCode Scatter_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 625 { 626 PetscInt i; 627 628 PetscFunctionBegin; 629 for (i = 0; i < nt; i++) { 630 PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 631 PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 632 } 633 PetscFunctionReturn(PETSC_SUCCESS); 634 } 635 636 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 637 { 638 PetscFunctionBegin; 639 PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 640 PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 641 PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 642 PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 643 PetscFunctionReturn(PETSC_SUCCESS); 644 } 645 646 PetscErrorCode Gather_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 647 { 648 PetscInt i; 649 650 PetscFunctionBegin; 651 for (i = 0; i < nt; i++) { 652 PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 653 PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 654 } 655 PetscFunctionReturn(PETSC_SUCCESS); 656 } 657 658 PetscErrorCode ParabolicInitialize(AppCtx *user) 659 { 660 PetscInt m, n, i, j, k, linear_index, istart, iend, iblock, lo, hi, lo2, hi2; 661 Vec XX, YY, ZZ, XXwork, YYwork, ZZwork, UTwork, yi, di, bc; 662 PetscReal *x, *y, *z; 663 PetscReal h, stime; 664 PetscScalar hinv, neg_hinv, half = 0.5, sqrt_beta; 665 PetscInt im, indx1, indx2, indy1, indy2, indz1, indz2, nx, ny, nz; 666 PetscReal xri, yri, zri, xim, yim, zim, dx1, dx2, dy1, dy2, dz1, dz2, Dx, Dy, Dz; 667 PetscScalar v, vx, vy, vz; 668 IS is_from_y, is_to_yi, is_from_d, is_to_di; 669 /* Data locations */ 670 PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160, 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774, 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997, 671 0.5198, 0.8326, 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728, 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273, 0.5724, 0.4331, 0.5136, 0.3547, 672 0.4413, 0.2602, 0.5698, 0.7278, 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594, 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756}; 673 674 PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256, 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792, 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437, 675 0.5938, 0.6137, 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985, 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288, 0.5761, 0.1129, 0.3841, 0.6150, 676 0.6904, 0.6686, 0.1361, 0.4601, 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480, 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658}; 677 678 PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295, 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819, 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236, 679 0.8800, 0.2939, 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712, 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078, 0.1987, 0.2297, 0.4321, 0.8115, 680 0.4915, 0.7764, 0.4657, 0.4627, 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313, 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711}; 681 682 PetscFunctionBegin; 683 PetscCall(PetscMalloc1(user->mx, &x)); 684 PetscCall(PetscMalloc1(user->mx, &y)); 685 PetscCall(PetscMalloc1(user->mx, &z)); 686 user->jformed = PETSC_FALSE; 687 user->dsg_formed = PETSC_FALSE; 688 689 n = user->mx * user->mx * user->mx; 690 m = 3 * user->mx * user->mx * (user->mx - 1); 691 sqrt_beta = PetscSqrtScalar(user->beta); 692 693 user->ksp_its = 0; 694 user->ksp_its_initial = 0; 695 696 stime = (PetscReal)user->nt / user->ns; 697 PetscCall(PetscMalloc1(user->ns, &user->sample_times)); 698 for (i = 0; i < user->ns; i++) user->sample_times[i] = (PetscInt)(stime * ((PetscReal)i + 1.0) - 0.5); 699 700 PetscCall(VecCreate(PETSC_COMM_WORLD, &XX)); 701 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q)); 702 PetscCall(VecSetSizes(XX, PETSC_DECIDE, n)); 703 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt)); 704 PetscCall(VecSetFromOptions(XX)); 705 PetscCall(VecSetFromOptions(user->q)); 706 707 PetscCall(VecDuplicate(XX, &YY)); 708 PetscCall(VecDuplicate(XX, &ZZ)); 709 PetscCall(VecDuplicate(XX, &XXwork)); 710 PetscCall(VecDuplicate(XX, &YYwork)); 711 PetscCall(VecDuplicate(XX, &ZZwork)); 712 PetscCall(VecDuplicate(XX, &UTwork)); 713 PetscCall(VecDuplicate(XX, &user->utrue)); 714 PetscCall(VecDuplicate(XX, &bc)); 715 716 /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */ 717 h = 1.0 / user->mx; 718 hinv = user->mx; 719 neg_hinv = -hinv; 720 721 PetscCall(VecGetOwnershipRange(XX, &istart, &iend)); 722 for (linear_index = istart; linear_index < iend; linear_index++) { 723 i = linear_index % user->mx; 724 j = ((linear_index - i) / user->mx) % user->mx; 725 k = ((linear_index - i) / user->mx - j) / user->mx; 726 vx = h * (i + 0.5); 727 vy = h * (j + 0.5); 728 vz = h * (k + 0.5); 729 PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES)); 730 PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES)); 731 PetscCall(VecSetValues(ZZ, 1, &linear_index, &vz, INSERT_VALUES)); 732 if ((vx < 0.6) && (vx > 0.4) && (vy < 0.6) && (vy > 0.4) && (vy < 0.6) && (vz < 0.6) && (vz > 0.4)) { 733 v = 1000.0; 734 PetscCall(VecSetValues(bc, 1, &linear_index, &v, INSERT_VALUES)); 735 } 736 } 737 738 PetscCall(VecAssemblyBegin(XX)); 739 PetscCall(VecAssemblyEnd(XX)); 740 PetscCall(VecAssemblyBegin(YY)); 741 PetscCall(VecAssemblyEnd(YY)); 742 PetscCall(VecAssemblyBegin(ZZ)); 743 PetscCall(VecAssemblyEnd(ZZ)); 744 PetscCall(VecAssemblyBegin(bc)); 745 PetscCall(VecAssemblyEnd(bc)); 746 747 /* Compute true parameter function 748 ut = 0.5 + exp(-10*((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)) */ 749 PetscCall(VecCopy(XX, XXwork)); 750 PetscCall(VecCopy(YY, YYwork)); 751 PetscCall(VecCopy(ZZ, ZZwork)); 752 753 PetscCall(VecShift(XXwork, -0.5)); 754 PetscCall(VecShift(YYwork, -0.5)); 755 PetscCall(VecShift(ZZwork, -0.5)); 756 757 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 758 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 759 PetscCall(VecPointwiseMult(ZZwork, ZZwork, ZZwork)); 760 761 PetscCall(VecCopy(XXwork, user->utrue)); 762 PetscCall(VecAXPY(user->utrue, 1.0, YYwork)); 763 PetscCall(VecAXPY(user->utrue, 1.0, ZZwork)); 764 PetscCall(VecScale(user->utrue, -10.0)); 765 PetscCall(VecExp(user->utrue)); 766 PetscCall(VecShift(user->utrue, 0.5)); 767 768 PetscCall(VecDestroy(&XX)); 769 PetscCall(VecDestroy(&YY)); 770 PetscCall(VecDestroy(&ZZ)); 771 PetscCall(VecDestroy(&XXwork)); 772 PetscCall(VecDestroy(&YYwork)); 773 PetscCall(VecDestroy(&ZZwork)); 774 PetscCall(VecDestroy(&UTwork)); 775 776 /* Initial guess and reference model */ 777 PetscCall(VecDuplicate(user->utrue, &user->ur)); 778 PetscCall(VecSet(user->ur, 0.0)); 779 780 /* Generate Grad matrix */ 781 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad)); 782 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, m, n)); 783 PetscCall(MatSetFromOptions(user->Grad)); 784 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 2, NULL, 2, NULL)); 785 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 2, NULL)); 786 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend)); 787 788 for (i = istart; i < iend; i++) { 789 if (i < m / 3) { 790 iblock = i / (user->mx - 1); 791 j = iblock * user->mx + (i % (user->mx - 1)); 792 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 793 j = j + 1; 794 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 795 } 796 if (i >= m / 3 && i < 2 * m / 3) { 797 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); 798 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); 799 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 800 j = j + user->mx; 801 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 802 } 803 if (i >= 2 * m / 3) { 804 j = i - 2 * m / 3; 805 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 806 j = j + user->mx * user->mx; 807 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 808 } 809 } 810 811 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY)); 812 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY)); 813 814 /* Generate arithmetic averaging matrix Av */ 815 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Av)); 816 PetscCall(MatSetSizes(user->Av, PETSC_DECIDE, PETSC_DECIDE, m, n)); 817 PetscCall(MatSetFromOptions(user->Av)); 818 PetscCall(MatMPIAIJSetPreallocation(user->Av, 2, NULL, 2, NULL)); 819 PetscCall(MatSeqAIJSetPreallocation(user->Av, 2, NULL)); 820 PetscCall(MatGetOwnershipRange(user->Av, &istart, &iend)); 821 822 for (i = istart; i < iend; i++) { 823 if (i < m / 3) { 824 iblock = i / (user->mx - 1); 825 j = iblock * user->mx + (i % (user->mx - 1)); 826 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 827 j = j + 1; 828 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 829 } 830 if (i >= m / 3 && i < 2 * m / 3) { 831 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); 832 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); 833 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 834 j = j + user->mx; 835 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 836 } 837 if (i >= 2 * m / 3) { 838 j = i - 2 * m / 3; 839 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 840 j = j + user->mx * user->mx; 841 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 842 } 843 } 844 845 PetscCall(MatAssemblyBegin(user->Av, MAT_FINAL_ASSEMBLY)); 846 PetscCall(MatAssemblyEnd(user->Av, MAT_FINAL_ASSEMBLY)); 847 848 /* Generate transpose of averaging matrix Av */ 849 PetscCall(MatTranspose(user->Av, MAT_INITIAL_MATRIX, &user->AvT)); 850 851 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->L)); 852 PetscCall(MatSetSizes(user->L, PETSC_DECIDE, PETSC_DECIDE, m + n, n)); 853 PetscCall(MatSetFromOptions(user->L)); 854 PetscCall(MatMPIAIJSetPreallocation(user->L, 2, NULL, 2, NULL)); 855 PetscCall(MatSeqAIJSetPreallocation(user->L, 2, NULL)); 856 PetscCall(MatGetOwnershipRange(user->L, &istart, &iend)); 857 858 for (i = istart; i < iend; i++) { 859 if (i < m / 3) { 860 iblock = i / (user->mx - 1); 861 j = iblock * user->mx + (i % (user->mx - 1)); 862 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 863 j = j + 1; 864 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 865 } 866 if (i >= m / 3 && i < 2 * m / 3) { 867 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); 868 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); 869 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 870 j = j + user->mx; 871 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 872 } 873 if (i >= 2 * m / 3 && i < m) { 874 j = i - 2 * m / 3; 875 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 876 j = j + user->mx * user->mx; 877 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 878 } 879 if (i >= m) { 880 j = i - m; 881 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &sqrt_beta, INSERT_VALUES)); 882 } 883 } 884 885 PetscCall(MatAssemblyBegin(user->L, MAT_FINAL_ASSEMBLY)); 886 PetscCall(MatAssemblyEnd(user->L, MAT_FINAL_ASSEMBLY)); 887 888 PetscCall(MatScale(user->L, PetscPowScalar(h, 1.5))); 889 890 /* Generate Div matrix */ 891 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div)); 892 893 /* Build work vectors and matrices */ 894 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->S)); 895 PetscCall(VecSetSizes(user->S, PETSC_DECIDE, user->mx * user->mx * (user->mx - 1) * 3)); 896 PetscCall(VecSetFromOptions(user->S)); 897 898 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork)); 899 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, m + user->mx * user->mx * user->mx)); 900 PetscCall(VecSetFromOptions(user->lwork)); 901 902 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork)); 903 PetscCall(MatDuplicate(user->Av, MAT_SHARE_NONZERO_PATTERN, &user->Avwork)); 904 905 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->d)); 906 PetscCall(VecSetSizes(user->d, PETSC_DECIDE, user->ndata * user->nt)); 907 PetscCall(VecSetFromOptions(user->d)); 908 909 PetscCall(VecDuplicate(user->S, &user->Swork)); 910 PetscCall(VecDuplicate(user->S, &user->Av_u)); 911 PetscCall(VecDuplicate(user->S, &user->Twork)); 912 PetscCall(VecDuplicate(user->S, &user->Rwork)); 913 PetscCall(VecDuplicate(user->y, &user->ywork)); 914 PetscCall(VecDuplicate(user->u, &user->uwork)); 915 PetscCall(VecDuplicate(user->u, &user->js_diag)); 916 PetscCall(VecDuplicate(user->c, &user->cwork)); 917 PetscCall(VecDuplicate(user->d, &user->dwork)); 918 919 /* Create matrix-free shell user->Js for computing A*x */ 920 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &user->Js)); 921 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (void (*)(void))StateMatMult)); 922 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate)); 923 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose)); 924 PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal)); 925 926 /* Diagonal blocks of user->Js */ 927 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlock)); 928 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (void (*)(void))StateMatBlockMult)); 929 /* Blocks are symmetric */ 930 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockMult)); 931 932 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U, 933 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular. 934 This is an SSOR preconditioner for user->JsBlock. */ 935 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlockPrec)); 936 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (void (*)(void))StateMatBlockPrecMult)); 937 /* JsBlockPrec is symmetric */ 938 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockPrecMult)); 939 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRIC, PETSC_TRUE)); 940 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRY_ETERNAL, PETSC_TRUE)); 941 942 /* Create a matrix-free shell user->Jd for computing B*x */ 943 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m, user, &user->Jd)); 944 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (void (*)(void))DesignMatMult)); 945 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (void (*)(void))DesignMatMultTranspose)); 946 947 /* User-defined routines for computing user->Js\x and user->Js^T\x*/ 948 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &user->JsInv)); 949 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (void (*)(void))StateMatInvMult)); 950 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatInvTransposeMult)); 951 952 /* Solver options and tolerances */ 953 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver)); 954 PetscCall(KSPSetType(user->solver, KSPCG)); 955 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec)); 956 PetscCall(KSPSetInitialGuessNonzero(user->solver, PETSC_FALSE)); 957 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500)); 958 PetscCall(KSPSetFromOptions(user->solver)); 959 PetscCall(KSPGetPC(user->solver, &user->prec)); 960 PetscCall(PCSetType(user->prec, PCSHELL)); 961 962 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult)); 963 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMult)); 964 PetscCall(PCShellSetContext(user->prec, user)); 965 966 /* Create scatter from y to y_1,y_2,...,y_nt */ 967 PetscCall(PetscMalloc1(user->nt * user->m, &user->yi_scatter)); 968 PetscCall(VecCreate(PETSC_COMM_WORLD, &yi)); 969 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx * user->mx)); 970 PetscCall(VecSetFromOptions(yi)); 971 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi)); 972 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork)); 973 974 PetscCall(VecGetOwnershipRange(user->y, &lo2, &hi2)); 975 istart = 0; 976 for (i = 0; i < user->nt; i++) { 977 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi)); 978 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi)); 979 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo2 + istart, 1, &is_from_y)); 980 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i])); 981 istart = istart + hi - lo; 982 PetscCall(ISDestroy(&is_to_yi)); 983 PetscCall(ISDestroy(&is_from_y)); 984 } 985 PetscCall(VecDestroy(&yi)); 986 987 /* Create scatter from d to d_1,d_2,...,d_ns */ 988 PetscCall(PetscMalloc1(user->ns * user->ndata, &user->di_scatter)); 989 PetscCall(VecCreate(PETSC_COMM_WORLD, &di)); 990 PetscCall(VecSetSizes(di, PETSC_DECIDE, user->ndata)); 991 PetscCall(VecSetFromOptions(di)); 992 PetscCall(VecDuplicateVecs(di, user->ns, &user->di)); 993 PetscCall(VecGetOwnershipRange(user->d, &lo2, &hi2)); 994 istart = 0; 995 for (i = 0; i < user->ns; i++) { 996 PetscCall(VecGetOwnershipRange(user->di[i], &lo, &hi)); 997 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_di)); 998 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo2 + istart, 1, &is_from_d)); 999 PetscCall(VecScatterCreate(user->d, is_from_d, user->di[i], is_to_di, &user->di_scatter[i])); 1000 istart = istart + hi - lo; 1001 PetscCall(ISDestroy(&is_to_di)); 1002 PetscCall(ISDestroy(&is_from_d)); 1003 } 1004 PetscCall(VecDestroy(&di)); 1005 1006 /* Assemble RHS of forward problem */ 1007 PetscCall(VecCopy(bc, user->yiwork[0])); 1008 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0)); 1009 PetscCall(Gather_i(user->q, user->yiwork, user->yi_scatter, user->nt)); 1010 PetscCall(VecDestroy(&bc)); 1011 1012 /* Compute true state function yt given ut */ 1013 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue)); 1014 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt)); 1015 PetscCall(VecSetFromOptions(user->ytrue)); 1016 1017 /* First compute Av_u = Av*exp(-u) */ 1018 PetscCall(VecSet(user->uwork, 0)); 1019 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); /* Note: user->utrue */ 1020 PetscCall(VecExp(user->uwork)); 1021 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); 1022 1023 /* Symbolic DSG = Div * Grad */ 1024 PetscCall(MatProductCreate(user->Div, user->Grad, NULL, &user->DSG)); 1025 PetscCall(MatProductSetType(user->DSG, MATPRODUCT_AB)); 1026 PetscCall(MatProductSetAlgorithm(user->DSG, "default")); 1027 PetscCall(MatProductSetFill(user->DSG, PETSC_DEFAULT)); 1028 PetscCall(MatProductSetFromOptions(user->DSG)); 1029 PetscCall(MatProductSymbolic(user->DSG)); 1030 1031 user->dsg_formed = PETSC_TRUE; 1032 1033 /* Next form DSG = Div*Grad */ 1034 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); 1035 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u)); 1036 if (user->dsg_formed) { 1037 PetscCall(MatProductNumeric(user->DSG)); 1038 } else { 1039 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &user->DSG)); 1040 user->dsg_formed = PETSC_TRUE; 1041 } 1042 /* B = speye(nx^3) + ht*DSG; */ 1043 PetscCall(MatScale(user->DSG, user->ht)); 1044 PetscCall(MatShift(user->DSG, 1.0)); 1045 1046 /* Now solve for ytrue */ 1047 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue)); 1048 1049 /* Initial guess y0 for state given u0 */ 1050 1051 /* First compute Av_u = Av*exp(-u) */ 1052 PetscCall(VecSet(user->uwork, 0)); 1053 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); /* Note: user->u */ 1054 PetscCall(VecExp(user->uwork)); 1055 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); 1056 1057 /* Next form DSG = Div*S*Grad */ 1058 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); 1059 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u)); 1060 if (user->dsg_formed) { 1061 PetscCall(MatProductNumeric(user->DSG)); 1062 } else { 1063 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &user->DSG)); 1064 1065 user->dsg_formed = PETSC_TRUE; 1066 } 1067 /* B = speye(nx^3) + ht*DSG; */ 1068 PetscCall(MatScale(user->DSG, user->ht)); 1069 PetscCall(MatShift(user->DSG, 1.0)); 1070 1071 /* Now solve for y */ 1072 PetscCall(StateMatInvMult(user->Js, user->q, user->y)); 1073 1074 /* Construct projection matrix Q, a block diagonal matrix consisting of nt copies of Qblock along the diagonal */ 1075 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Qblock)); 1076 PetscCall(MatSetSizes(user->Qblock, PETSC_DECIDE, PETSC_DECIDE, user->ndata, n)); 1077 PetscCall(MatSetFromOptions(user->Qblock)); 1078 PetscCall(MatMPIAIJSetPreallocation(user->Qblock, 8, NULL, 8, NULL)); 1079 PetscCall(MatSeqAIJSetPreallocation(user->Qblock, 8, NULL)); 1080 1081 for (i = 0; i < user->mx; i++) { 1082 x[i] = h * (i + 0.5); 1083 y[i] = h * (i + 0.5); 1084 z[i] = h * (i + 0.5); 1085 } 1086 1087 PetscCall(MatGetOwnershipRange(user->Qblock, &istart, &iend)); 1088 nx = user->mx; 1089 ny = user->mx; 1090 nz = user->mx; 1091 for (i = istart; i < iend; i++) { 1092 xri = xr[i]; 1093 im = 0; 1094 xim = x[im]; 1095 while (xri > xim && im < nx) { 1096 im = im + 1; 1097 xim = x[im]; 1098 } 1099 indx1 = im - 1; 1100 indx2 = im; 1101 dx1 = xri - x[indx1]; 1102 dx2 = x[indx2] - xri; 1103 1104 yri = yr[i]; 1105 im = 0; 1106 yim = y[im]; 1107 while (yri > yim && im < ny) { 1108 im = im + 1; 1109 yim = y[im]; 1110 } 1111 indy1 = im - 1; 1112 indy2 = im; 1113 dy1 = yri - y[indy1]; 1114 dy2 = y[indy2] - yri; 1115 1116 zri = zr[i]; 1117 im = 0; 1118 zim = z[im]; 1119 while (zri > zim && im < nz) { 1120 im = im + 1; 1121 zim = z[im]; 1122 } 1123 indz1 = im - 1; 1124 indz2 = im; 1125 dz1 = zri - z[indz1]; 1126 dz2 = z[indz2] - zri; 1127 1128 Dx = x[indx2] - x[indx1]; 1129 Dy = y[indy2] - y[indy1]; 1130 Dz = z[indz2] - z[indz1]; 1131 1132 j = indx1 + indy1 * nx + indz1 * nx * ny; 1133 v = (1 - dx1 / Dx) * (1 - dy1 / Dy) * (1 - dz1 / Dz); 1134 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1135 1136 j = indx1 + indy1 * nx + indz2 * nx * ny; 1137 v = (1 - dx1 / Dx) * (1 - dy1 / Dy) * (1 - dz2 / Dz); 1138 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1139 1140 j = indx1 + indy2 * nx + indz1 * nx * ny; 1141 v = (1 - dx1 / Dx) * (1 - dy2 / Dy) * (1 - dz1 / Dz); 1142 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1143 1144 j = indx1 + indy2 * nx + indz2 * nx * ny; 1145 v = (1 - dx1 / Dx) * (1 - dy2 / Dy) * (1 - dz2 / Dz); 1146 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1147 1148 j = indx2 + indy1 * nx + indz1 * nx * ny; 1149 v = (1 - dx2 / Dx) * (1 - dy1 / Dy) * (1 - dz1 / Dz); 1150 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1151 1152 j = indx2 + indy1 * nx + indz2 * nx * ny; 1153 v = (1 - dx2 / Dx) * (1 - dy1 / Dy) * (1 - dz2 / Dz); 1154 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1155 1156 j = indx2 + indy2 * nx + indz1 * nx * ny; 1157 v = (1 - dx2 / Dx) * (1 - dy2 / Dy) * (1 - dz1 / Dz); 1158 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1159 1160 j = indx2 + indy2 * nx + indz2 * nx * ny; 1161 v = (1 - dx2 / Dx) * (1 - dy2 / Dy) * (1 - dz2 / Dz); 1162 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1163 } 1164 PetscCall(MatAssemblyBegin(user->Qblock, MAT_FINAL_ASSEMBLY)); 1165 PetscCall(MatAssemblyEnd(user->Qblock, MAT_FINAL_ASSEMBLY)); 1166 1167 PetscCall(MatTranspose(user->Qblock, MAT_INITIAL_MATRIX, &user->QblockT)); 1168 PetscCall(MatTranspose(user->L, MAT_INITIAL_MATRIX, &user->LT)); 1169 1170 /* Add noise to the measurement data */ 1171 PetscCall(VecSet(user->ywork, 1.0)); 1172 PetscCall(VecAYPX(user->ywork, user->noise, user->ytrue)); 1173 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); 1174 for (j = 0; j < user->ns; j++) { 1175 i = user->sample_times[j]; 1176 PetscCall(MatMult(user->Qblock, user->yiwork[i], user->di[j])); 1177 } 1178 PetscCall(Gather_i(user->d, user->di, user->di_scatter, user->ns)); 1179 1180 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */ 1181 PetscCall(KSPSetFromOptions(user->solver)); 1182 PetscCall(PetscFree(x)); 1183 PetscCall(PetscFree(y)); 1184 PetscCall(PetscFree(z)); 1185 PetscFunctionReturn(PETSC_SUCCESS); 1186 } 1187 1188 PetscErrorCode ParabolicDestroy(AppCtx *user) 1189 { 1190 PetscInt i; 1191 1192 PetscFunctionBegin; 1193 PetscCall(MatDestroy(&user->Qblock)); 1194 PetscCall(MatDestroy(&user->QblockT)); 1195 PetscCall(MatDestroy(&user->Div)); 1196 PetscCall(MatDestroy(&user->Divwork)); 1197 PetscCall(MatDestroy(&user->Grad)); 1198 PetscCall(MatDestroy(&user->Av)); 1199 PetscCall(MatDestroy(&user->Avwork)); 1200 PetscCall(MatDestroy(&user->AvT)); 1201 PetscCall(MatDestroy(&user->DSG)); 1202 PetscCall(MatDestroy(&user->L)); 1203 PetscCall(MatDestroy(&user->LT)); 1204 PetscCall(KSPDestroy(&user->solver)); 1205 PetscCall(MatDestroy(&user->Js)); 1206 PetscCall(MatDestroy(&user->Jd)); 1207 PetscCall(MatDestroy(&user->JsInv)); 1208 PetscCall(MatDestroy(&user->JsBlock)); 1209 PetscCall(MatDestroy(&user->JsBlockPrec)); 1210 PetscCall(VecDestroy(&user->u)); 1211 PetscCall(VecDestroy(&user->uwork)); 1212 PetscCall(VecDestroy(&user->utrue)); 1213 PetscCall(VecDestroy(&user->y)); 1214 PetscCall(VecDestroy(&user->ywork)); 1215 PetscCall(VecDestroy(&user->ytrue)); 1216 PetscCall(VecDestroyVecs(user->nt, &user->yi)); 1217 PetscCall(VecDestroyVecs(user->nt, &user->yiwork)); 1218 PetscCall(VecDestroyVecs(user->ns, &user->di)); 1219 PetscCall(PetscFree(user->yi)); 1220 PetscCall(PetscFree(user->yiwork)); 1221 PetscCall(PetscFree(user->di)); 1222 PetscCall(VecDestroy(&user->c)); 1223 PetscCall(VecDestroy(&user->cwork)); 1224 PetscCall(VecDestroy(&user->ur)); 1225 PetscCall(VecDestroy(&user->q)); 1226 PetscCall(VecDestroy(&user->d)); 1227 PetscCall(VecDestroy(&user->dwork)); 1228 PetscCall(VecDestroy(&user->lwork)); 1229 PetscCall(VecDestroy(&user->S)); 1230 PetscCall(VecDestroy(&user->Swork)); 1231 PetscCall(VecDestroy(&user->Av_u)); 1232 PetscCall(VecDestroy(&user->Twork)); 1233 PetscCall(VecDestroy(&user->Rwork)); 1234 PetscCall(VecDestroy(&user->js_diag)); 1235 PetscCall(ISDestroy(&user->s_is)); 1236 PetscCall(ISDestroy(&user->d_is)); 1237 PetscCall(VecScatterDestroy(&user->state_scatter)); 1238 PetscCall(VecScatterDestroy(&user->design_scatter)); 1239 for (i = 0; i < user->nt; i++) PetscCall(VecScatterDestroy(&user->yi_scatter[i])); 1240 for (i = 0; i < user->ns; i++) PetscCall(VecScatterDestroy(&user->di_scatter[i])); 1241 PetscCall(PetscFree(user->yi_scatter)); 1242 PetscCall(PetscFree(user->di_scatter)); 1243 PetscCall(PetscFree(user->sample_times)); 1244 PetscFunctionReturn(PETSC_SUCCESS); 1245 } 1246 1247 PetscErrorCode ParabolicMonitor(Tao tao, void *ptr) 1248 { 1249 Vec X; 1250 PetscReal unorm, ynorm; 1251 AppCtx *user = (AppCtx *)ptr; 1252 1253 PetscFunctionBegin; 1254 PetscCall(TaoGetSolution(tao, &X)); 1255 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 1256 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue)); 1257 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); 1258 PetscCall(VecNorm(user->uwork, NORM_2, &unorm)); 1259 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm)); 1260 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm)); 1261 PetscFunctionReturn(PETSC_SUCCESS); 1262 } 1263 1264 /*TEST 1265 1266 build: 1267 requires: !complex 1268 1269 test: 1270 args: -tao_cmonitor -tao_type lcl -ns 1 -tao_gatol 1.e-4 -ksp_max_it 30 1271 requires: !single 1272 1273 TEST*/ 1274