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