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