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