Actual source code: basicsymplectic.c

petsc-master 2019-05-24
Report Typos and Errors
  1: /*
  2:   Code for Timestepping with basic symplectic integrators for separable Hamiltonian systems
  3: */
  4:  #include <petsc/private/tsimpl.h>
  5:  #include <petscdm.h>

  7: static TSBasicSymplecticType TSBasicSymplecticDefault = TSBASICSYMPLECTICSIEULER;
  8: static PetscBool TSBasicSymplecticRegisterAllCalled;
  9: static PetscBool TSBasicSymplecticPackageInitialized;

 11: typedef struct _BasicSymplecticScheme *BasicSymplecticScheme;
 12: typedef struct _BasicSymplecticSchemeLink *BasicSymplecticSchemeLink;

 14: struct _BasicSymplecticScheme {
 15:   char      *name;
 16:   PetscInt  order;
 17:   PetscInt  s;       /* number of stages */
 18:   PetscReal *c,*d;
 19: };
 20: struct _BasicSymplecticSchemeLink {
 21:   struct _BasicSymplecticScheme sch;
 22:   BasicSymplecticSchemeLink     next;
 23: };
 24: static BasicSymplecticSchemeLink BasicSymplecticSchemeList;
 25: typedef struct {
 26:   TS                    subts_p,subts_q; /* sub TS contexts that holds the RHSFunction pointers */
 27:   IS                    is_p,is_q; /* IS sets for position and momentum respectively */
 28:   Vec                   update;    /* a nest work vector for generalized coordinates */
 29:   BasicSymplecticScheme scheme;
 30: } TS_BasicSymplectic;

 32: /*MC
 33:   TSBASICSYMPLECTICSIEULER - first order semi-implicit Euler method
 34:   Level: intermediate
 35: .seealso: TSBASICSYMPLECTIC
 36: M*/

 38: /*MC
 39:   TSBASICSYMPLECTICVELVERLET - second order Velocity Verlet method (leapfrog method with starting process and determing velocity and position at the same time)
 40: Level: intermediate
 41: .seealso: TSBASICSYMPLECTIC
 42: M*/

 44: /*@C
 45:   TSBasicSymplecticRegisterAll - Registers all of the basic symplectic integration methods in TSBasicSymplectic

 47:   Not Collective, but should be called by all processes which will need the schemes to be registered

 49:   Level: advanced

 51: .keywords: TS, TSBASICSYMPLECTIC, register, all

 53: .seealso:  TSBasicSymplecticRegisterDestroy()
 54: @*/
 55: PetscErrorCode TSBasicSymplecticRegisterAll(void)
 56: {

 60:   if (TSBasicSymplecticRegisterAllCalled) return(0);
 61:   TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
 62:   {
 63:     PetscReal c[1] = {1.0},d[1] = {1.0};
 64:     TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER,1,1,c,d);
 65:   }
 66:   {
 67:     PetscReal c[2] = {0,1.0},d[2] = {0.5,0.5};
 68:     TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET,2,2,c,d);
 69:   }
 70:   {
 71:     PetscReal c[3] = {1,-2.0/3.0,2.0/3.0},d[3] = {-1.0/24.0,3.0/4.0,7.0/24.0};
 72:     TSBasicSymplecticRegister(TSBASICSYMPLECTIC3,3,3,c,d);
 73:   }
 74:   {
 75: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
 76:     PetscReal c[4] = {1.0/2.0/(2.0-CUBE../../../../..OFTWO),(1.0-CUBE../../../../..OFTWO)/2.0/(2.0-CUBE../../../../..OFTWO),(1.0-CUBE../../../../..OFTWO)/2.0/(2.0-CUBE../../../../..OFTWO),1.0/2.0/(2.0-CUBE../../../../..OFTWO)},d[4] = {1.0/(2.0-CUBE../../../../..OFTWO),-CUBE../../../../..OFTWO/(2.0-CUBE../../../../..OFTWO),1.0/(2.0-CUBE../../../../..OFTWO),0};
 77:     TSBasicSymplecticRegister(TSBASICSYMPLECTIC4,4,4,c,d);
 78:   }
 79:   return(0);
 80: }

 82: /*@C
 83:    TSBasicSymplecticRegisterDestroy - Frees the list of schemes that were registered by TSBasicSymplecticRegister().

 85:    Not Collective

 87:    Level: advanced

 89: .keywords: TSBasicSymplectic, register, destroy
 90: .seealso: TSBasicSymplecticRegister(), TSBasicSymplecticRegisterAll()
 91: @*/
 92: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
 93: {
 94:   PetscErrorCode            ierr;
 95:   BasicSymplecticSchemeLink link;

 98:   while ((link = BasicSymplecticSchemeList)) {
 99:     BasicSymplecticScheme scheme = &link->sch;
100:     BasicSymplecticSchemeList = link->next;
101:     PetscFree2(scheme->c,scheme->d);
102:     PetscFree(scheme->name);
103:     PetscFree(link);
104:   }
105:   TSBasicSymplecticRegisterAllCalled = PETSC_FALSE;
106:   return(0);
107: }

