1 /* 2 Code for timestepping with discrete gradient integrators 3 */ 4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 5 #include <petscdm.h> 6 7 PetscBool DGCite = PETSC_FALSE; 8 const char DGCitation[] = "@article{Gonzalez1996,\n" 9 " title = {Time integration and discrete Hamiltonian systems},\n" 10 " author = {Oscar Gonzalez},\n" 11 " journal = {Journal of Nonlinear Science},\n" 12 " volume = {6},\n" 13 " pages = {449--467},\n" 14 " doi = {10.1007/978-1-4612-1246-1_10},\n" 15 " year = {1996}\n}\n"; 16 17 typedef struct { 18 PetscReal stage_time; 19 Vec X0, X, Xdot; 20 void *funcCtx; 21 PetscBool gonzalez; 22 PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *); 23 PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *); 24 PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *); 25 } TS_DiscGrad; 26 27 static PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot) 28 { 29 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 if (X0) { 34 if (dm && dm != ts->dm) {ierr = DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0);CHKERRQ(ierr);} 35 else {*X0 = ts->vec_sol;} 36 } 37 if (Xdot) { 38 if (dm && dm != ts->dm) {ierr = DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);CHKERRQ(ierr);} 39 else {*Xdot = dg->Xdot;} 40 } 41 PetscFunctionReturn(0); 42 } 43 44 static PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot) 45 { 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 if (X0) { 50 if (dm && dm != ts->dm) {ierr = DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0);CHKERRQ(ierr);} 51 } 52 if (Xdot) { 53 if (dm && dm != ts->dm) {ierr = DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);CHKERRQ(ierr);} 54 } 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, void *ctx) 59 { 60 PetscFunctionBegin; 61 PetscFunctionReturn(0); 62 } 63 64 static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 65 { 66 TS ts = (TS) ctx; 67 Vec X0, Xdot, X0_c, Xdot_c; 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 ierr = TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot);CHKERRQ(ierr); 72 ierr = TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c);CHKERRQ(ierr); 73 ierr = MatRestrict(restrct, X0, X0_c);CHKERRQ(ierr); 74 ierr = MatRestrict(restrct, Xdot, Xdot_c);CHKERRQ(ierr); 75 ierr = VecPointwiseMult(X0_c, rscale, X0_c);CHKERRQ(ierr); 76 ierr = VecPointwiseMult(Xdot_c, rscale, Xdot_c);CHKERRQ(ierr); 77 ierr = TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot);CHKERRQ(ierr); 78 ierr = TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, void *ctx) 83 { 84 PetscFunctionBegin; 85 PetscFunctionReturn(0); 86 } 87 88 static PetscErrorCode DMSubDomainRestrictHook_TSDiscGrad(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 89 { 90 TS ts = (TS) ctx; 91 Vec X0, Xdot, X0_sub, Xdot_sub; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 ierr = TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr); 96 ierr = TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);CHKERRQ(ierr); 97 98 ierr = VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 99 ierr = VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 100 101 ierr = VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 102 ierr = VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 103 104 ierr = TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr); 105 ierr = TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 static PetscErrorCode TSSetUp_DiscGrad(TS ts) 110 { 111 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 112 DM dm; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 if (!dg->X) {ierr = VecDuplicate(ts->vec_sol, &dg->X);CHKERRQ(ierr);} 117 if (!dg->X0) {ierr = VecDuplicate(ts->vec_sol, &dg->X0);CHKERRQ(ierr);} 118 if (!dg->Xdot) {ierr = VecDuplicate(ts->vec_sol, &dg->Xdot);CHKERRQ(ierr);} 119 120 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 121 ierr = DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr); 122 ierr = DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode TSSetFromOptions_DiscGrad(PetscOptionItems *PetscOptionsObject, TS ts) 127 { 128 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 ierr = PetscOptionsHead(PetscOptionsObject, "Discrete Gradients ODE solver options");CHKERRQ(ierr); 133 { 134 ierr = PetscOptionsBool("-ts_discgrad_gonzalez","Use Gonzalez term in discrete gradients formulation","TSDiscGradUseGonzalez",dg->gonzalez,&dg->gonzalez,NULL);CHKERRQ(ierr); 135 } 136 ierr = PetscOptionsTail();CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 static PetscErrorCode TSView_DiscGrad(TS ts,PetscViewer viewer) 141 { 142 PetscBool iascii; 143 PetscErrorCode ierr; 144 145 PetscFunctionBegin; 146 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 147 if (iascii) { 148 ierr = PetscViewerASCIIPrintf(viewer," Discrete Gradients\n");CHKERRQ(ierr); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode TSDiscGradIsGonzalez_DiscGrad(TS ts,PetscBool *gonzalez) 154 { 155 TS_DiscGrad *dg = (TS_DiscGrad*)ts->data; 156 157 PetscFunctionBegin; 158 *gonzalez = dg->gonzalez; 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode TSDiscGradUseGonzalez_DiscGrad(TS ts,PetscBool flg) 163 { 164 TS_DiscGrad *dg = (TS_DiscGrad*)ts->data; 165 166 PetscFunctionBegin; 167 dg->gonzalez = flg; 168 PetscFunctionReturn(0); 169 } 170 171 static PetscErrorCode TSReset_DiscGrad(TS ts) 172 { 173 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 ierr = VecDestroy(&dg->X);CHKERRQ(ierr); 178 ierr = VecDestroy(&dg->X0);CHKERRQ(ierr); 179 ierr = VecDestroy(&dg->Xdot);CHKERRQ(ierr); 180 PetscFunctionReturn(0); 181 } 182 183 static PetscErrorCode TSDestroy_DiscGrad(TS ts) 184 { 185 DM dm; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = TSReset_DiscGrad(ts);CHKERRQ(ierr); 190 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 191 if (dm) { 192 ierr = DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr); 193 ierr = DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);CHKERRQ(ierr); 194 } 195 ierr = PetscFree(ts->data);CHKERRQ(ierr); 196 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",NULL);CHKERRQ(ierr); 197 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",NULL);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X) 202 { 203 TS_DiscGrad *dg = (TS_DiscGrad*)ts->data; 204 PetscReal dt = t - ts->ptime; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = VecCopy(ts->vec_sol, dg->X);CHKERRQ(ierr); 209 ierr = VecWAXPY(X, dt, dg->Xdot, dg->X);CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 213 static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x) 214 { 215 SNES snes; 216 PetscInt nits, lits; 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 221 ierr = SNESSolve(snes, b, x);CHKERRQ(ierr); 222 ierr = SNESGetIterationNumber(snes, &nits);CHKERRQ(ierr); 223 ierr = SNESGetLinearSolveIterations(snes, &lits);CHKERRQ(ierr); 224 ts->snes_its += nits; 225 ts->ksp_its += lits; 226 PetscFunctionReturn(0); 227 } 228 229 static PetscErrorCode TSStep_DiscGrad(TS ts) 230 { 231 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 232 TSAdapt adapt; 233 TSStepStatus status = TS_STEP_INCOMPLETE; 234 PetscInt rejections = 0; 235 PetscBool stageok, accept = PETSC_TRUE; 236 PetscReal next_time_step = ts->time_step; 237 PetscErrorCode ierr; 238 239 PetscFunctionBegin; 240 ierr = TSGetAdapt(ts, &adapt);CHKERRQ(ierr); 241 if (!ts->steprollback) {ierr = VecCopy(ts->vec_sol, dg->X0);CHKERRQ(ierr);} 242 243 while (!ts->reason && status != TS_STEP_COMPLETE) { 244 PetscReal shift = 1/(0.5*ts->time_step); 245 246 dg->stage_time = ts->ptime + 0.5*ts->time_step; 247 248 ierr = VecCopy(dg->X0, dg->X);CHKERRQ(ierr); 249 ierr = TSPreStage(ts, dg->stage_time);CHKERRQ(ierr); 250 ierr = TSDiscGrad_SNESSolve(ts, NULL, dg->X);CHKERRQ(ierr); 251 ierr = TSPostStage(ts, dg->stage_time, 0, &dg->X);CHKERRQ(ierr); 252 ierr = TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok);CHKERRQ(ierr); 253 if (!stageok) goto reject_step; 254 255 status = TS_STEP_PENDING; 256 ierr = VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X);CHKERRQ(ierr); 257 ierr = VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot);CHKERRQ(ierr); 258 ierr = TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept);CHKERRQ(ierr); 259 status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 260 if (!accept) { 261 ierr = VecCopy(dg->X0, ts->vec_sol);CHKERRQ(ierr); 262 ts->time_step = next_time_step; 263 goto reject_step; 264 } 265 ts->ptime += ts->time_step; 266 ts->time_step = next_time_step; 267 break; 268 269 reject_step: 270 ts->reject++; accept = PETSC_FALSE; 271 if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) { 272 ts->reason = TS_DIVERGED_STEP_REJECTED; 273 ierr = PetscInfo2(ts, "Step=%D, step rejections %D greater than current TS allowed, stopping solve\n", ts->steps, rejections);CHKERRQ(ierr); 274 } 275 } 276 PetscFunctionReturn(0); 277 } 278 279 static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y) 280 { 281 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 282 283 PetscFunctionBegin; 284 if (ns) *ns = 1; 285 if (Y) *Y = &(dg->X); 286 PetscFunctionReturn(0); 287 } 288 289 /* 290 This defines the nonlinear equation that is to be solved with SNES 291 G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0 292 */ 293 294 /* x = (x+x')/2 */ 295 /* NEED TO CALCULATE x_{n+1} from x and x_{n}*/ 296 static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts) 297 { 298 299 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 300 PetscReal norm, shift = 1/(0.5*ts->time_step); 301 PetscInt n; 302 Vec X0, Xdot, Xp, Xdiff; 303 Mat S; 304 PetscScalar F=0, F0=0, Gp; 305 Vec G, SgF; 306 DM dm, dmsave; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 311 312 ierr = VecDuplicate(y, &Xp);CHKERRQ(ierr); 313 ierr = VecDuplicate(y, &Xdiff);CHKERRQ(ierr); 314 ierr = VecDuplicate(y, &SgF);CHKERRQ(ierr); 315 ierr = VecDuplicate(y, &G);CHKERRQ(ierr); 316 317 ierr = VecGetLocalSize(y, &n);CHKERRQ(ierr); 318 ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr); 319 ierr = MatSetSizes(S, PETSC_DECIDE, PETSC_DECIDE, n, n);CHKERRQ(ierr); 320 ierr = MatSetFromOptions(S);CHKERRQ(ierr); 321 ierr = MatSetUp(S);CHKERRQ(ierr); 322 ierr = MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 323 ierr = MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 324 325 ierr = TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr); 326 ierr = VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x);CHKERRQ(ierr); /* Xdot = shift (x - X0) */ 327 328 ierr = VecAXPBYPCZ(Xp, -1, 2, 0, X0, x);CHKERRQ(ierr); /* Xp = 2*x - X0 + (0)*Xmid */ 329 ierr = VecAXPBYPCZ(Xdiff, -1, 1, 0, X0, Xp);CHKERRQ(ierr); /* Xdiff = xp - X0 + (0)*Xdiff */ 330 331 if (dg->gonzalez) { 332 ierr = (*dg->Sfunc)(ts, dg->stage_time, x , S, dg->funcCtx);CHKERRQ(ierr); 333 ierr = (*dg->Ffunc)(ts, dg->stage_time, Xp, &F, dg->funcCtx);CHKERRQ(ierr); 334 ierr = (*dg->Ffunc)(ts, dg->stage_time, X0, &F0, dg->funcCtx);CHKERRQ(ierr); 335 ierr = (*dg->Gfunc)(ts, dg->stage_time, x , G, dg->funcCtx);CHKERRQ(ierr); 336 337 /* Adding Extra Gonzalez Term */ 338 ierr = VecDot(Xdiff, G, &Gp);CHKERRQ(ierr); 339 ierr = VecNorm(Xdiff, NORM_2, &norm);CHKERRQ(ierr); 340 if (norm < PETSC_SQRT_MACHINE_EPSILON) { 341 Gp = 0; 342 } else { 343 /* Gp = (1/|xn+1 - xn|^2) * (F(xn+1) - F(xn) - Gp) */ 344 Gp = (F - F0 - Gp)/PetscSqr(norm); 345 } 346 ierr = VecAXPY(G, Gp, Xdiff);CHKERRQ(ierr); 347 ierr = MatMult(S, G , SgF);CHKERRQ(ierr); /* S*gradF */ 348 349 } else { 350 ierr = (*dg->Sfunc)(ts, dg->stage_time, x, S, dg->funcCtx);CHKERRQ(ierr); 351 ierr = (*dg->Gfunc)(ts, dg->stage_time, x, G, dg->funcCtx);CHKERRQ(ierr); 352 353 ierr = MatMult(S, G , SgF);CHKERRQ(ierr); /* Xdot = S*gradF */ 354 } 355 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 356 dmsave = ts->dm; 357 ts->dm = dm; 358 ierr = VecAXPBYPCZ(y, 1, -1, 0, Xdot, SgF);CHKERRQ(ierr); 359 ts->dm = dmsave; 360 ierr = TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);CHKERRQ(ierr); 361 362 ierr = VecDestroy(&Xp);CHKERRQ(ierr); 363 ierr = VecDestroy(&Xdiff);CHKERRQ(ierr); 364 ierr = VecDestroy(&SgF);CHKERRQ(ierr); 365 ierr = VecDestroy(&G);CHKERRQ(ierr); 366 ierr = MatDestroy(&S);CHKERRQ(ierr); 367 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts) 372 { 373 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 374 PetscReal shift = 1/(0.5*ts->time_step); 375 Vec Xdot; 376 DM dm,dmsave; 377 PetscErrorCode ierr; 378 379 PetscFunctionBegin; 380 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 381 /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */ 382 ierr = TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot);CHKERRQ(ierr); 383 384 dmsave = ts->dm; 385 ts->dm = dm; 386 ierr = TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE);CHKERRQ(ierr); 387 ts->dm = dmsave; 388 ierr = TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot);CHKERRQ(ierr); 389 PetscFunctionReturn(0); 390 } 391 392 static PetscErrorCode TSDiscGradGetFormulation_DiscGrad(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx) 393 { 394 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 395 396 PetscFunctionBegin; 397 *Sfunc = dg->Sfunc; 398 *Ffunc = dg->Ffunc; 399 *Gfunc = dg->Gfunc; 400 PetscFunctionReturn(0); 401 } 402 403 static PetscErrorCode TSDiscGradSetFormulation_DiscGrad(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx) 404 { 405 TS_DiscGrad *dg = (TS_DiscGrad *) ts->data; 406 407 PetscFunctionBegin; 408 dg->Sfunc = Sfunc; 409 dg->Ffunc = Ffunc; 410 dg->Gfunc = Gfunc; 411 dg->funcCtx = ctx; 412 PetscFunctionReturn(0); 413 } 414 415 /*MC 416 TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method 417 418 Notes: 419 This is the implicit midpoint rule, with an optional term that guarantees the discrete gradient property. This 420 timestepper applies to systems of the form 421 $ u_t = S(u) grad F(u) 422 where S(u) is a linear operator, and F is a functional of u. 423 424 .seealso: TSCreate(), TSSetType(), TS, TSDISCGRAD, TSDiscGradSetFormulation() 425 M*/ 426 PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts) 427 { 428 TS_DiscGrad *th; 429 PetscErrorCode ierr; 430 431 PetscFunctionBegin; 432 ierr = PetscCitationsRegister(DGCitation, &DGCite);CHKERRQ(ierr); 433 ts->ops->reset = TSReset_DiscGrad; 434 ts->ops->destroy = TSDestroy_DiscGrad; 435 ts->ops->view = TSView_DiscGrad; 436 ts->ops->setfromoptions = TSSetFromOptions_DiscGrad; 437 ts->ops->setup = TSSetUp_DiscGrad; 438 ts->ops->step = TSStep_DiscGrad; 439 ts->ops->interpolate = TSInterpolate_DiscGrad; 440 ts->ops->getstages = TSGetStages_DiscGrad; 441 ts->ops->snesfunction = SNESTSFormFunction_DiscGrad; 442 ts->ops->snesjacobian = SNESTSFormJacobian_DiscGrad; 443 ts->default_adapt_type = TSADAPTNONE; 444 445 ts->usessnes = PETSC_TRUE; 446 447 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 448 ts->data = (void*)th; 449 450 th->gonzalez = PETSC_FALSE; 451 452 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",TSDiscGradGetFormulation_DiscGrad);CHKERRQ(ierr); 453 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",TSDiscGradSetFormulation_DiscGrad);CHKERRQ(ierr); 454 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradIsGonzalez_C",TSDiscGradIsGonzalez_DiscGrad);CHKERRQ(ierr); 455 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradUseGonzalez_C",TSDiscGradUseGonzalez_DiscGrad);CHKERRQ(ierr); 456 PetscFunctionReturn(0); 457 } 458 459 /*@C 460 TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the formulation u_t = S grad F 461 462 Not Collective 463 464 Input Parameter: 465 . ts - timestepping context 466 467 Output Parameters: 468 + Sfunc - constructor for the S matrix from the formulation 469 . Ffunc - functional F from the formulation 470 - Gfunc - constructor for the gradient of F from the formulation 471 472 Calling sequence of Sfunc: 473 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *) 474 475 Calling sequence of Ffunc: 476 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *) 477 478 Calling sequence of Gfunc: 479 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *) 480 481 Level: intermediate 482 483 .seealso: TSDiscGradSetFormulation() 484 @*/ 485 PetscErrorCode TSDiscGradGetFormulation(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx) 486 { 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 491 PetscValidPointer(Sfunc, 2); 492 PetscValidPointer(Ffunc, 3); 493 PetscValidPointer(Gfunc, 4); 494 ierr = PetscUseMethod(ts,"TSDiscGradGetFormulation_C",(TS,PetscErrorCode(**Sfunc)(TS,PetscReal,Vec,Mat,void*),PetscErrorCode(**Ffunc)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode(**Gfunc)(TS,PetscReal,Vec,Vec,void*), void*),(ts,Sfunc,Ffunc,Gfunc,ctx));CHKERRQ(ierr); 495 PetscFunctionReturn(0); 496 } 497 498 /*@C 499 TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the formulation u_t = S(u) grad F(u) 500 501 Not Collective 502 503 Input Parameters: 504 + ts - timestepping context 505 . Sfunc - constructor for the S matrix from the formulation 506 . Ffunc - functional F from the formulation 507 - Gfunc - constructor for the gradient of F from the formulation 508 Calling sequence of Sfunc: 509 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *) 510 511 Calling sequence of Ffunc: 512 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *) 513 514 Calling sequence of Gfunc: 515 $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *) 516 517 Level: Intermediate 518 519 .seealso: TSDiscGradGetFormulation() 520 @*/ 521 PetscErrorCode TSDiscGradSetFormulation(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec , PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx) 522 { 523 PetscErrorCode ierr; 524 525 PetscFunctionBegin; 526 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 527 PetscValidFunction(Sfunc, 2); 528 PetscValidFunction(Ffunc, 3); 529 PetscValidFunction(Gfunc, 4); 530 ierr = PetscTryMethod(ts,"TSDiscGradSetFormulation_C",(TS,PetscErrorCode(*Sfunc)(TS,PetscReal,Vec,Mat,void*),PetscErrorCode(*Ffunc)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode(*Gfunc)(TS,PetscReal,Vec,Vec,void*), void*),(ts,Sfunc,Ffunc,Gfunc,ctx));CHKERRQ(ierr); 531 PetscFunctionReturn(0); 532 } 533 534 /*@ 535 TSDiscGradIsGonzalez - Checks flag for whether to use additional conservative terms in discrete gradient formulation. 536 537 Not Collective 538 539 Input Parameter: 540 . ts - timestepping context 541 542 Output Parameter: 543 . gonzalez - PETSC_TRUE when using the Gonzalez term 544 545 Level: Advanced 546 547 .seealso: TSDiscGradUseGonzalez(), TSDISCGRAD 548 @*/ 549 PetscErrorCode TSDiscGradIsGonzalez(TS ts,PetscBool *gonzalez) 550 { 551 PetscErrorCode ierr; 552 553 PetscFunctionBegin; 554 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 555 PetscValidPointer(gonzalez,2); 556 ierr = PetscUseMethod(ts,"TSDiscGradIsGonzalez_C",(TS,PetscBool*),(ts,gonzalez));CHKERRQ(ierr); 557 PetscFunctionReturn(0); 558 } 559 560 /*@ 561 TSDiscGradUseGonzalez - Sets discrete gradient formulation with or without additional conservative terms. Without flag, the discrete gradients timestepper is just backwards euler 562 563 Not Collective 564 565 Input Parameter: 566 + ts - timestepping context 567 - flg - PETSC_TRUE to use the Gonzalez term 568 569 Options Database: 570 . -ts_discgrad_gonzalez <flg> 571 572 Level: Intermediate 573 574 .seealso: TSDISCGRAD 575 @*/ 576 PetscErrorCode TSDiscGradUseGonzalez(TS ts,PetscBool flg) 577 { 578 PetscErrorCode ierr; 579 580 PetscFunctionBegin; 581 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 582 ierr = PetscTryMethod(ts,"TSDiscGradUseGonzalez_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 583 PetscFunctionReturn(0); 584 } 585