Actual source code: elliptic.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/taoimpl.h>
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: TaoSetHistory(); TaoGetHistory();
16: Routines: TaoSolve();
17: Routines: TaoDestroy();
18: Processors: n
19: T*/
21: typedef struct {
22: PetscInt n; /* Number of total variables */
23: PetscInt m; /* Number of constraints */
24: PetscInt nstate;
25: PetscInt ndesign;
26: PetscInt mx; /* grid points in each direction */
27: PetscInt ns; /* Number of data samples (1<=ns<=8)
28: Currently only ns=1 is supported */
29: PetscInt ndata; /* Number of data points per sample */
30: IS s_is;
31: IS d_is;
33: VecScatter state_scatter;
34: VecScatter design_scatter;
35: VecScatter *yi_scatter, *di_scatter;
36: Vec suby,subq,subd;
37: Mat Js,Jd,JsPrec,JsInv,JsBlock;
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 *ones;
43: Mat Q;
44: Mat MQ;
45: Mat L;
47: Mat Grad;
48: Mat Av,Avwork;
49: Mat Div, Divwork;
50: Mat DSG;
51: Mat Diag,Ones;
54: Vec q;
55: Vec ur; /* reference */
57: Vec d;
58: Vec dwork;
60: Vec x; /* super vec of y,u */
62: Vec y; /* state variables */
63: Vec ywork;
65: Vec ytrue;
67: Vec u; /* design variables */
68: Vec uwork;
70: Vec utrue;
72: Vec js_diag;
74: Vec c; /* constraint vector */
75: Vec cwork;
77: Vec lwork;
78: Vec S;
79: Vec Swork,Twork,Sdiag,Ywork;
80: Vec Av_u;
82: KSP solver;
83: PC prec;
85: PetscReal tola,tolb,tolc,told;
86: PetscInt ksp_its;
87: PetscInt ksp_its_initial;
88: PetscLogStage stages[10];
89: PetscBool use_ptap;
90: PetscBool use_lrc;
91: } AppCtx;
93: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
94: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
95: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
96: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
97: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
98: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
99: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
100: PetscErrorCode Gather(Vec, Vec, VecScatter, Vec, VecScatter);
101: PetscErrorCode Scatter(Vec, Vec, VecScatter, Vec, VecScatter);
102: PetscErrorCode EllipticInitialize(AppCtx*);
103: PetscErrorCode EllipticDestroy(AppCtx*);
104: PetscErrorCode EllipticMonitor(Tao, void*);
106: PetscErrorCode StateBlockMatMult(Mat,Vec,Vec);
107: PetscErrorCode StateMatMult(Mat,Vec,Vec);
109: PetscErrorCode StateInvMatMult(Mat,Vec,Vec);
110: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
111: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
113: PetscErrorCode QMatMult(Mat,Vec,Vec);
114: PetscErrorCode QMatMultTranspose(Mat,Vec,Vec);
116: static char help[]="";
118: int main(int argc, char **argv)
119: {
120: PetscErrorCode ierr;
121: Vec x0;
122: Tao tao;
123: AppCtx user;
124: PetscInt ntests = 1;
125: PetscInt i;
127: PetscInitialize(&argc, &argv, (char*)0,help);if (ierr) return ierr;
128: user.mx = 8;
129: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
130: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
131: user.ns = 6;
132: PetscOptionsInt("-ns","Number of data samples (1<=ns<=8)","",user.ns,&user.ns,NULL);
133: user.ndata = 64;
134: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
135: user.alpha = 0.1;
136: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
137: user.beta = 0.00001;
138: PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL);
139: user.noise = 0.01;
140: PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL);
142: user.use_ptap = PETSC_FALSE;
143: PetscOptionsBool("-use_ptap","Use ptap matrix for DSG","",user.use_ptap,&user.use_ptap,NULL);
144: user.use_lrc = PETSC_FALSE;
145: PetscOptionsBool("-use_lrc","Use lrc matrix for Js","",user.use_lrc,&user.use_lrc,NULL);
146: user.m = user.ns*user.mx*user.mx*user.mx; /* number of constraints */
147: user.nstate = user.m;
148: user.ndesign = user.mx*user.mx*user.mx;
149: user.n = user.nstate + user.ndesign; /* number of variables */
151: /* Create TAO solver and set desired solution method */
152: TaoCreate(PETSC_COMM_WORLD,&tao);
153: TaoSetType(tao,TAOLCL);
155: /* Set up initial vectors and matrices */
156: EllipticInitialize(&user);
158: Gather(user.x,user.y,user.state_scatter,user.u,user.design_scatter);
159: VecDuplicate(user.x,&x0);
160: VecCopy(user.x,x0);
162: /* Set solution vector with an initial guess */
163: TaoSetInitialVector(tao,user.x);
164: TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
165: TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
166: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);
168: TaoSetJacobianStateRoutine(tao, user.Js, NULL, user.JsInv, FormJacobianState, (void *)&user);
169: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);
171: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
172: TaoSetFromOptions(tao);
174: /* SOLVE THE APPLICATION */
175: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
176: PetscOptionsEnd();
177: PetscLogStageRegister("Trials",&user.stages[1]);
178: PetscLogStagePush(user.stages[1]);
179: for (i=0; i<ntests; i++){
180: TaoSolve(tao);
181: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its);
182: VecCopy(x0,user.x);
183: }
184: PetscLogStagePop();
185: PetscBarrier((PetscObject)user.x);
186: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
187: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
189: TaoDestroy(&tao);
190: VecDestroy(&x0);
191: EllipticDestroy(&user);
192: PetscFinalize();
193: return ierr;
194: }
195: /* ------------------------------------------------------------------- */
196: /*
197: dwork = Qy - d
198: lwork = L*(u-ur)
199: f = 1/2 * (dwork.dwork + alpha*lwork.lwork)
200: */
201: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
202: {
204: PetscReal d1=0,d2=0;
205: AppCtx *user = (AppCtx*)ptr;
208: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
209: MatMult(user->MQ,user->y,user->dwork);
210: VecAXPY(user->dwork,-1.0,user->d);
211: VecDot(user->dwork,user->dwork,&d1);
212: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
213: MatMult(user->L,user->uwork,user->lwork);
214: VecDot(user->lwork,user->lwork,&d2);
215: *f = 0.5 * (d1 + user->alpha*d2);
216: return(0);
217: }
219: /* ------------------------------------------------------------------- */
220: /*
221: state: g_s = Q' *(Qy - d)
222: design: g_d = alpha*L'*L*(u-ur)
223: */
224: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
225: {
227: AppCtx *user = (AppCtx*)ptr;
230: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
231: MatMult(user->MQ,user->y,user->dwork);
232: VecAXPY(user->dwork,-1.0,user->d);
233: MatMultTranspose(user->MQ,user->dwork,user->ywork);
234: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
235: MatMult(user->L,user->uwork,user->lwork);
236: MatMultTranspose(user->L,user->lwork,user->uwork);
237: VecScale(user->uwork, user->alpha);
238: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
239: return(0);
240: }
242: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
243: {
245: PetscReal d1,d2;
246: AppCtx *user = (AppCtx*)ptr;
249: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
250: MatMult(user->MQ,user->y,user->dwork);
251: VecAXPY(user->dwork,-1.0,user->d);
252: VecDot(user->dwork,user->dwork,&d1);
253: MatMultTranspose(user->MQ,user->dwork,user->ywork);
255: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
256: MatMult(user->L,user->uwork,user->lwork);
257: VecDot(user->lwork,user->lwork,&d2);
258: MatMultTranspose(user->L,user->lwork,user->uwork);
259: VecScale(user->uwork, user->alpha);
260: *f = 0.5 * (d1 + user->alpha*d2);
261: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
262: return(0);
263: }
265: /* ------------------------------------------------------------------- */
266: /* A
267: MatShell object
268: */
269: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
270: {
272: AppCtx *user = (AppCtx*)ptr;
275: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
276: /* DSG = Div * (1/Av_u) * Grad */
277: VecSet(user->uwork,0);
278: VecAXPY(user->uwork,-1.0,user->u);
279: VecExp(user->uwork);
280: MatMult(user->Av,user->uwork,user->Av_u);
281: VecCopy(user->Av_u,user->Swork);
282: VecReciprocal(user->Swork);
283: if (user->use_ptap) {
284: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
285: MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
286: } else {
287: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
288: MatDiagonalScale(user->Divwork,NULL,user->Swork);
289: MatMatMultNumeric(user->Divwork,user->Grad,user->DSG);
290: }
291: return(0);
292: }
293: /* ------------------------------------------------------------------- */
294: /* B */
295: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
296: {
298: AppCtx *user = (AppCtx*)ptr;
301: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
302: return(0);
303: }
305: PetscErrorCode StateBlockMatMult(Mat J_shell, Vec X, Vec Y)
306: {
308: PetscReal sum;
309: AppCtx *user;
312: MatShellGetContext(J_shell,(void**)&user);
313: MatMult(user->DSG,X,Y);
314: VecSum(X,&sum);
315: sum /= user->ndesign;
316: VecShift(Y,sum);
317: return(0);
318: }
320: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
321: {
323: PetscInt i;
324: AppCtx *user;
327: MatShellGetContext(J_shell,(void**)&user);
328: if (user->ns == 1) {
329: MatMult(user->JsBlock,X,Y);
330: } else {
331: for (i=0;i<user->ns;i++) {
332: Scatter(X,user->subq,user->yi_scatter[i],0,0);
333: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
334: MatMult(user->JsBlock,user->subq,user->suby);
335: Gather(Y,user->suby,user->yi_scatter[i],0,0);
336: }
337: }
338: return(0);
339: }
341: PetscErrorCode StateInvMatMult(Mat J_shell, Vec X, Vec Y)
342: {
344: PetscInt its,i;
345: AppCtx *user;
348: MatShellGetContext(J_shell,(void**)&user);
349: KSPSetOperators(user->solver,user->JsBlock,user->DSG);
350: if (Y == user->ytrue) {
351: /* First solve is done using true solution to set up problem */
352: KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
353: } else {
354: KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
355: }
356: if (user->ns == 1) {
357: KSPSolve(user->solver,X,Y);
358: KSPGetIterationNumber(user->solver,&its);
359: user->ksp_its+=its;
360: } else {
361: for (i=0;i<user->ns;i++) {
362: Scatter(X,user->subq,user->yi_scatter[i],0,0);
363: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
364: KSPSolve(user->solver,user->subq,user->suby);
365: KSPGetIterationNumber(user->solver,&its);
366: user->ksp_its+=its;
367: Gather(Y,user->suby,user->yi_scatter[i],0,0);
368: }
369: }
370: return(0);
371: }
372: PetscErrorCode QMatMult(Mat J_shell, Vec X, Vec Y)
373: {
375: AppCtx *user;
376: PetscInt i;
379: MatShellGetContext(J_shell,(void**)&user);
380: if (user->ns == 1) {
381: MatMult(user->Q,X,Y);
382: } else {
383: for (i=0;i<user->ns;i++) {
384: Scatter(X,user->subq,user->yi_scatter[i],0,0);
385: Scatter(Y,user->subd,user->di_scatter[i],0,0);
386: MatMult(user->Q,user->subq,user->subd);
387: Gather(Y,user->subd,user->di_scatter[i],0,0);
388: }
389: }
390: return(0);
391: }
393: PetscErrorCode QMatMultTranspose(Mat J_shell, Vec X, Vec Y)
394: {
396: AppCtx *user;
397: PetscInt i;
400: MatShellGetContext(J_shell,(void**)&user);
401: if (user->ns == 1) {
402: MatMultTranspose(user->Q,X,Y);
403: } else {
404: for (i=0;i<user->ns;i++) {
405: Scatter(X,user->subd,user->di_scatter[i],0,0);
406: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
407: MatMultTranspose(user->Q,user->subd,user->suby);
408: Gather(Y,user->suby,user->yi_scatter[i],0,0);
409: }
410: }
411: return(0);
412: }
414: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
415: {
417: PetscInt i;
418: AppCtx *user;
421: MatShellGetContext(J_shell,(void**)&user);
423: /* sdiag(1./v) */
424: VecSet(user->uwork,0);
425: VecAXPY(user->uwork,-1.0,user->u);
426: VecExp(user->uwork);
428: /* sdiag(1./((Av*(1./v)).^2)) */
429: MatMult(user->Av,user->uwork,user->Swork);
430: VecPointwiseMult(user->Swork,user->Swork,user->Swork);
431: VecReciprocal(user->Swork);
433: /* (Av * (sdiag(1./v) * b)) */
434: VecPointwiseMult(user->uwork,user->uwork,X);
435: MatMult(user->Av,user->uwork,user->Twork);
437: /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */
438: VecPointwiseMult(user->Swork,user->Twork,user->Swork);
440: if (user->ns == 1) {
441: /* (sdiag(Grad*y(:,i)) */
442: MatMult(user->Grad,user->y,user->Twork);
444: /* Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */
445: VecPointwiseMult(user->Swork,user->Twork,user->Swork);
446: MatMultTranspose(user->Grad,user->Swork,Y);
447: } else {
448: for (i=0;i<user->ns;i++) {
449: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
450: Scatter(Y,user->subq,user->yi_scatter[i],0,0);
452: MatMult(user->Grad,user->suby,user->Twork);
453: VecPointwiseMult(user->Twork,user->Twork,user->Swork);
454: MatMultTranspose(user->Grad,user->Twork,user->subq);
455: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
456: Gather(Y,user->subq,user->yi_scatter[i],0,0);
457: }
458: }
459: return(0);
460: }
462: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
463: {
465: PetscInt i;
466: AppCtx *user;
469: MatShellGetContext(J_shell,(void**)&user);
470: VecZeroEntries(Y);
472: /* Sdiag = 1./((Av*(1./v)).^2) */
473: VecSet(user->uwork,0);
474: VecAXPY(user->uwork,-1.0,user->u);
475: VecExp(user->uwork);
476: MatMult(user->Av,user->uwork,user->Swork);
477: VecPointwiseMult(user->Sdiag,user->Swork,user->Swork);
478: VecReciprocal(user->Sdiag);
480: for (i=0;i<user->ns;i++) {
481: Scatter(X,user->subq,user->yi_scatter[i],0,0);
482: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
484: /* Swork = (Div' * b(:,i)) */
485: MatMult(user->Grad,user->subq,user->Swork);
487: /* Twork = Grad*y(:,i) */
488: MatMult(user->Grad,user->suby,user->Twork);
490: /* Twork = sdiag(Twork) * Swork */
491: VecPointwiseMult(user->Twork,user->Swork,user->Twork);
494: /* Swork = pointwisemult(Sdiag,Twork) */
495: VecPointwiseMult(user->Swork,user->Twork,user->Sdiag);
497: /* Ywork = Av' * Swork */
498: MatMultTranspose(user->Av,user->Swork,user->Ywork);
500: /* Ywork = pointwisemult(uwork,Ywork) */
501: VecPointwiseMult(user->Ywork,user->uwork,user->Ywork);
502: VecAXPY(Y,1.0,user->Ywork);
503: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
504: }
505: return(0);
506: }
508: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
509: {
510: /* C=Ay - q A = Div * Sigma * Grad + hx*hx*hx*ones(n,n) */
512: PetscReal sum;
513: PetscInt i;
514: AppCtx *user = (AppCtx*)ptr;
517: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
518: if (user->ns == 1) {
519: MatMult(user->Grad,user->y,user->Swork);
520: VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
521: MatMultTranspose(user->Grad,user->Swork,C);
522: VecSum(user->y,&sum);
523: sum /= user->ndesign;
524: VecShift(C,sum);
525: } else {
526: for (i=0;i<user->ns;i++) {
527: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
528: Scatter(C,user->subq,user->yi_scatter[i],0,0);
529: MatMult(user->Grad,user->suby,user->Swork);
530: VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
531: MatMultTranspose(user->Grad,user->Swork,user->subq);
533: VecSum(user->suby,&sum);
534: sum /= user->ndesign;
535: VecShift(user->subq,sum);
537: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
538: Gather(C,user->subq,user->yi_scatter[i],0,0);
539: }
540: }
541: VecAXPY(C,-1.0,user->q);
542: return(0);
543: }
545: PetscErrorCode Scatter(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
546: {
550: VecScatterBegin(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
551: VecScatterEnd(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
552: if (sub2) {
553: VecScatterBegin(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
554: VecScatterEnd(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
555: }
556: return(0);
557: }
559: PetscErrorCode Gather(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
560: {
564: VecScatterBegin(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
565: VecScatterEnd(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
566: if (sub2) {
567: VecScatterBegin(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
568: VecScatterEnd(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
569: }
570: return(0);
571: }
573: PetscErrorCode EllipticInitialize(AppCtx *user)
574: {
576: PetscInt m,n,i,j,k,l,linear_index,is,js,ks,ls,istart,iend,iblock;
577: Vec XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork;
578: PetscReal *x,*y,*z;
579: PetscReal h,meanut;
580: PetscScalar hinv,neg_hinv,half = 0.5,sqrt_beta;
581: PetscInt im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz;
582: IS is_alldesign,is_allstate;
583: IS is_from_d;
584: IS is_from_y;
585: PetscInt lo,hi,hi2,lo2,ysubnlocal,dsubnlocal;
586: const PetscInt *ranges, *subranges;
587: PetscMPIInt size;
588: PetscReal xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz;
589: PetscScalar v,vx,vy,vz;
590: PetscInt offset,subindex,subvec,nrank,kk;
592: PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160,
593: 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774,
594: 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997, 0.5198, 0.8326,
595: 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728,
596: 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273,
597: 0.5724, 0.4331, 0.5136, 0.3547, 0.4413, 0.2602, 0.5698, 0.7278,
598: 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594,
599: 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756};
601: PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256,
602: 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792,
603: 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437, 0.5938, 0.6137,
604: 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985,
605: 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288,
606: 0.5761, 0.1129, 0.3841, 0.6150, 0.6904, 0.6686, 0.1361, 0.4601,
607: 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480,
608: 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658};
610: PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295,
611: 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819,
612: 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236, 0.8800, 0.2939,
613: 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712,
614: 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078,
615: 0.1987, 0.2297, 0.4321, 0.8115, 0.4915, 0.7764, 0.4657, 0.4627,
616: 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313,
617: 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711};
620: MPI_Comm_size(PETSC_COMM_WORLD,&size);
621: PetscLogStageRegister("Elliptic Setup",&user->stages[0]);
622: PetscLogStagePush(user->stages[0]);
624: /* Create u,y,c,x */
625: VecCreate(PETSC_COMM_WORLD,&user->u);
626: VecCreate(PETSC_COMM_WORLD,&user->y);
627: VecCreate(PETSC_COMM_WORLD,&user->c);
628: VecSetSizes(user->u,PETSC_DECIDE,user->ndesign);
629: VecSetFromOptions(user->u);
630: VecGetLocalSize(user->u,&ysubnlocal);
631: VecSetSizes(user->y,ysubnlocal*user->ns,user->nstate);
632: VecSetSizes(user->c,ysubnlocal*user->ns,user->m);
633: VecSetFromOptions(user->y);
634: VecSetFromOptions(user->c);
636: /*
637: *******************************
638: Create scatters for x <-> y,u
639: *******************************
641: If the state vector y and design vector u are partitioned as
642: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
643: then the solution vector x is organized as
644: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
645: The index sets user->s_is and user->d_is correspond to the indices of the
646: state and design variables owned by the current processor.
647: */
648: VecCreate(PETSC_COMM_WORLD,&user->x);
650: VecGetOwnershipRange(user->y,&lo,&hi);
651: VecGetOwnershipRange(user->u,&lo2,&hi2);
653: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
654: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user->s_is);
655: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
656: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user->d_is);
658: VecSetSizes(user->x,hi-lo+hi2-lo2,user->n);
659: VecSetFromOptions(user->x);
661: VecScatterCreate(user->x,user->s_is,user->y,is_allstate,&user->state_scatter);
662: VecScatterCreate(user->x,user->d_is,user->u,is_alldesign,&user->design_scatter);
663: ISDestroy(&is_alldesign);
664: ISDestroy(&is_allstate);
665: /*
666: *******************************
667: Create scatter from y to y_1,y_2,...,y_ns
668: *******************************
669: */
670: PetscMalloc1(user->ns,&user->yi_scatter);
671: VecDuplicate(user->u,&user->suby);
672: VecDuplicate(user->u,&user->subq);
674: VecGetOwnershipRange(user->y,&lo2,&hi2);
675: istart = 0;
676: for (i=0; i<user->ns; i++){
677: VecGetOwnershipRange(user->suby,&lo,&hi);
678: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y);
679: VecScatterCreate(user->y,is_from_y,user->suby,NULL,&user->yi_scatter[i]);
680: istart = istart + hi-lo;
681: ISDestroy(&is_from_y);
682: }
683: /*
684: *******************************
685: Create scatter from d to d_1,d_2,...,d_ns
686: *******************************
687: */
688: VecCreate(PETSC_COMM_WORLD,&user->subd);
689: VecSetSizes(user->subd,PETSC_DECIDE,user->ndata);
690: VecSetFromOptions(user->subd);
691: VecCreate(PETSC_COMM_WORLD,&user->d);
692: VecGetLocalSize(user->subd,&dsubnlocal);
693: VecSetSizes(user->d,dsubnlocal*user->ns,user->ndata*user->ns);
694: VecSetFromOptions(user->d);
695: PetscMalloc1(user->ns,&user->di_scatter);
697: VecGetOwnershipRange(user->d,&lo2,&hi2);
698: istart = 0;
699: for (i=0; i<user->ns; i++){
700: VecGetOwnershipRange(user->subd,&lo,&hi);
701: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d);
702: VecScatterCreate(user->d,is_from_d,user->subd,NULL,&user->di_scatter[i]);
703: istart = istart + hi-lo;
704: ISDestroy(&is_from_d);
705: }
707: PetscMalloc1(user->mx,&x);
708: PetscMalloc1(user->mx,&y);
709: PetscMalloc1(user->mx,&z);
711: user->ksp_its = 0;
712: user->ksp_its_initial = 0;
714: n = user->mx * user->mx * user->mx;
715: m = 3 * user->mx * user->mx * (user->mx-1);
716: sqrt_beta = PetscSqrtScalar(user->beta);
718: VecCreate(PETSC_COMM_WORLD,&XX);
719: VecCreate(PETSC_COMM_WORLD,&user->q);
720: VecSetSizes(XX,ysubnlocal,n);
721: VecSetSizes(user->q,ysubnlocal*user->ns,user->m);
722: VecSetFromOptions(XX);
723: VecSetFromOptions(user->q);
725: VecDuplicate(XX,&YY);
726: VecDuplicate(XX,&ZZ);
727: VecDuplicate(XX,&XXwork);
728: VecDuplicate(XX,&YYwork);
729: VecDuplicate(XX,&ZZwork);
730: VecDuplicate(XX,&UTwork);
731: VecDuplicate(XX,&user->utrue);
733: /* map for striding q */
734: VecGetOwnershipRanges(user->q,&ranges);
735: VecGetOwnershipRanges(user->u,&subranges);
737: VecGetOwnershipRange(user->q,&lo2,&hi2);
738: VecGetOwnershipRange(user->u,&lo,&hi);
739: /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */
740: h = 1.0/user->mx;
741: hinv = user->mx;
742: neg_hinv = -hinv;
744: VecGetOwnershipRange(XX,&istart,&iend);
745: for (linear_index=istart; linear_index<iend; linear_index++){
746: i = linear_index % user->mx;
747: j = ((linear_index-i)/user->mx) % user->mx;
748: k = ((linear_index-i)/user->mx-j) / user->mx;
749: vx = h*(i+0.5);
750: vy = h*(j+0.5);
751: vz = h*(k+0.5);
752: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
753: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
754: VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES);
755: for (is=0; is<2; is++){
756: for (js=0; js<2; js++){
757: for (ks=0; ks<2; ks++){
758: ls = is*4 + js*2 + ks;
759: if (ls<user->ns){
760: l =ls*n + linear_index;
761: /* remap */
762: subindex = l%n;
763: subvec = l/n;
764: nrank=0;
765: while (subindex >= subranges[nrank+1]) nrank++;
766: offset = subindex - subranges[nrank];
767: istart=0;
768: for (kk=0;kk<nrank;kk++) istart+=user->ns*(subranges[kk+1]-subranges[kk]);
769: istart += (subranges[nrank+1]-subranges[nrank])*subvec;
770: l = istart+offset;
771: v = 100*PetscSinScalar(2*PETSC_PI*(vx+0.25*is))*PetscSinScalar(2*PETSC_PI*(vy+0.25*js))*PetscSinScalar(2*PETSC_PI*(vz+0.25*ks));
772: VecSetValues(user->q,1,&l,&v,INSERT_VALUES);
773: }
774: }
775: }
776: }
777: }
779: VecAssemblyBegin(XX);
780: VecAssemblyEnd(XX);
781: VecAssemblyBegin(YY);
782: VecAssemblyEnd(YY);
783: VecAssemblyBegin(ZZ);
784: VecAssemblyEnd(ZZ);
785: VecAssemblyBegin(user->q);
786: VecAssemblyEnd(user->q);
788: /* Compute true parameter function
789: ut = exp(-((x-0.25)^2+(y-0.25)^2+(z-0.25)^2)/0.05) - exp((x-0.75)^2-(y-0.75)^2-(z-0.75))^2/0.05) */
790: VecCopy(XX,XXwork);
791: VecCopy(YY,YYwork);
792: VecCopy(ZZ,ZZwork);
794: VecShift(XXwork,-0.25);
795: VecShift(YYwork,-0.25);
796: VecShift(ZZwork,-0.25);
798: VecPointwiseMult(XXwork,XXwork,XXwork);
799: VecPointwiseMult(YYwork,YYwork,YYwork);
800: VecPointwiseMult(ZZwork,ZZwork,ZZwork);
802: VecCopy(XXwork,UTwork);
803: VecAXPY(UTwork,1.0,YYwork);
804: VecAXPY(UTwork,1.0,ZZwork);
805: VecScale(UTwork,-20.0);
806: VecExp(UTwork);
807: VecCopy(UTwork,user->utrue);
809: VecCopy(XX,XXwork);
810: VecCopy(YY,YYwork);
811: VecCopy(ZZ,ZZwork);
813: VecShift(XXwork,-0.75);
814: VecShift(YYwork,-0.75);
815: VecShift(ZZwork,-0.75);
817: VecPointwiseMult(XXwork,XXwork,XXwork);
818: VecPointwiseMult(YYwork,YYwork,YYwork);
819: VecPointwiseMult(ZZwork,ZZwork,ZZwork);
821: VecCopy(XXwork,UTwork);
822: VecAXPY(UTwork,1.0,YYwork);
823: VecAXPY(UTwork,1.0,ZZwork);
824: VecScale(UTwork,-20.0);
825: VecExp(UTwork);
827: VecAXPY(user->utrue,-1.0,UTwork);
829: VecDestroy(&XX);
830: VecDestroy(&YY);
831: VecDestroy(&ZZ);
832: VecDestroy(&XXwork);
833: VecDestroy(&YYwork);
834: VecDestroy(&ZZwork);
835: VecDestroy(&UTwork);
837: /* Initial guess and reference model */
838: VecDuplicate(user->utrue,&user->ur);
839: VecSum(user->utrue,&meanut);
840: meanut = meanut / n;
841: VecSet(user->ur,meanut);
842: VecCopy(user->ur,user->u);
844: /* Generate Grad matrix */
845: MatCreate(PETSC_COMM_WORLD,&user->Grad);
846: MatSetSizes(user->Grad,PETSC_DECIDE,ysubnlocal,m,n);
847: MatSetFromOptions(user->Grad);
848: MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL);
849: MatSeqAIJSetPreallocation(user->Grad,2,NULL);
850: MatGetOwnershipRange(user->Grad,&istart,&iend);
852: for (i=istart; i<iend; i++){
853: if (i<m/3){
854: iblock = i / (user->mx-1);
855: j = iblock*user->mx + (i % (user->mx-1));
856: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
857: j = j+1;
858: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
859: }
860: if (i>=m/3 && i<2*m/3){
861: iblock = (i-m/3) / (user->mx*(user->mx-1));
862: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
863: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
864: j = j + user->mx;
865: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
866: }
867: if (i>=2*m/3){
868: j = i-2*m/3;
869: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
870: j = j + user->mx*user->mx;
871: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
872: }
873: }
875: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
876: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
878: /* Generate arithmetic averaging matrix Av */
879: MatCreate(PETSC_COMM_WORLD,&user->Av);
880: MatSetSizes(user->Av,PETSC_DECIDE,ysubnlocal,m,n);
881: MatSetFromOptions(user->Av);
882: MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL);
883: MatSeqAIJSetPreallocation(user->Av,2,NULL);
884: MatGetOwnershipRange(user->Av,&istart,&iend);
886: for (i=istart; i<iend; i++){
887: if (i<m/3){
888: iblock = i / (user->mx-1);
889: j = iblock*user->mx + (i % (user->mx-1));
890: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
891: j = j+1;
892: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
893: }
894: if (i>=m/3 && i<2*m/3){
895: iblock = (i-m/3) / (user->mx*(user->mx-1));
896: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
897: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
898: j = j + user->mx;
899: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
900: }
901: if (i>=2*m/3){
902: j = i-2*m/3;
903: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
904: j = j + user->mx*user->mx;
905: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
906: }
907: }
909: MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY);
910: MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY);
912: MatCreate(PETSC_COMM_WORLD,&user->L);
913: MatSetSizes(user->L,PETSC_DECIDE,ysubnlocal,m+n,n);
914: MatSetFromOptions(user->L);
915: MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL);
916: MatSeqAIJSetPreallocation(user->L,2,NULL);
917: MatGetOwnershipRange(user->L,&istart,&iend);
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: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
924: j = j+1;
925: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
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: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
931: j = j + user->mx;
932: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
933: }
934: if (i>=2*m/3 && i<m){
935: j = i-2*m/3;
936: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
937: j = j + user->mx*user->mx;
938: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
939: }
940: if (i>=m){
941: j = i - m;
942: MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES);
943: }
944: }
945: MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY);
946: MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY);
947: MatScale(user->L,PetscPowScalar(h,1.5));
949: /* Generate Div matrix */
950: if (!user->use_ptap) {
951: /* Generate Div matrix */
952: MatCreate(PETSC_COMM_WORLD,&user->Div);
953: MatSetSizes(user->Div,ysubnlocal,PETSC_DECIDE,n,m);
954: MatSetFromOptions(user->Div);
955: MatMPIAIJSetPreallocation(user->Div,4,NULL,4,NULL);
956: MatSeqAIJSetPreallocation(user->Div,6,NULL);
957: MatGetOwnershipRange(user->Grad,&istart,&iend);
959: for (i=istart; i<iend; i++){
960: if (i<m/3){
961: iblock = i / (user->mx-1);
962: j = iblock*user->mx + (i % (user->mx-1));
963: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
964: j = j+1;
965: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
966: }
967: if (i>=m/3 && i<2*m/3){
968: iblock = (i-m/3) / (user->mx*(user->mx-1));
969: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
970: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
971: j = j + user->mx;
972: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
973: }
974: if (i>=2*m/3){
975: j = i-2*m/3;
976: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
977: j = j + user->mx*user->mx;
978: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
979: }
980: }
982: MatAssemblyBegin(user->Div,MAT_FINAL_ASSEMBLY);
983: MatAssemblyEnd(user->Div,MAT_FINAL_ASSEMBLY);
984: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
985: } else {
986: MatCreate(PETSC_COMM_WORLD,&user->Diag);
987: MatSetSizes(user->Diag,PETSC_DECIDE,PETSC_DECIDE,m,m);
988: MatSetFromOptions(user->Diag);
989: MatMPIAIJSetPreallocation(user->Diag,1,NULL,0,NULL);
990: MatSeqAIJSetPreallocation(user->Diag,1,NULL);
991: }
993: /* Build work vectors and matrices */
994: VecCreate(PETSC_COMM_WORLD,&user->S);
995: VecSetSizes(user->S, PETSC_DECIDE, m);
996: VecSetFromOptions(user->S);
998: VecCreate(PETSC_COMM_WORLD,&user->lwork);
999: VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);
1000: VecSetFromOptions(user->lwork);
1002: MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork);
1004: VecDuplicate(user->S,&user->Swork);
1005: VecDuplicate(user->S,&user->Sdiag);
1006: VecDuplicate(user->S,&user->Av_u);
1007: VecDuplicate(user->S,&user->Twork);
1008: VecDuplicate(user->y,&user->ywork);
1009: VecDuplicate(user->u,&user->Ywork);
1010: VecDuplicate(user->u,&user->uwork);
1011: VecDuplicate(user->u,&user->js_diag);
1012: VecDuplicate(user->c,&user->cwork);
1013: VecDuplicate(user->d,&user->dwork);
1015: /* Create a matrix-free shell user->Jd for computing B*x */
1016: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal,user->nstate,user->ndesign,user,&user->Jd);
1017: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1018: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1020: /* Compute true state function ytrue given utrue */
1021: VecDuplicate(user->y,&user->ytrue);
1023: /* First compute Av_u = Av*exp(-u) */
1024: VecSet(user->uwork, 0);
1025: VecAXPY(user->uwork,-1.0,user->utrue); /* Note: user->utrue */
1026: VecExp(user->uwork);
1027: MatMult(user->Av,user->uwork,user->Av_u);
1029: /* Next form DSG = Div*S*Grad */
1030: VecCopy(user->Av_u,user->Swork);
1031: VecReciprocal(user->Swork);
1032: if (user->use_ptap) {
1033: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1034: MatPtAP(user->Diag,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);
1035: } else {
1036: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1037: MatDiagonalScale(user->Divwork,NULL,user->Swork);
1038: MatMatMultSymbolic(user->Divwork,user->Grad,1.0,&user->DSG);
1039: MatMatMultNumeric(user->Divwork,user->Grad,user->DSG);
1040: }
1042: MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1043: MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1045: if (user->use_lrc == PETSC_TRUE) {
1046: v=PetscSqrtReal(1.0 /user->ndesign);
1047: PetscMalloc1(user->ndesign,&user->ones);
1049: for (i=0;i<user->ndesign;i++) {
1050: user->ones[i]=v;
1051: }
1052: MatCreateDense(PETSC_COMM_WORLD,ysubnlocal,PETSC_DECIDE,user->ndesign,1,user->ones,&user->Ones);
1053: MatAssemblyBegin(user->Ones, MAT_FINAL_ASSEMBLY);
1054: MatAssemblyEnd(user->Ones, MAT_FINAL_ASSEMBLY);
1055: MatCreateLRC(user->DSG,user->Ones,NULL,user->Ones,&user->JsBlock);
1056: MatSetUp(user->JsBlock);
1057: } else {
1058: /* Create matrix-free shell user->Js for computing (A + h^3*e*e^T)*x */
1059: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal,ysubnlocal,user->ndesign,user->ndesign,user,&user->JsBlock);
1060: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateBlockMatMult);
1061: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateBlockMatMult);
1062: }
1063: MatSetOption(user->JsBlock,MAT_SYMMETRIC,PETSC_TRUE);
1064: MatSetOption(user->JsBlock,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1065: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->Js);
1066: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1067: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMult);
1068: MatSetOption(user->Js,MAT_SYMMETRIC,PETSC_TRUE);
1069: MatSetOption(user->Js,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1071: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->JsInv);
1072: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateInvMatMult);
1073: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateInvMatMult);
1074: MatSetOption(user->JsInv,MAT_SYMMETRIC,PETSC_TRUE);
1075: MatSetOption(user->JsInv,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1077: MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1078: MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1079: /* Now solve for ytrue */
1080: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1081: KSPSetFromOptions(user->solver);
1083: KSPSetOperators(user->solver,user->JsBlock,user->DSG);
1085: MatMult(user->JsInv,user->q,user->ytrue);
1086: /* First compute Av_u = Av*exp(-u) */
1087: VecSet(user->uwork,0);
1088: VecAXPY(user->uwork,-1.0,user->u); /* Note: user->u */
1089: VecExp(user->uwork);
1090: MatMult(user->Av,user->uwork,user->Av_u);
1092: /* Next update DSG = Div*S*Grad with user->u */
1093: VecCopy(user->Av_u,user->Swork);
1094: VecReciprocal(user->Swork);
1095: if (user->use_ptap) {
1096: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1097: MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
1098: } else {
1099: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1100: MatDiagonalScale(user->Divwork,NULL,user->Av_u);
1101: MatMatMultNumeric(user->Divwork,user->Grad,user->DSG);
1102: }
1104: /* Now solve for y */
1106: MatMult(user->JsInv,user->q,user->y);
1108: user->ksp_its_initial = user->ksp_its;
1109: user->ksp_its = 0;
1110: /* Construct projection matrix Q (blocks) */
1111: MatCreate(PETSC_COMM_WORLD,&user->Q);
1112: MatSetSizes(user->Q,dsubnlocal,ysubnlocal,user->ndata,user->ndesign);
1113: MatSetFromOptions(user->Q);
1114: MatMPIAIJSetPreallocation(user->Q,8,NULL,8,NULL);
1115: MatSeqAIJSetPreallocation(user->Q,8,NULL);
1117: for (i=0; i<user->mx; i++){
1118: x[i] = h*(i+0.5);
1119: y[i] = h*(i+0.5);
1120: z[i] = h*(i+0.5);
1121: }
1122: MatGetOwnershipRange(user->Q,&istart,&iend);
1124: nx = user->mx; ny = user->mx; nz = user->mx;
1125: for (i=istart; i<iend; i++){
1127: xri = xr[i];
1128: im = 0;
1129: xim = x[im];
1130: while (xri>xim && im<nx){
1131: im = im+1;
1132: xim = x[im];
1133: }
1134: indx1 = im-1;
1135: indx2 = im;
1136: dx1 = xri - x[indx1];
1137: dx2 = x[indx2] - xri;
1139: yri = yr[i];
1140: im = 0;
1141: yim = y[im];
1142: while (yri>yim && im<ny){
1143: im = im+1;
1144: yim = y[im];
1145: }
1146: indy1 = im-1;
1147: indy2 = im;
1148: dy1 = yri - y[indy1];
1149: dy2 = y[indy2] - yri;
1151: zri = zr[i];
1152: im = 0;
1153: zim = z[im];
1154: while (zri>zim && im<nz){
1155: im = im+1;
1156: zim = z[im];
1157: }
1158: indz1 = im-1;
1159: indz2 = im;
1160: dz1 = zri - z[indz1];
1161: dz2 = z[indz2] - zri;
1163: Dx = x[indx2] - x[indx1];
1164: Dy = y[indy2] - y[indy1];
1165: Dz = z[indz2] - z[indz1];
1167: j = indx1 + indy1*nx + indz1*nx*ny;
1168: v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1169: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1171: j = indx1 + indy1*nx + indz2*nx*ny;
1172: v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1173: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1175: j = indx1 + indy2*nx + indz1*nx*ny;
1176: v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1177: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1179: j = indx1 + indy2*nx + indz2*nx*ny;
1180: v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1181: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1183: j = indx2 + indy1*nx + indz1*nx*ny;
1184: v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1185: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1187: j = indx2 + indy1*nx + indz2*nx*ny;
1188: v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1189: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1191: j = indx2 + indy2*nx + indz1*nx*ny;
1192: v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1193: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1195: j = indx2 + indy2*nx + indz2*nx*ny;
1196: v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1197: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1198: }
1200: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1201: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1202: /* Create MQ (composed of blocks of Q */
1203: MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->MQ);
1204: MatShellSetOperation(user->MQ,MATOP_MULT,(void(*)(void))QMatMult);
1205: MatShellSetOperation(user->MQ,MATOP_MULT_TRANSPOSE,(void(*)(void))QMatMultTranspose);
1207: /* Add noise to the measurement data */
1208: VecSet(user->ywork,1.0);
1209: VecAYPX(user->ywork,user->noise,user->ytrue);
1210: MatMult(user->MQ,user->ywork,user->d);
1212: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1213: PetscFree(x);
1214: PetscFree(y);
1215: PetscFree(z);
1216: PetscLogStagePop();
1217: return(0);
1218: }
1220: PetscErrorCode EllipticDestroy(AppCtx *user)
1221: {
1223: PetscInt i;
1226: MatDestroy(&user->DSG);
1227: KSPDestroy(&user->solver);
1228: MatDestroy(&user->Q);
1229: MatDestroy(&user->MQ);
1230: if (!user->use_ptap) {
1231: MatDestroy(&user->Div);
1232: MatDestroy(&user->Divwork);
1233: } else {
1234: MatDestroy(&user->Diag);
1235: }
1236: if (user->use_lrc) {
1237: MatDestroy(&user->Ones);
1238: }
1240: MatDestroy(&user->Grad);
1241: MatDestroy(&user->Av);
1242: MatDestroy(&user->Avwork);
1243: MatDestroy(&user->L);
1244: MatDestroy(&user->Js);
1245: MatDestroy(&user->Jd);
1246: MatDestroy(&user->JsBlock);
1247: MatDestroy(&user->JsInv);
1249: VecDestroy(&user->x);
1250: VecDestroy(&user->u);
1251: VecDestroy(&user->uwork);
1252: VecDestroy(&user->utrue);
1253: VecDestroy(&user->y);
1254: VecDestroy(&user->ywork);
1255: VecDestroy(&user->ytrue);
1256: VecDestroy(&user->c);
1257: VecDestroy(&user->cwork);
1258: VecDestroy(&user->ur);
1259: VecDestroy(&user->q);
1260: VecDestroy(&user->d);
1261: VecDestroy(&user->dwork);
1262: VecDestroy(&user->lwork);
1263: VecDestroy(&user->S);
1264: VecDestroy(&user->Swork);
1265: VecDestroy(&user->Sdiag);
1266: VecDestroy(&user->Ywork);
1267: VecDestroy(&user->Twork);
1268: VecDestroy(&user->Av_u);
1269: VecDestroy(&user->js_diag);
1270: ISDestroy(&user->s_is);
1271: ISDestroy(&user->d_is);
1272: VecDestroy(&user->suby);
1273: VecDestroy(&user->subd);
1274: VecDestroy(&user->subq);
1275: VecScatterDestroy(&user->state_scatter);
1276: VecScatterDestroy(&user->design_scatter);
1277: for (i=0;i<user->ns;i++) {
1278: VecScatterDestroy(&user->yi_scatter[i]);
1279: VecScatterDestroy(&user->di_scatter[i]);
1280: }
1281: PetscFree(user->yi_scatter);
1282: PetscFree(user->di_scatter);
1283: if (user->use_lrc) {
1284: PetscFree(user->ones);
1285: MatDestroy(&user->Ones);
1286: }
1287: return(0);
1288: }
1290: PetscErrorCode EllipticMonitor(Tao tao, void *ptr)
1291: {
1293: Vec X;
1294: PetscReal unorm,ynorm;
1295: AppCtx *user = (AppCtx*)ptr;
1298: TaoGetSolutionVector(tao,&X);
1299: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1300: VecAXPY(user->ywork,-1.0,user->ytrue);
1301: VecAXPY(user->uwork,-1.0,user->utrue);
1302: VecNorm(user->uwork,NORM_2,&unorm);
1303: VecNorm(user->ywork,NORM_2,&ynorm);
1304: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1305: return(0);
1306: }