109: /*@C
110:   TSBasicSymplecticInitializePackage - This function initializes everything in the TSBasicSymplectic package. It is called
111:   from TSInitializePackage().

113:   Level: developer

115: .keywords: TS, TSBasicSymplectic, initialize, package
116: .seealso: PetscInitialize()
117: @*/
118: PetscErrorCode TSBasicSymplecticInitializePackage(void)
119: {

123:   if (TSBasicSymplecticPackageInitialized) return(0);
124:   TSBasicSymplecticPackageInitialized = PETSC_TRUE;
125:   TSBasicSymplecticRegisterAll();
126:   PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
127:   return(0);
128: }

130: /*@C
131:   TSBasicSymplecticFinalizePackage - This function destroys everything in the TSBasicSymplectic package. It is
132:   called from PetscFinalize().

134:   Level: developer

136: .keywords: Petsc, destroy, package
137: .seealso: PetscFinalize()
138: @*/
139: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
140: {

144:   TSBasicSymplecticPackageInitialized = PETSC_FALSE;
145:   TSBasicSymplecticRegisterDestroy();
146:   return(0);
147: }

149: /*@C
150:    TSBasicSymplecticRegister - register a basic symplectic integration scheme by providing the coefficients.

152:    Not Collective, but the same schemes should be registered on all processes on which they will be used

154:    Input Parameters:
155: +  name - identifier for method
156: .  order - approximation order of method
157: .  s - number of stages, this is the dimension of the matrices below
158: .  c - coefficients for updating generalized position (dimension s)
159: -  d - coefficients for updating generalized momentum (dimension s)

161:    Notes:
162:    Several symplectic methods are provided, this function is only needed to create new methods.

164:    Level: advanced

166: .keywords: TS, register

168: .seealso: TSBasicSymplectic
169: @*/
170: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name,PetscInt order,PetscInt s,PetscReal c[],PetscReal d[])
171: {
172:   BasicSymplecticSchemeLink link;
173:   BasicSymplecticScheme     scheme;
174:   PetscErrorCode  ierr;


181:   TSBasicSymplecticInitializePackage();
182:   PetscCalloc1(1,&link);
183:   scheme = &link->sch;
184:   PetscStrallocpy(name,&scheme->name);
185:   scheme->order = order;
186:   scheme->s = s;
187:   PetscMalloc2(s,&scheme->c,s,&scheme->d);
188:   PetscMemcpy(scheme->c,c,s*sizeof(c[0]));
189:   PetscMemcpy(scheme->d,d,s*sizeof(d[0]));
190:   link->next = BasicSymplecticSchemeList;
191:   BasicSymplecticSchemeList = link;
192:   return(0);
193: }

195: /*
196: The simplified form of the equations are:

198: $ p_{i+1} = p_i + c_i*g(q_i)*h
199: $ q_{i+1} = q_i + d_i*f(p_{i+1},t_{i+1})*h

201: Several symplectic integrators are given below. An illustrative way to use them is to consider a particle with position q and velocity p.

203: To apply a timestep with values c_{1,2},d_{1,2} to the particle, carry out the following steps:

205: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
206: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
207: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
208: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2

210: */
211: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
212: {
213:   TS_BasicSymplectic    *bsymp = (TS_BasicSymplectic*)ts->data;
214:   BasicSymplecticScheme scheme = bsymp->scheme;
215:   Vec                   solution = ts->vec_sol,update = bsymp->update,q,p,q_update,p_update;
216:   IS                    is_q = bsymp->is_q,is_p = bsymp->is_p;
217:   TS                    subts_q = bsymp->subts_q,subts_p = bsymp->subts_p;
218:   PetscBool             stageok;
219:   PetscReal             next_time_step = ts->time_step;
220:   PetscInt              iter;
221:   PetscErrorCode        ierr;

224:   VecGetSubVector(solution,is_q,&q);
225:   VecGetSubVector(solution,is_p,&p);
226:   VecGetSubVector(update,is_q,&q_update);
227:   VecGetSubVector(update,is_p,&p_update);

229:   for (iter = 0;iter<scheme->s;iter++) {
230:     TSPreStage(ts,ts->ptime);
231:     /* update velocity p */
232:     if (scheme->c[iter]) {
233:       TSComputeRHSFunction(subts_p,ts->ptime,q,p_update);
234:       VecAXPY(p,scheme->c[iter]*ts->time_step,p_update);
235:     }
236:     /* update position q */
237:     if (scheme->d[iter]) {
238:       TSComputeRHSFunction(subts_q,ts->ptime,p,q_update);
239:       VecAXPY(q,scheme->d[iter]*ts->time_step,q_update);
240:       ts->ptime = ts->ptime+scheme->d[iter]*ts->time_step;
241:     }
242:     TSPostStage(ts,ts->ptime,0,&solution);
243:     TSAdaptCheckStage(ts->adapt,ts,ts->ptime,solution,&stageok);
244:     if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
245:     TSFunctionDomainError(ts,ts->ptime+ts->time_step,update,&stageok);
246:     if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
247:   }

249:   ts->time_step = next_time_step;
250:   VecRestoreSubVector(solution,is_q,&q);
251:   VecRestoreSubVector(solution,is_p,&p);
252:   VecRestoreSubVector(update,is_q,&q_update);
253:   VecRestoreSubVector(update,is_p,&p_update);
254:   return(0);
255: }

