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