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