257: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine,DM coarse,void *ctx)
258: {
260:   return(0);
261: }

263: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
264: {
266:   return(0);
267: }

269: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm,DM subdm,void *ctx)
270: {
272:   return(0);
273: }

275: static PetscErrorCode DMSubDomainRestrictHook_BasicSymplectic(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
276: {

279:   return(0);
280: }

282: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
283: {
284:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
285:   DM                 dm;
286:   PetscErrorCode     ierr;

289:   TSRHSSplitGetIS(ts,"position",&bsymp->is_q);
290:   TSRHSSplitGetIS(ts,"momentum",&bsymp->is_p);
291:   if (!bsymp->is_q || !bsymp->is_p) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names positon and momentum respectively in order to use -ts_type basicsymplectic");
292:   TSRHSSplitGetSubTS(ts,"position",&bsymp->subts_q);
293:   TSRHSSplitGetSubTS(ts,"momentum",&bsymp->subts_p);
294:   if (!bsymp->subts_q || !bsymp->subts_p) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up the RHSFunctions for position and momentum using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS");

296:   VecDuplicate(ts->vec_sol,&bsymp->update);

298:   TSGetAdapt(ts,&ts->adapt);
299:   TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
300:   TSGetDM(ts,&dm);
301:   if (dm) {
302:     DMCoarsenHookAdd(dm,DMCoarsenHook_BasicSymplectic,DMRestrictHook_BasicSymplectic,ts);
303:     DMSubDomainHookAdd(dm,DMSubDomainHook_BasicSymplectic,DMSubDomainRestrictHook_BasicSymplectic,ts);
304:   }
305:   return(0);
306: }

308: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
309: {
310:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
311:   PetscErrorCode     ierr;

314:   VecDestroy(&bsymp->update);
315:   return(0);
316: }

318: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
319: {

323:   TSReset_BasicSymplectic(ts);
324:   PetscFree(ts->data);
325:   return(0);
326: }

328: static PetscErrorCode TSSetFromOptions_BasicSymplectic(PetscOptionItems *PetscOptionsObject,TS ts)
329: {
330:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
331:   PetscErrorCode     ierr;

334:   PetscOptionsHead(PetscOptionsObject,"Basic symplectic integrator options");
335:   {
336:     BasicSymplecticSchemeLink link;
337:     PetscInt                  count,choice;
338:     PetscBool                 flg;
339:     const char                **namelist;

341:     for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) ;
342:     PetscMalloc1(count,(char***)&namelist);
343:     for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) namelist[count] = link->sch.name;
344:     PetscOptionsEList("-ts_basicsymplectic_type","Family of basic symplectic integration method","TSBasicSymplecticSetType",(const char*const*)namelist,count,bsymp->scheme->name,&choice,&flg);
345:     if (flg) {TSBasicSymplecticSetType(ts,namelist[choice]);}
346:     PetscFree(namelist);
347:   }
348:   PetscOptionsTail();
349:   return(0);
350: }

352: static PetscErrorCode TSView_BasicSymplectic(TS ts,PetscViewer viewer)
353: {
355:   return(0);
356: }

358: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts,PetscReal t,Vec X)
359: {
360:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
361:   Vec                update = bsymp->update;
362:   PetscReal          alpha = (ts->ptime - t)/ts->time_step;
363:   PetscErrorCode     ierr;

366:   VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
367:   VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
368:   return(0);
369: }

371: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
372: {
374:   *yr = 1.0 + xr;
375:   *yi = xi;
376:   return(0);
377: }

