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