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