379: /*@C
380:   TSBasicSymplecticSetType - Set the type of the basic symplectic method

382:   Logically Collective on TS

384:   Input Parameter:
385: +  ts - timestepping context
386: -  bsymptype - type of the symplectic scheme

388:   Options Database:
389: .  -ts_basicsymplectic_type <scheme>

391:   Notes:
392:   The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum". Each split is associated with an IS object and a sub-TS that is intended to store the user-provided RHS function.

394:   Level: intermediate
395: @*/
396: PetscErrorCode TSBasicSymplecticSetType(TS ts,TSBasicSymplecticType bsymptype)
397: {

402:   PetscTryMethod(ts,"TSBasicSymplecticSetType_C",(TS,TSBasicSymplecticType),(ts,bsymptype));
403:   return(0);
404: }

406: /*@C
407:   TSBasicSymplecticGetType - Get the type of the basic symplectic method

409:   Logically Collective on TS

411:   Input Parameter:
412: +  ts - timestepping context
413: -  bsymptype - type of the basic symplectic scheme

415:   Level: intermediate
416: @*/
417: PetscErrorCode TSBasicSymplecticGetType(TS ts,TSBasicSymplecticType *bsymptype)
418: {

423:   PetscUseMethod(ts,"TSBasicSymplecticGetType_C",(TS,TSBasicSymplecticType*),(ts,bsymptype));
424:   return(0);
425: }

427: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts,TSBasicSymplecticType bsymptype)
428: {
429:   TS_BasicSymplectic        *bsymp = (TS_BasicSymplectic*)ts->data;
430:   BasicSymplecticSchemeLink link;
431:   PetscBool                 match;
432:   PetscErrorCode            ierr;

435:   if (bsymp->scheme) {
436:     PetscStrcmp(bsymp->scheme->name,bsymptype,&match);
437:     if (match) return(0);
438:   }
439:   for (link = BasicSymplecticSchemeList; link; link=link->next) {
440:     PetscStrcmp(link->sch.name,bsymptype,&match);
441:     if (match) {
442:       bsymp->scheme = &link->sch;
443:       return(0);
444:     }
445:   }
446:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",bsymptype);
447:   return(0);
448: }

450: static PetscErrorCode  TSBasicSymplecticGetType_BasicSymplectic(TS ts,TSBasicSymplecticType *bsymptype)
451: {
452:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;

455:   *bsymptype = bsymp->scheme->name;
456:   return(0);
457: }

459: /*MC
460:   TSBasicSymplectic - ODE solver using basic symplectic integration schemes

462:   These methods are intened for separable Hamiltonian systems

464: $  qdot = dH(q,p,t)/dp
465: $  pdot = -dH(q,p,t)/dq

467:   where the Hamiltonian can be split into the sum of kinetic energy and potential energy

469: $  H(q,p,t) = T(p,t) + V(q,t).

471:   As a result, the system can be genearlly represented by

473: $  qdot = f(p,t) = dT(p,t)/dp
474: $  pdot = g(q,t) = -dV(q,t)/dq

476:   and solved iteratively with

478: $  q_new = q_old + d_i*h*f(p_old,t_old)
479: $  t_new = t_old + d_i*h
480: $  p_new = p_old + c_i*h*g(p_new,t_new)
481: $  i=0,1,...,n.

483:   The solution vector should contain both q and p, which correspond to (generalized) position and momentum respectively. Note that the momentum component could simply be velocity in some representations.
484:   The symplectic solver always expects a two-way splitting with the split names being "position" and "momentum". Each split is associated with an IS object and a sub-TS that is intended to store the user-provided RHS function.

486:   Reference: wikipedia (https://en.wikipedia.org/wiki/Symplectic_integrator)

488:   Level: beginner

490: .seealso:  TSCreate(), TSSetType(), TSRHSSplitSetIS(), TSRHSSplitSetRHSFunction()

492: M*/
493: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
494: {
495:   TS_BasicSymplectic *bsymp;
496:   PetscErrorCode     ierr;

499:   TSBasicSymplecticInitializePackage();
500:   PetscNewLog(ts,&bsymp);
501:   ts->data = (void*)bsymp;

503:   ts->ops->setup           = TSSetUp_BasicSymplectic;
504:   ts->ops->step            = TSStep_BasicSymplectic;
505:   ts->ops->reset           = TSReset_BasicSymplectic;
506:   ts->ops->destroy         = TSDestroy_BasicSymplectic;
507:   ts->ops->setfromoptions  = TSSetFromOptions_BasicSymplectic;
508:   ts->ops->view            = TSView_BasicSymplectic;
509:   ts->ops->interpolate     = TSInterpolate_BasicSymplectic;
510:   ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;

512:   PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticSetType_C",TSBasicSymplecticSetType_BasicSymplectic);
513:   PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticGetType_C",TSBasicSymplecticGetType_BasicSymplectic);

515:   TSBasicSymplecticSetType(ts,TSBasicSymplecticDefault);
516:   return(0);
517: }