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