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 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 452 { 453 PetscInt i; 454 AppCtx *user; 455 456 PetscFunctionBegin; 457 PetscCall(MatShellGetContext(J_shell, &user)); 458 459 /* sdiag(1./((Av*(1./v)).^2)) */ 460 PetscCall(VecSet(user->uwork, 0)); 461 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); 462 PetscCall(VecExp(user->uwork)); 463 PetscCall(MatMult(user->Av, user->uwork, user->Rwork)); 464 PetscCall(VecPointwiseMult(user->Rwork, user->Rwork, user->Rwork)); 465 PetscCall(VecReciprocal(user->Rwork)); 466 467 PetscCall(VecSet(Y, 0.0)); 468 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 469 PetscCall(Scatter_i(X, user->yiwork, user->yi_scatter, user->nt)); 470 for (i = 0; i < user->nt; i++) { 471 /* (Div' * b(:,i)) */ 472 PetscCall(MatMult(user->Grad, user->yiwork[i], user->Swork)); 473 474 /* sdiag(Grad*y(:,i)) */ 475 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork)); 476 477 /* (sdiag(Grad*y(:,i)) * (Div' * b(:,i))) */ 478 PetscCall(VecPointwiseMult(user->Twork, user->Swork, user->Twork)); 479 480 /* (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i)))) */ 481 PetscCall(VecPointwiseMult(user->Twork, user->Rwork, user->Twork)); 482 483 /* (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 484 PetscCall(MatMult(user->AvT, user->Twork, user->yiwork[i])); 485 486 /* sdiag(1./v) * (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 487 PetscCall(VecPointwiseMult(user->yiwork[i], user->uwork, user->yiwork[i])); 488 PetscCall(VecAXPY(Y, user->ht, user->yiwork[i])); 489 } 490 PetscFunctionReturn(PETSC_SUCCESS); 491 } 492 493 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y) 494 { 495 AppCtx *user; 496 497 PetscFunctionBegin; 498 PetscCall(PCShellGetContext(PC_shell, &user)); 499 500 if (user->dsg_formed) { 501 PetscCall(MatSOR(user->DSG, X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y)); 502 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "DSG not formed"); 503 PetscFunctionReturn(PETSC_SUCCESS); 504 } 505 506 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y) 507 { 508 AppCtx *user; 509 PetscInt its, i; 510 511 PetscFunctionBegin; 512 PetscCall(MatShellGetContext(J_shell, &user)); 513 514 if (Y == user->ytrue) { 515 /* First solve is done with true solution to set up problem */ 516 PetscCall(KSPSetTolerances(user->solver, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 517 } else { 518 PetscCall(KSPSetTolerances(user->solver, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 519 } 520 521 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); 522 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0])); 523 PetscCall(KSPGetIterationNumber(user->solver, &its)); 524 user->ksp_its = user->ksp_its + its; 525 526 for (i = 1; i < user->nt; i++) { 527 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i - 1])); 528 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 529 PetscCall(KSPGetIterationNumber(user->solver, &its)); 530 user->ksp_its = user->ksp_its + its; 531 } 532 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 533 PetscFunctionReturn(PETSC_SUCCESS); 534 } 535 536 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y) 537 { 538 AppCtx *user; 539 PetscInt its, i; 540 541 PetscFunctionBegin; 542 PetscCall(MatShellGetContext(J_shell, &user)); 543 544 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); 545 546 i = user->nt - 1; 547 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 548 549 PetscCall(KSPGetIterationNumber(user->solver, &its)); 550 user->ksp_its = user->ksp_its + its; 551 552 for (i = user->nt - 2; i >= 0; i--) { 553 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i + 1])); 554 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); 555 556 PetscCall(KSPGetIterationNumber(user->solver, &its)); 557 user->ksp_its = user->ksp_its + its; 558 } 559 560 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); 561 PetscFunctionReturn(PETSC_SUCCESS); 562 } 563 564 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell) 565 { 566 AppCtx *user; 567 568 PetscFunctionBegin; 569 PetscCall(MatShellGetContext(J_shell, &user)); 570 571 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell)); 572 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (void (*)(void))StateMatMult)); 573 PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate)); 574 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose)); 575 PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal)); 576 PetscFunctionReturn(PETSC_SUCCESS); 577 } 578 579 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X) 580 { 581 AppCtx *user; 582 583 PetscFunctionBegin; 584 PetscCall(MatShellGetContext(J_shell, &user)); 585 PetscCall(VecCopy(user->js_diag, X)); 586 PetscFunctionReturn(PETSC_SUCCESS); 587 } 588 589 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr) 590 { 591 /* con = Ay - q, A = [B 0 0 ... 0; 592 -I B 0 ... 0; 593 0 -I B ... 0; 594 ... ; 595 0 ... -I B] 596 B = ht * Div * Sigma * Grad + eye */ 597 PetscInt i; 598 AppCtx *user = (AppCtx *)ptr; 599 600 PetscFunctionBegin; 601 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); 602 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); 603 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); 604 for (i = 1; i < user->nt; i++) { 605 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); 606 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1])); 607 } 608 PetscCall(Gather_i(C, user->yiwork, user->yi_scatter, user->nt)); 609 PetscCall(VecAXPY(C, -1.0, user->q)); 610 PetscFunctionReturn(PETSC_SUCCESS); 611 } 612 613 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 614 { 615 PetscFunctionBegin; 616 PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 617 PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD)); 618 PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 619 PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD)); 620 PetscFunctionReturn(PETSC_SUCCESS); 621 } 622 623 PetscErrorCode Scatter_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 624 { 625 PetscInt i; 626 627 PetscFunctionBegin; 628 for (i = 0; i < nt; i++) { 629 PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 630 PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD)); 631 } 632 PetscFunctionReturn(PETSC_SUCCESS); 633 } 634 635 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 636 { 637 PetscFunctionBegin; 638 PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 639 PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE)); 640 PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 641 PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE)); 642 PetscFunctionReturn(PETSC_SUCCESS); 643 } 644 645 PetscErrorCode Gather_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 646 { 647 PetscInt i; 648 649 PetscFunctionBegin; 650 for (i = 0; i < nt; i++) { 651 PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 652 PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE)); 653 } 654 PetscFunctionReturn(PETSC_SUCCESS); 655 } 656 657 PetscErrorCode ParabolicInitialize(AppCtx *user) 658 { 659 PetscInt m, n, i, j, k, linear_index, istart, iend, iblock, lo, hi, lo2, hi2; 660 Vec XX, YY, ZZ, XXwork, YYwork, ZZwork, UTwork, yi, di, bc; 661 PetscReal *x, *y, *z; 662 PetscReal h, stime; 663 PetscScalar hinv, neg_hinv, half = 0.5, sqrt_beta; 664 PetscInt im, indx1, indx2, indy1, indy2, indz1, indz2, nx, ny, nz; 665 PetscReal xri, yri, zri, xim, yim, zim, dx1, dx2, dy1, dy2, dz1, dz2, Dx, Dy, Dz; 666 PetscScalar v, vx, vy, vz; 667 IS is_from_y, is_to_yi, is_from_d, is_to_di; 668 /* Data locations */ 669 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, 670 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, 671 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}; 672 673 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, 674 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, 675 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}; 676 677 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, 678 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, 679 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}; 680 681 PetscFunctionBegin; 682 PetscCall(PetscMalloc1(user->mx, &x)); 683 PetscCall(PetscMalloc1(user->mx, &y)); 684 PetscCall(PetscMalloc1(user->mx, &z)); 685 user->jformed = PETSC_FALSE; 686 user->dsg_formed = PETSC_FALSE; 687 688 n = user->mx * user->mx * user->mx; 689 m = 3 * user->mx * user->mx * (user->mx - 1); 690 sqrt_beta = PetscSqrtScalar(user->beta); 691 692 user->ksp_its = 0; 693 user->ksp_its_initial = 0; 694 695 stime = (PetscReal)user->nt / user->ns; 696 PetscCall(PetscMalloc1(user->ns, &user->sample_times)); 697 for (i = 0; i < user->ns; i++) user->sample_times[i] = (PetscInt)(stime * ((PetscReal)i + 1.0) - 0.5); 698 699 PetscCall(VecCreate(PETSC_COMM_WORLD, &XX)); 700 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q)); 701 PetscCall(VecSetSizes(XX, PETSC_DECIDE, n)); 702 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt)); 703 PetscCall(VecSetFromOptions(XX)); 704 PetscCall(VecSetFromOptions(user->q)); 705 706 PetscCall(VecDuplicate(XX, &YY)); 707 PetscCall(VecDuplicate(XX, &ZZ)); 708 PetscCall(VecDuplicate(XX, &XXwork)); 709 PetscCall(VecDuplicate(XX, &YYwork)); 710 PetscCall(VecDuplicate(XX, &ZZwork)); 711 PetscCall(VecDuplicate(XX, &UTwork)); 712 PetscCall(VecDuplicate(XX, &user->utrue)); 713 PetscCall(VecDuplicate(XX, &bc)); 714 715 /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */ 716 h = 1.0 / user->mx; 717 hinv = user->mx; 718 neg_hinv = -hinv; 719 720 PetscCall(VecGetOwnershipRange(XX, &istart, &iend)); 721 for (linear_index = istart; linear_index < iend; linear_index++) { 722 i = linear_index % user->mx; 723 j = ((linear_index - i) / user->mx) % user->mx; 724 k = ((linear_index - i) / user->mx - j) / user->mx; 725 vx = h * (i + 0.5); 726 vy = h * (j + 0.5); 727 vz = h * (k + 0.5); 728 PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES)); 729 PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES)); 730 PetscCall(VecSetValues(ZZ, 1, &linear_index, &vz, INSERT_VALUES)); 731 if ((vx < 0.6) && (vx > 0.4) && (vy < 0.6) && (vy > 0.4) && (vy < 0.6) && (vz < 0.6) && (vz > 0.4)) { 732 v = 1000.0; 733 PetscCall(VecSetValues(bc, 1, &linear_index, &v, INSERT_VALUES)); 734 } 735 } 736 737 PetscCall(VecAssemblyBegin(XX)); 738 PetscCall(VecAssemblyEnd(XX)); 739 PetscCall(VecAssemblyBegin(YY)); 740 PetscCall(VecAssemblyEnd(YY)); 741 PetscCall(VecAssemblyBegin(ZZ)); 742 PetscCall(VecAssemblyEnd(ZZ)); 743 PetscCall(VecAssemblyBegin(bc)); 744 PetscCall(VecAssemblyEnd(bc)); 745 746 /* Compute true parameter function 747 ut = 0.5 + exp(-10*((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)) */ 748 PetscCall(VecCopy(XX, XXwork)); 749 PetscCall(VecCopy(YY, YYwork)); 750 PetscCall(VecCopy(ZZ, ZZwork)); 751 752 PetscCall(VecShift(XXwork, -0.5)); 753 PetscCall(VecShift(YYwork, -0.5)); 754 PetscCall(VecShift(ZZwork, -0.5)); 755 756 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork)); 757 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork)); 758 PetscCall(VecPointwiseMult(ZZwork, ZZwork, ZZwork)); 759 760 PetscCall(VecCopy(XXwork, user->utrue)); 761 PetscCall(VecAXPY(user->utrue, 1.0, YYwork)); 762 PetscCall(VecAXPY(user->utrue, 1.0, ZZwork)); 763 PetscCall(VecScale(user->utrue, -10.0)); 764 PetscCall(VecExp(user->utrue)); 765 PetscCall(VecShift(user->utrue, 0.5)); 766 767 PetscCall(VecDestroy(&XX)); 768 PetscCall(VecDestroy(&YY)); 769 PetscCall(VecDestroy(&ZZ)); 770 PetscCall(VecDestroy(&XXwork)); 771 PetscCall(VecDestroy(&YYwork)); 772 PetscCall(VecDestroy(&ZZwork)); 773 PetscCall(VecDestroy(&UTwork)); 774 775 /* Initial guess and reference model */ 776 PetscCall(VecDuplicate(user->utrue, &user->ur)); 777 PetscCall(VecSet(user->ur, 0.0)); 778 779 /* Generate Grad matrix */ 780 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad)); 781 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, m, n)); 782 PetscCall(MatSetFromOptions(user->Grad)); 783 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 2, NULL, 2, NULL)); 784 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 2, NULL)); 785 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend)); 786 787 for (i = istart; i < iend; i++) { 788 if (i < m / 3) { 789 iblock = i / (user->mx - 1); 790 j = iblock * user->mx + (i % (user->mx - 1)); 791 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 792 j = j + 1; 793 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 794 } 795 if (i >= m / 3 && i < 2 * m / 3) { 796 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); 797 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); 798 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 799 j = j + user->mx; 800 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 801 } 802 if (i >= 2 * m / 3) { 803 j = i - 2 * m / 3; 804 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 805 j = j + user->mx * user->mx; 806 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 807 } 808 } 809 810 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY)); 811 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY)); 812 813 /* Generate arithmetic averaging matrix Av */ 814 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Av)); 815 PetscCall(MatSetSizes(user->Av, PETSC_DECIDE, PETSC_DECIDE, m, n)); 816 PetscCall(MatSetFromOptions(user->Av)); 817 PetscCall(MatMPIAIJSetPreallocation(user->Av, 2, NULL, 2, NULL)); 818 PetscCall(MatSeqAIJSetPreallocation(user->Av, 2, NULL)); 819 PetscCall(MatGetOwnershipRange(user->Av, &istart, &iend)); 820 821 for (i = istart; i < iend; i++) { 822 if (i < m / 3) { 823 iblock = i / (user->mx - 1); 824 j = iblock * user->mx + (i % (user->mx - 1)); 825 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 826 j = j + 1; 827 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 828 } 829 if (i >= m / 3 && i < 2 * m / 3) { 830 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); 831 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); 832 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 833 j = j + user->mx; 834 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 835 } 836 if (i >= 2 * m / 3) { 837 j = i - 2 * m / 3; 838 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 839 j = j + user->mx * user->mx; 840 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); 841 } 842 } 843 844 PetscCall(MatAssemblyBegin(user->Av, MAT_FINAL_ASSEMBLY)); 845 PetscCall(MatAssemblyEnd(user->Av, MAT_FINAL_ASSEMBLY)); 846 847 /* Generate transpose of averaging matrix Av */ 848 PetscCall(MatTranspose(user->Av, MAT_INITIAL_MATRIX, &user->AvT)); 849 850 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->L)); 851 PetscCall(MatSetSizes(user->L, PETSC_DECIDE, PETSC_DECIDE, m + n, n)); 852 PetscCall(MatSetFromOptions(user->L)); 853 PetscCall(MatMPIAIJSetPreallocation(user->L, 2, NULL, 2, NULL)); 854 PetscCall(MatSeqAIJSetPreallocation(user->L, 2, NULL)); 855 PetscCall(MatGetOwnershipRange(user->L, &istart, &iend)); 856 857 for (i = istart; i < iend; i++) { 858 if (i < m / 3) { 859 iblock = i / (user->mx - 1); 860 j = iblock * user->mx + (i % (user->mx - 1)); 861 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 862 j = j + 1; 863 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 864 } 865 if (i >= m / 3 && i < 2 * m / 3) { 866 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); 867 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); 868 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 869 j = j + user->mx; 870 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 871 } 872 if (i >= 2 * m / 3 && i < m) { 873 j = i - 2 * m / 3; 874 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); 875 j = j + user->mx * user->mx; 876 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); 877 } 878 if (i >= m) { 879 j = i - m; 880 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &sqrt_beta, INSERT_VALUES)); 881 } 882 } 883 884 PetscCall(MatAssemblyBegin(user->L, MAT_FINAL_ASSEMBLY)); 885 PetscCall(MatAssemblyEnd(user->L, MAT_FINAL_ASSEMBLY)); 886 887 PetscCall(MatScale(user->L, PetscPowScalar(h, 1.5))); 888 889 /* Generate Div matrix */ 890 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div)); 891 892 /* Build work vectors and matrices */ 893 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->S)); 894 PetscCall(VecSetSizes(user->S, PETSC_DECIDE, user->mx * user->mx * (user->mx - 1) * 3)); 895 PetscCall(VecSetFromOptions(user->S)); 896 897 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork)); 898 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, m + user->mx * user->mx * user->mx)); 899 PetscCall(VecSetFromOptions(user->lwork)); 900 901 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork)); 902 PetscCall(MatDuplicate(user->Av, MAT_SHARE_NONZERO_PATTERN, &user->Avwork)); 903 904 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->d)); 905 PetscCall(VecSetSizes(user->d, PETSC_DECIDE, user->ndata * user->nt)); 906 PetscCall(VecSetFromOptions(user->d)); 907 908 PetscCall(VecDuplicate(user->S, &user->Swork)); 909 PetscCall(VecDuplicate(user->S, &user->Av_u)); 910 PetscCall(VecDuplicate(user->S, &user->Twork)); 911 PetscCall(VecDuplicate(user->S, &user->Rwork)); 912 PetscCall(VecDuplicate(user->y, &user->ywork)); 913 PetscCall(VecDuplicate(user->u, &user->uwork)); 914 PetscCall(VecDuplicate(user->u, &user->js_diag)); 915 PetscCall(VecDuplicate(user->c, &user->cwork)); 916 PetscCall(VecDuplicate(user->d, &user->dwork)); 917 918 /* Create matrix-free shell user->Js for computing A*x */ 919 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &user->Js)); 920 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (void (*)(void))StateMatMult)); 921 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate)); 922 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose)); 923 PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal)); 924 925 /* Diagonal blocks of user->Js */ 926 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlock)); 927 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (void (*)(void))StateMatBlockMult)); 928 /* Blocks are symmetric */ 929 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockMult)); 930 931 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U, 932 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular. 933 This is an SSOR preconditioner for user->JsBlock. */ 934 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlockPrec)); 935 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (void (*)(void))StateMatBlockPrecMult)); 936 /* JsBlockPrec is symmetric */ 937 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockPrecMult)); 938 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRIC, PETSC_TRUE)); 939 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRY_ETERNAL, PETSC_TRUE)); 940 941 /* Create a matrix-free shell user->Jd for computing B*x */ 942 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m, user, &user->Jd)); 943 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (void (*)(void))DesignMatMult)); 944 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (void (*)(void))DesignMatMultTranspose)); 945 946 /* User-defined routines for computing user->Js\x and user->Js^T\x*/ 947 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &user->JsInv)); 948 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (void (*)(void))StateMatInvMult)); 949 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatInvTransposeMult)); 950 951 /* Solver options and tolerances */ 952 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver)); 953 PetscCall(KSPSetType(user->solver, KSPCG)); 954 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec)); 955 PetscCall(KSPSetInitialGuessNonzero(user->solver, PETSC_FALSE)); 956 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500)); 957 PetscCall(KSPSetFromOptions(user->solver)); 958 PetscCall(KSPGetPC(user->solver, &user->prec)); 959 PetscCall(PCSetType(user->prec, PCSHELL)); 960 961 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult)); 962 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMult)); 963 PetscCall(PCShellSetContext(user->prec, user)); 964 965 /* Create scatter from y to y_1,y_2,...,y_nt */ 966 PetscCall(PetscMalloc1(user->nt * user->m, &user->yi_scatter)); 967 PetscCall(VecCreate(PETSC_COMM_WORLD, &yi)); 968 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx * user->mx)); 969 PetscCall(VecSetFromOptions(yi)); 970 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi)); 971 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork)); 972 973 PetscCall(VecGetOwnershipRange(user->y, &lo2, &hi2)); 974 istart = 0; 975 for (i = 0; i < user->nt; i++) { 976 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi)); 977 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi)); 978 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo2 + istart, 1, &is_from_y)); 979 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i])); 980 istart = istart + hi - lo; 981 PetscCall(ISDestroy(&is_to_yi)); 982 PetscCall(ISDestroy(&is_from_y)); 983 } 984 PetscCall(VecDestroy(&yi)); 985 986 /* Create scatter from d to d_1,d_2,...,d_ns */ 987 PetscCall(PetscMalloc1(user->ns * user->ndata, &user->di_scatter)); 988 PetscCall(VecCreate(PETSC_COMM_WORLD, &di)); 989 PetscCall(VecSetSizes(di, PETSC_DECIDE, user->ndata)); 990 PetscCall(VecSetFromOptions(di)); 991 PetscCall(VecDuplicateVecs(di, user->ns, &user->di)); 992 PetscCall(VecGetOwnershipRange(user->d, &lo2, &hi2)); 993 istart = 0; 994 for (i = 0; i < user->ns; i++) { 995 PetscCall(VecGetOwnershipRange(user->di[i], &lo, &hi)); 996 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_di)); 997 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo2 + istart, 1, &is_from_d)); 998 PetscCall(VecScatterCreate(user->d, is_from_d, user->di[i], is_to_di, &user->di_scatter[i])); 999 istart = istart + hi - lo; 1000 PetscCall(ISDestroy(&is_to_di)); 1001 PetscCall(ISDestroy(&is_from_d)); 1002 } 1003 PetscCall(VecDestroy(&di)); 1004 1005 /* Assemble RHS of forward problem */ 1006 PetscCall(VecCopy(bc, user->yiwork[0])); 1007 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0)); 1008 PetscCall(Gather_i(user->q, user->yiwork, user->yi_scatter, user->nt)); 1009 PetscCall(VecDestroy(&bc)); 1010 1011 /* Compute true state function yt given ut */ 1012 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue)); 1013 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt)); 1014 PetscCall(VecSetFromOptions(user->ytrue)); 1015 1016 /* First compute Av_u = Av*exp(-u) */ 1017 PetscCall(VecSet(user->uwork, 0)); 1018 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); /* Note: user->utrue */ 1019 PetscCall(VecExp(user->uwork)); 1020 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); 1021 1022 /* Symbolic DSG = Div * Grad */ 1023 PetscCall(MatProductCreate(user->Div, user->Grad, NULL, &user->DSG)); 1024 PetscCall(MatProductSetType(user->DSG, MATPRODUCT_AB)); 1025 PetscCall(MatProductSetAlgorithm(user->DSG, "default")); 1026 PetscCall(MatProductSetFill(user->DSG, PETSC_DEFAULT)); 1027 PetscCall(MatProductSetFromOptions(user->DSG)); 1028 PetscCall(MatProductSymbolic(user->DSG)); 1029 1030 user->dsg_formed = PETSC_TRUE; 1031 1032 /* Next form DSG = Div*Grad */ 1033 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); 1034 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u)); 1035 if (user->dsg_formed) { 1036 PetscCall(MatProductNumeric(user->DSG)); 1037 } else { 1038 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &user->DSG)); 1039 user->dsg_formed = PETSC_TRUE; 1040 } 1041 /* B = speye(nx^3) + ht*DSG; */ 1042 PetscCall(MatScale(user->DSG, user->ht)); 1043 PetscCall(MatShift(user->DSG, 1.0)); 1044 1045 /* Now solve for ytrue */ 1046 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue)); 1047 1048 /* Initial guess y0 for state given u0 */ 1049 1050 /* First compute Av_u = Av*exp(-u) */ 1051 PetscCall(VecSet(user->uwork, 0)); 1052 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); /* Note: user->u */ 1053 PetscCall(VecExp(user->uwork)); 1054 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); 1055 1056 /* Next form DSG = Div*S*Grad */ 1057 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); 1058 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u)); 1059 if (user->dsg_formed) { 1060 PetscCall(MatProductNumeric(user->DSG)); 1061 } else { 1062 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &user->DSG)); 1063 1064 user->dsg_formed = PETSC_TRUE; 1065 } 1066 /* B = speye(nx^3) + ht*DSG; */ 1067 PetscCall(MatScale(user->DSG, user->ht)); 1068 PetscCall(MatShift(user->DSG, 1.0)); 1069 1070 /* Now solve for y */ 1071 PetscCall(StateMatInvMult(user->Js, user->q, user->y)); 1072 1073 /* Construct projection matrix Q, a block diagonal matrix consisting of nt copies of Qblock along the diagonal */ 1074 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Qblock)); 1075 PetscCall(MatSetSizes(user->Qblock, PETSC_DECIDE, PETSC_DECIDE, user->ndata, n)); 1076 PetscCall(MatSetFromOptions(user->Qblock)); 1077 PetscCall(MatMPIAIJSetPreallocation(user->Qblock, 8, NULL, 8, NULL)); 1078 PetscCall(MatSeqAIJSetPreallocation(user->Qblock, 8, NULL)); 1079 1080 for (i = 0; i < user->mx; i++) { 1081 x[i] = h * (i + 0.5); 1082 y[i] = h * (i + 0.5); 1083 z[i] = h * (i + 0.5); 1084 } 1085 1086 PetscCall(MatGetOwnershipRange(user->Qblock, &istart, &iend)); 1087 nx = user->mx; 1088 ny = user->mx; 1089 nz = user->mx; 1090 for (i = istart; i < iend; i++) { 1091 xri = xr[i]; 1092 im = 0; 1093 xim = x[im]; 1094 while (xri > xim && im < nx) { 1095 im = im + 1; 1096 xim = x[im]; 1097 } 1098 indx1 = im - 1; 1099 indx2 = im; 1100 dx1 = xri - x[indx1]; 1101 dx2 = x[indx2] - xri; 1102 1103 yri = yr[i]; 1104 im = 0; 1105 yim = y[im]; 1106 while (yri > yim && im < ny) { 1107 im = im + 1; 1108 yim = y[im]; 1109 } 1110 indy1 = im - 1; 1111 indy2 = im; 1112 dy1 = yri - y[indy1]; 1113 dy2 = y[indy2] - yri; 1114 1115 zri = zr[i]; 1116 im = 0; 1117 zim = z[im]; 1118 while (zri > zim && im < nz) { 1119 im = im + 1; 1120 zim = z[im]; 1121 } 1122 indz1 = im - 1; 1123 indz2 = im; 1124 dz1 = zri - z[indz1]; 1125 dz2 = z[indz2] - zri; 1126 1127 Dx = x[indx2] - x[indx1]; 1128 Dy = y[indy2] - y[indy1]; 1129 Dz = z[indz2] - z[indz1]; 1130 1131 j = indx1 + indy1 * nx + indz1 * nx * ny; 1132 v = (1 - dx1 / Dx) * (1 - dy1 / Dy) * (1 - dz1 / Dz); 1133 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1134 1135 j = indx1 + indy1 * nx + indz2 * nx * ny; 1136 v = (1 - dx1 / Dx) * (1 - dy1 / Dy) * (1 - dz2 / Dz); 1137 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1138 1139 j = indx1 + indy2 * nx + indz1 * nx * ny; 1140 v = (1 - dx1 / Dx) * (1 - dy2 / Dy) * (1 - dz1 / Dz); 1141 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1142 1143 j = indx1 + indy2 * nx + indz2 * nx * ny; 1144 v = (1 - dx1 / Dx) * (1 - dy2 / Dy) * (1 - dz2 / Dz); 1145 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1146 1147 j = indx2 + indy1 * nx + indz1 * nx * ny; 1148 v = (1 - dx2 / Dx) * (1 - dy1 / Dy) * (1 - dz1 / Dz); 1149 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1150 1151 j = indx2 + indy1 * nx + indz2 * nx * ny; 1152 v = (1 - dx2 / Dx) * (1 - dy1 / Dy) * (1 - dz2 / Dz); 1153 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1154 1155 j = indx2 + indy2 * nx + indz1 * nx * ny; 1156 v = (1 - dx2 / Dx) * (1 - dy2 / Dy) * (1 - dz1 / Dz); 1157 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1158 1159 j = indx2 + indy2 * nx + indz2 * nx * ny; 1160 v = (1 - dx2 / Dx) * (1 - dy2 / Dy) * (1 - dz2 / Dz); 1161 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); 1162 } 1163 PetscCall(MatAssemblyBegin(user->Qblock, MAT_FINAL_ASSEMBLY)); 1164 PetscCall(MatAssemblyEnd(user->Qblock, MAT_FINAL_ASSEMBLY)); 1165 1166 PetscCall(MatTranspose(user->Qblock, MAT_INITIAL_MATRIX, &user->QblockT)); 1167 PetscCall(MatTranspose(user->L, MAT_INITIAL_MATRIX, &user->LT)); 1168 1169 /* Add noise to the measurement data */ 1170 PetscCall(VecSet(user->ywork, 1.0)); 1171 PetscCall(VecAYPX(user->ywork, user->noise, user->ytrue)); 1172 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); 1173 for (j = 0; j < user->ns; j++) { 1174 i = user->sample_times[j]; 1175 PetscCall(MatMult(user->Qblock, user->yiwork[i], user->di[j])); 1176 } 1177 PetscCall(Gather_i(user->d, user->di, user->di_scatter, user->ns)); 1178 1179 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */ 1180 PetscCall(KSPSetFromOptions(user->solver)); 1181 PetscCall(PetscFree(x)); 1182 PetscCall(PetscFree(y)); 1183 PetscCall(PetscFree(z)); 1184 PetscFunctionReturn(PETSC_SUCCESS); 1185 } 1186 1187 PetscErrorCode ParabolicDestroy(AppCtx *user) 1188 { 1189 PetscInt i; 1190 1191 PetscFunctionBegin; 1192 PetscCall(MatDestroy(&user->Qblock)); 1193 PetscCall(MatDestroy(&user->QblockT)); 1194 PetscCall(MatDestroy(&user->Div)); 1195 PetscCall(MatDestroy(&user->Divwork)); 1196 PetscCall(MatDestroy(&user->Grad)); 1197 PetscCall(MatDestroy(&user->Av)); 1198 PetscCall(MatDestroy(&user->Avwork)); 1199 PetscCall(MatDestroy(&user->AvT)); 1200 PetscCall(MatDestroy(&user->DSG)); 1201 PetscCall(MatDestroy(&user->L)); 1202 PetscCall(MatDestroy(&user->LT)); 1203 PetscCall(KSPDestroy(&user->solver)); 1204 PetscCall(MatDestroy(&user->Js)); 1205 PetscCall(MatDestroy(&user->Jd)); 1206 PetscCall(MatDestroy(&user->JsInv)); 1207 PetscCall(MatDestroy(&user->JsBlock)); 1208 PetscCall(MatDestroy(&user->JsBlockPrec)); 1209 PetscCall(VecDestroy(&user->u)); 1210 PetscCall(VecDestroy(&user->uwork)); 1211 PetscCall(VecDestroy(&user->utrue)); 1212 PetscCall(VecDestroy(&user->y)); 1213 PetscCall(VecDestroy(&user->ywork)); 1214 PetscCall(VecDestroy(&user->ytrue)); 1215 PetscCall(VecDestroyVecs(user->nt, &user->yi)); 1216 PetscCall(VecDestroyVecs(user->nt, &user->yiwork)); 1217 PetscCall(VecDestroyVecs(user->ns, &user->di)); 1218 PetscCall(PetscFree(user->yi)); 1219 PetscCall(PetscFree(user->yiwork)); 1220 PetscCall(PetscFree(user->di)); 1221 PetscCall(VecDestroy(&user->c)); 1222 PetscCall(VecDestroy(&user->cwork)); 1223 PetscCall(VecDestroy(&user->ur)); 1224 PetscCall(VecDestroy(&user->q)); 1225 PetscCall(VecDestroy(&user->d)); 1226 PetscCall(VecDestroy(&user->dwork)); 1227 PetscCall(VecDestroy(&user->lwork)); 1228 PetscCall(VecDestroy(&user->S)); 1229 PetscCall(VecDestroy(&user->Swork)); 1230 PetscCall(VecDestroy(&user->Av_u)); 1231 PetscCall(VecDestroy(&user->Twork)); 1232 PetscCall(VecDestroy(&user->Rwork)); 1233 PetscCall(VecDestroy(&user->js_diag)); 1234 PetscCall(ISDestroy(&user->s_is)); 1235 PetscCall(ISDestroy(&user->d_is)); 1236 PetscCall(VecScatterDestroy(&user->state_scatter)); 1237 PetscCall(VecScatterDestroy(&user->design_scatter)); 1238 for (i = 0; i < user->nt; i++) PetscCall(VecScatterDestroy(&user->yi_scatter[i])); 1239 for (i = 0; i < user->ns; i++) PetscCall(VecScatterDestroy(&user->di_scatter[i])); 1240 PetscCall(PetscFree(user->yi_scatter)); 1241 PetscCall(PetscFree(user->di_scatter)); 1242 PetscCall(PetscFree(user->sample_times)); 1243 PetscFunctionReturn(PETSC_SUCCESS); 1244 } 1245 1246 PetscErrorCode ParabolicMonitor(Tao tao, void *ptr) 1247 { 1248 Vec X; 1249 PetscReal unorm, ynorm; 1250 AppCtx *user = (AppCtx *)ptr; 1251 1252 PetscFunctionBegin; 1253 PetscCall(TaoGetSolution(tao, &X)); 1254 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); 1255 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue)); 1256 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); 1257 PetscCall(VecNorm(user->uwork, NORM_2, &unorm)); 1258 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm)); 1259 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm)); 1260 PetscFunctionReturn(PETSC_SUCCESS); 1261 } 1262 1263 /*TEST 1264 1265 build: 1266 requires: !complex 1267 1268 test: 1269 args: -tao_monitor_constraint_norm -tao_type lcl -ns 1 -tao_gatol 1.e-4 -ksp_max_it 30 1270 requires: !single 1271 1272 TEST*/ 1273