1 2 static char help[] = "Solves 1D heat equation with FEM formulation.\n\ 3 Input arguments are\n\ 4 -useAlhs: solve Alhs*U' = (Arhs*U + g) \n\ 5 otherwise, solve U' = inv(Alhs)*(Arhs*U + g) \n\n"; 6 7 /*-------------------------------------------------------------------------- 8 Solves 1D heat equation U_t = U_xx with FEM formulation: 9 Alhs*U' = rhs (= Arhs*U + g) 10 We thank Chris Cox <clcox@clemson.edu> for contributing the original code 11 ----------------------------------------------------------------------------*/ 12 13 #include <petscksp.h> 14 #include <petscts.h> 15 16 /* special variable - max size of all arrays */ 17 #define num_z 10 18 19 /* 20 User-defined application context - contains data needed by the 21 application-provided call-back routines. 22 */ 23 typedef struct { 24 Mat Amat; /* left hand side matrix */ 25 Vec ksp_rhs,ksp_sol; /* working vectors for formulating inv(Alhs)*(Arhs*U+g) */ 26 int max_probsz; /* max size of the problem */ 27 PetscBool useAlhs; /* flag (1 indicates solving Alhs*U' = Arhs*U+g */ 28 int nz; /* total number of grid points */ 29 PetscInt m; /* total number of interio grid points */ 30 Vec solution; /* global exact ts solution vector */ 31 PetscScalar *z; /* array of grid points */ 32 PetscBool debug; /* flag (1 indicates activation of debugging printouts) */ 33 } AppCtx; 34 35 extern PetscScalar exact(PetscScalar,PetscReal); 36 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); 37 extern PetscErrorCode Petsc_KSPSolve(AppCtx*); 38 extern PetscScalar bspl(PetscScalar*,PetscScalar,PetscInt,PetscInt,PetscInt[][2],PetscInt); 39 extern PetscErrorCode femBg(PetscScalar[][3],PetscScalar*,PetscInt,PetscScalar*,PetscReal); 40 extern PetscErrorCode femA(AppCtx*,PetscInt,PetscScalar*); 41 extern PetscErrorCode rhs(AppCtx*,PetscScalar*, PetscInt, PetscScalar*,PetscReal); 42 extern PetscErrorCode RHSfunction(TS,PetscReal,Vec,Vec,void*); 43 44 int main(int argc,char **argv) 45 { 46 PetscInt i,m,nz,steps,max_steps,k,nphase=1; 47 PetscScalar zInitial,zFinal,val,*z; 48 PetscReal stepsz[4],T,ftime; 49 TS ts; 50 SNES snes; 51 Mat Jmat; 52 AppCtx appctx; /* user-defined application context */ 53 Vec init_sol; /* ts solution vector */ 54 PetscMPIInt size; 55 56 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 57 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 58 PetscCheck(size == 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only"); 59 60 /* initializations */ 61 zInitial = 0.0; 62 zFinal = 1.0; 63 nz = num_z; 64 m = nz-2; 65 appctx.nz = nz; 66 max_steps = (PetscInt)10000; 67 68 appctx.m = m; 69 appctx.max_probsz = nz; 70 appctx.debug = PETSC_FALSE; 71 appctx.useAlhs = PETSC_FALSE; 72 73 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"",""); 74 PetscCall(PetscOptionsName("-debug",NULL,NULL,&appctx.debug)); 75 PetscCall(PetscOptionsName("-useAlhs",NULL,NULL,&appctx.useAlhs)); 76 PetscCall(PetscOptionsRangeInt("-nphase",NULL,NULL,nphase,&nphase,NULL,1,3)); 77 PetscOptionsEnd(); 78 T = 0.014/nphase; 79 80 /* create vector to hold ts solution */ 81 /*-----------------------------------*/ 82 PetscCall(VecCreate(PETSC_COMM_WORLD, &init_sol)); 83 PetscCall(VecSetSizes(init_sol, PETSC_DECIDE, m)); 84 PetscCall(VecSetFromOptions(init_sol)); 85 86 /* create vector to hold true ts soln for comparison */ 87 PetscCall(VecDuplicate(init_sol, &appctx.solution)); 88 89 /* create LHS matrix Amat */ 90 /*------------------------*/ 91 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat)); 92 PetscCall(MatSetFromOptions(appctx.Amat)); 93 PetscCall(MatSetUp(appctx.Amat)); 94 /* set space grid points - interio points only! */ 95 PetscCall(PetscMalloc1(nz+1,&z)); 96 for (i=0; i<nz; i++) z[i]=(i)*((zFinal-zInitial)/(nz-1)); 97 appctx.z = z; 98 femA(&appctx,nz,z); 99 100 /* create the jacobian matrix */ 101 /*----------------------------*/ 102 PetscCall(MatCreate(PETSC_COMM_WORLD, &Jmat)); 103 PetscCall(MatSetSizes(Jmat,PETSC_DECIDE,PETSC_DECIDE,m,m)); 104 PetscCall(MatSetFromOptions(Jmat)); 105 PetscCall(MatSetUp(Jmat)); 106 107 /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */ 108 PetscCall(VecDuplicate(init_sol,&appctx.ksp_rhs)); 109 PetscCall(VecDuplicate(init_sol,&appctx.ksp_sol)); 110 111 /* set initial guess */ 112 /*-------------------*/ 113 for (i=0; i<nz-2; i++) { 114 val = exact(z[i+1], 0.0); 115 PetscCall(VecSetValue(init_sol,i,(PetscScalar)val,INSERT_VALUES)); 116 } 117 PetscCall(VecAssemblyBegin(init_sol)); 118 PetscCall(VecAssemblyEnd(init_sol)); 119 120 /*create a time-stepping context and set the problem type */ 121 /*--------------------------------------------------------*/ 122 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 123 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 124 125 /* set time-step method */ 126 PetscCall(TSSetType(ts,TSCN)); 127 128 /* Set optional user-defined monitoring routine */ 129 PetscCall(TSMonitorSet(ts,Monitor,&appctx,NULL)); 130 /* set the right hand side of U_t = RHSfunction(U,t) */ 131 PetscCall(TSSetRHSFunction(ts,NULL,(PetscErrorCode (*)(TS,PetscScalar,Vec,Vec,void*))RHSfunction,&appctx)); 132 133 if (appctx.useAlhs) { 134 /* set the left hand side matrix of Amat*U_t = rhs(U,t) */ 135 136 /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the 137 * Alhs matrix without making a copy. Either finite difference the entire thing or use analytic Jacobians in both 138 * places. 139 */ 140 PetscCall(TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,&appctx)); 141 PetscCall(TSSetIJacobian(ts,appctx.Amat,appctx.Amat,TSComputeIJacobianConstant,&appctx)); 142 } 143 144 /* use petsc to compute the jacobian by finite differences */ 145 PetscCall(TSGetSNES(ts,&snes)); 146 PetscCall(SNESSetJacobian(snes,Jmat,Jmat,SNESComputeJacobianDefault,NULL)); 147 148 /* get the command line options if there are any and set them */ 149 PetscCall(TSSetFromOptions(ts)); 150 151 #if defined(PETSC_HAVE_SUNDIALS2) 152 { 153 TSType type; 154 PetscBool sundialstype=PETSC_FALSE; 155 PetscCall(TSGetType(ts,&type)); 156 PetscCall(PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundialstype)); 157 PetscCheck(!sundialstype || !appctx.useAlhs,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use Alhs formulation for TSSUNDIALS type"); 158 } 159 #endif 160 /* Sets the initial solution */ 161 PetscCall(TSSetSolution(ts,init_sol)); 162 163 stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */ 164 ftime = 0.0; 165 for (k=0; k<nphase; k++) { 166 if (nphase > 1) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Phase %" PetscInt_FMT " initial time %g, stepsz %g, duration: %g\n",k,(double)ftime,(double)stepsz[k],(double)((k+1)*T))); 167 PetscCall(TSSetTime(ts,ftime)); 168 PetscCall(TSSetTimeStep(ts,stepsz[k])); 169 PetscCall(TSSetMaxSteps(ts,max_steps)); 170 PetscCall(TSSetMaxTime(ts,(k+1)*T)); 171 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 172 173 /* loop over time steps */ 174 /*----------------------*/ 175 PetscCall(TSSolve(ts,init_sol)); 176 PetscCall(TSGetSolveTime(ts,&ftime)); 177 PetscCall(TSGetStepNumber(ts,&steps)); 178 stepsz[k+1] = stepsz[k]*1.5; /* change step size for the next phase */ 179 } 180 181 /* free space */ 182 PetscCall(TSDestroy(&ts)); 183 PetscCall(MatDestroy(&appctx.Amat)); 184 PetscCall(MatDestroy(&Jmat)); 185 PetscCall(VecDestroy(&appctx.ksp_rhs)); 186 PetscCall(VecDestroy(&appctx.ksp_sol)); 187 PetscCall(VecDestroy(&init_sol)); 188 PetscCall(VecDestroy(&appctx.solution)); 189 PetscCall(PetscFree(z)); 190 191 PetscCall(PetscFinalize()); 192 return 0; 193 } 194 195 /*------------------------------------------------------------------------ 196 Set exact solution 197 u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t) 198 --------------------------------------------------------------------------*/ 199 PetscScalar exact(PetscScalar z,PetscReal t) 200 { 201 PetscScalar val, ex1, ex2; 202 203 ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t); 204 ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t); 205 val = PetscSinScalar(6*PETSC_PI*z)*ex1 + 3.*PetscSinScalar(2*PETSC_PI*z)*ex2; 206 return val; 207 } 208 209 /* 210 Monitor - User-provided routine to monitor the solution computed at 211 each timestep. This example plots the solution and computes the 212 error in two different norms. 213 214 Input Parameters: 215 ts - the timestep context 216 step - the count of the current step (with 0 meaning the 217 initial condition) 218 time - the current time 219 u - the solution at this timestep 220 ctx - the user-provided context for this monitoring routine. 221 In this case we use the application context which contains 222 information about the problem size, workspace and the exact 223 solution. 224 */ 225 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx) 226 { 227 AppCtx *appctx = (AppCtx*)ctx; 228 PetscInt i,m=appctx->m; 229 PetscReal norm_2,norm_max,h=1.0/(m+1); 230 PetscScalar *u_exact; 231 232 /* Compute the exact solution */ 233 PetscCall(VecGetArrayWrite(appctx->solution,&u_exact)); 234 for (i=0; i<m; i++) u_exact[i] = exact(appctx->z[i+1],time); 235 PetscCall(VecRestoreArrayWrite(appctx->solution,&u_exact)); 236 237 /* Print debugging information if desired */ 238 if (appctx->debug) { 239 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Computed solution vector at time %g\n",(double)time)); 240 PetscCall(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 241 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n")); 242 PetscCall(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 243 } 244 245 /* Compute the 2-norm and max-norm of the error */ 246 PetscCall(VecAXPY(appctx->solution,-1.0,u)); 247 PetscCall(VecNorm(appctx->solution,NORM_2,&norm_2)); 248 249 norm_2 = PetscSqrtReal(h)*norm_2; 250 PetscCall(VecNorm(appctx->solution,NORM_MAX,&norm_max)); 251 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Timestep %" PetscInt_FMT ": time = %g, 2-norm error = %6.4f, max norm error = %6.4f\n",step,(double)time,(double)norm_2,(double)norm_max)); 252 253 /* 254 Print debugging information if desired 255 */ 256 if (appctx->debug) { 257 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error vector\n")); 258 PetscCall(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 259 } 260 return 0; 261 } 262 263 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 264 Function to solve a linear system using KSP 265 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ 266 267 PetscErrorCode Petsc_KSPSolve(AppCtx *obj) 268 { 269 KSP ksp; 270 PC pc; 271 272 /*create the ksp context and set the operators,that is, associate the system matrix with it*/ 273 PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp)); 274 PetscCall(KSPSetOperators(ksp,obj->Amat,obj->Amat)); 275 276 /*get the preconditioner context, set its type and the tolerances*/ 277 PetscCall(KSPGetPC(ksp,&pc)); 278 PetscCall(PCSetType(pc,PCLU)); 279 PetscCall(KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT)); 280 281 /*get the command line options if there are any and set them*/ 282 PetscCall(KSPSetFromOptions(ksp)); 283 284 /*get the linear system (ksp) solve*/ 285 PetscCall(KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol)); 286 287 PetscCall(KSPDestroy(&ksp)); 288 return 0; 289 } 290 291 /*********************************************************************** 292 Function to return value of basis function or derivative of basis function. 293 *********************************************************************** 294 295 Arguments: 296 x = array of xpoints or nodal values 297 xx = point at which the basis function is to be 298 evaluated. 299 il = interval containing xx. 300 iq = indicates which of the two basis functions in 301 interval intrvl should be used 302 nll = array containing the endpoints of each interval. 303 id = If id ~= 2, the value of the basis function 304 is calculated; if id = 2, the value of the 305 derivative of the basis function is returned. 306 ***********************************************************************/ 307 308 PetscScalar bspl(PetscScalar *x, PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id) 309 { 310 PetscScalar x1,x2,bfcn; 311 PetscInt i1,i2,iq1,iq2; 312 313 /* Determine which basis function in interval intrvl is to be used in */ 314 iq1 = iq; 315 if (iq1==0) iq2 = 1; 316 else iq2 = 0; 317 318 /* Determine endpoint of the interval intrvl */ 319 i1=nll[il][iq1]; 320 i2=nll[il][iq2]; 321 322 /* Determine nodal values at the endpoints of the interval intrvl */ 323 x1=x[i1]; 324 x2=x[i2]; 325 326 /* Evaluate basis function */ 327 if (id == 2) bfcn=(1.0)/(x1-x2); 328 else bfcn=(xx-x2)/(x1-x2); 329 return bfcn; 330 } 331 332 /*--------------------------------------------------------- 333 Function called by rhs function to get B and g 334 ---------------------------------------------------------*/ 335 PetscErrorCode femBg(PetscScalar btri[][3],PetscScalar *f,PetscInt nz,PetscScalar *z, PetscReal t) 336 { 337 PetscInt i,j,jj,il,ip,ipp,ipq,iq,iquad,iqq; 338 PetscInt nli[num_z][2],indx[num_z]; 339 PetscScalar dd,dl,zip,zipq,zz,b_z,bb_z,bij; 340 PetscScalar zquad[num_z][3],dlen[num_z],qdwt[3]; 341 342 /* initializing everything - btri and f are initialized in rhs.c */ 343 for (i=0; i < nz; i++) { 344 nli[i][0] = 0; 345 nli[i][1] = 0; 346 indx[i] = 0; 347 zquad[i][0] = 0.0; 348 zquad[i][1] = 0.0; 349 zquad[i][2] = 0.0; 350 dlen[i] = 0.0; 351 } /*end for (i)*/ 352 353 /* quadrature weights */ 354 qdwt[0] = 1.0/6.0; 355 qdwt[1] = 4.0/6.0; 356 qdwt[2] = 1.0/6.0; 357 358 /* 1st and last nodes have Dirichlet boundary condition - 359 set indices there to -1 */ 360 361 for (i=0; i < nz-1; i++) indx[i] = i-1; 362 indx[nz-1] = -1; 363 364 ipq = 0; 365 for (il=0; il < nz-1; il++) { 366 ip = ipq; 367 ipq = ip+1; 368 zip = z[ip]; 369 zipq = z[ipq]; 370 dl = zipq-zip; 371 zquad[il][0] = zip; 372 zquad[il][1] = (0.5)*(zip+zipq); 373 zquad[il][2] = zipq; 374 dlen[il] = PetscAbsScalar(dl); 375 nli[il][0] = ip; 376 nli[il][1] = ipq; 377 } 378 379 for (il=0; il < nz-1; il++) { 380 for (iquad=0; iquad < 3; iquad++) { 381 dd = (dlen[il])*(qdwt[iquad]); 382 zz = zquad[il][iquad]; 383 384 for (iq=0; iq < 2; iq++) { 385 ip = nli[il][iq]; 386 b_z = bspl(z,zz,il,iq,nli,2); 387 i = indx[ip]; 388 389 if (i > -1) { 390 for (iqq=0; iqq < 2; iqq++) { 391 ipp = nli[il][iqq]; 392 bb_z = bspl(z,zz,il,iqq,nli,2); 393 j = indx[ipp]; 394 bij = -b_z*bb_z; 395 396 if (j > -1) { 397 jj = 1+j-i; 398 btri[i][jj] += bij*dd; 399 } else { 400 f[i] += bij*dd*exact(z[ipp], t); 401 /* f[i] += 0.0; */ 402 /* if (il==0 && j==-1) { */ 403 /* f[i] += bij*dd*exact(zz,t); */ 404 /* }*/ /*end if*/ 405 } /*end else*/ 406 } /*end for (iqq)*/ 407 } /*end if (i>0)*/ 408 } /*end for (iq)*/ 409 } /*end for (iquad)*/ 410 } /*end for (il)*/ 411 return 0; 412 } 413 414 PetscErrorCode femA(AppCtx *obj,PetscInt nz,PetscScalar *z) 415 { 416 PetscInt i,j,il,ip,ipp,ipq,iq,iquad,iqq; 417 PetscInt nli[num_z][2],indx[num_z]; 418 PetscScalar dd,dl,zip,zipq,zz,bb,bbb,aij; 419 PetscScalar rquad[num_z][3],dlen[num_z],qdwt[3],add_term; 420 421 /* initializing everything */ 422 for (i=0; i < nz; i++) { 423 nli[i][0] = 0; 424 nli[i][1] = 0; 425 indx[i] = 0; 426 rquad[i][0] = 0.0; 427 rquad[i][1] = 0.0; 428 rquad[i][2] = 0.0; 429 dlen[i] = 0.0; 430 } /*end for (i)*/ 431 432 /* quadrature weights */ 433 qdwt[0] = 1.0/6.0; 434 qdwt[1] = 4.0/6.0; 435 qdwt[2] = 1.0/6.0; 436 437 /* 1st and last nodes have Dirichlet boundary condition - 438 set indices there to -1 */ 439 440 for (i=0; i < nz-1; i++) indx[i]=i-1; 441 indx[nz-1]=-1; 442 443 ipq = 0; 444 445 for (il=0; il < nz-1; il++) { 446 ip = ipq; 447 ipq = ip+1; 448 zip = z[ip]; 449 zipq = z[ipq]; 450 dl = zipq-zip; 451 rquad[il][0] = zip; 452 rquad[il][1] = (0.5)*(zip+zipq); 453 rquad[il][2] = zipq; 454 dlen[il] = PetscAbsScalar(dl); 455 nli[il][0] = ip; 456 nli[il][1] = ipq; 457 } /*end for (il)*/ 458 459 for (il=0; il < nz-1; il++) { 460 for (iquad=0; iquad < 3; iquad++) { 461 dd = (dlen[il])*(qdwt[iquad]); 462 zz = rquad[il][iquad]; 463 464 for (iq=0; iq < 2; iq++) { 465 ip = nli[il][iq]; 466 bb = bspl(z,zz,il,iq,nli,1); 467 i = indx[ip]; 468 if (i > -1) { 469 for (iqq=0; iqq < 2; iqq++) { 470 ipp = nli[il][iqq]; 471 bbb = bspl(z,zz,il,iqq,nli,1); 472 j = indx[ipp]; 473 aij = bb*bbb; 474 if (j > -1) { 475 add_term = aij*dd; 476 PetscCall(MatSetValue(obj->Amat,i,j,add_term,ADD_VALUES)); 477 }/*endif*/ 478 } /*end for (iqq)*/ 479 } /*end if (i>0)*/ 480 } /*end for (iq)*/ 481 } /*end for (iquad)*/ 482 } /*end for (il)*/ 483 PetscCall(MatAssemblyBegin(obj->Amat,MAT_FINAL_ASSEMBLY)); 484 PetscCall(MatAssemblyEnd(obj->Amat,MAT_FINAL_ASSEMBLY)); 485 return 0; 486 } 487 488 /*--------------------------------------------------------- 489 Function to fill the rhs vector with 490 By + g values **** 491 ---------------------------------------------------------*/ 492 PetscErrorCode rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t) 493 { 494 PetscInt i,j,js,je,jj; 495 PetscScalar val,g[num_z],btri[num_z][3],add_term; 496 497 for (i=0; i < nz-2; i++) { 498 for (j=0; j <= 2; j++) btri[i][j]=0.0; 499 g[i] = 0.0; 500 } 501 502 /* call femBg to set the tri-diagonal b matrix and vector g */ 503 femBg(btri,g,nz,z,t); 504 505 /* setting the entries of the right hand side vector */ 506 for (i=0; i < nz-2; i++) { 507 val = 0.0; 508 js = 0; 509 if (i == 0) js = 1; 510 je = 2; 511 if (i == nz-2) je = 1; 512 513 for (jj=js; jj <= je; jj++) { 514 j = i+jj-1; 515 val += (btri[i][jj])*(y[j]); 516 } 517 add_term = val + g[i]; 518 PetscCall(VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES)); 519 } 520 PetscCall(VecAssemblyBegin(obj->ksp_rhs)); 521 PetscCall(VecAssemblyEnd(obj->ksp_rhs)); 522 return 0; 523 } 524 525 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 526 %% Function to form the right hand side of the time-stepping problem. %% 527 %% -------------------------------------------------------------------------------------------%% 528 if (useAlhs): 529 globalout = By+g 530 else if (!useAlhs): 531 globalout = f(y,t)=Ainv(By+g), 532 in which the ksp solver to transform the problem A*ydot=By+g 533 to the problem ydot=f(y,t)=inv(A)*(By+g) 534 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ 535 536 PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 537 { 538 AppCtx *obj = (AppCtx*)ctx; 539 PetscScalar soln[num_z]; 540 const PetscScalar *soln_ptr; 541 PetscInt i,nz=obj->nz; 542 PetscReal time; 543 544 /* get the previous solution to compute updated system */ 545 PetscCall(VecGetArrayRead(globalin,&soln_ptr)); 546 for (i=0; i < num_z-2; i++) soln[i] = soln_ptr[i]; 547 PetscCall(VecRestoreArrayRead(globalin,&soln_ptr)); 548 soln[num_z-1] = 0.0; 549 soln[num_z-2] = 0.0; 550 551 /* clear out the matrix and rhs for ksp to keep things straight */ 552 PetscCall(VecSet(obj->ksp_rhs,(PetscScalar)0.0)); 553 554 time = t; 555 /* get the updated system */ 556 rhs(obj,soln,nz,obj->z,time); /* setup of the By+g rhs */ 557 558 /* do a ksp solve to get the rhs for the ts problem */ 559 if (obj->useAlhs) { 560 /* ksp_sol = ksp_rhs */ 561 PetscCall(VecCopy(obj->ksp_rhs,globalout)); 562 } else { 563 /* ksp_sol = inv(Amat)*ksp_rhs */ 564 PetscCall(Petsc_KSPSolve(obj)); 565 PetscCall(VecCopy(obj->ksp_sol,globalout)); 566 } 567 return 0; 568 } 569 570 /*TEST 571 572 build: 573 requires: !complex 574 575 test: 576 suffix: euler 577 output_file: output/ex3.out 578 579 test: 580 suffix: 2 581 args: -useAlhs 582 output_file: output/ex3.out 583 TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant 584 585 TEST*/ 586