Actual source code: basicsymplectic.c

petsc-3.14.6 2021-03-30
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: .seealso:  TSBasicSymplecticRegisterDestroy()
 52: @*/
 53: PetscErrorCode TSBasicSymplecticRegisterAll(void)
 54: {

 58:   if (TSBasicSymplecticRegisterAllCalled) return(0);
 59:   TSBasicSymplecticRegisterAllCalled = PETSC_TRUE;
 60:   {
 61:     PetscReal c[1] = {1.0},d[1] = {1.0};
 62:     TSBasicSymplecticRegister(TSBASICSYMPLECTICSIEULER,1,1,c,d);
 63:   }
 64:   {
 65:     PetscReal c[2] = {0,1.0},d[2] = {0.5,0.5};
 66:     TSBasicSymplecticRegister(TSBASICSYMPLECTICVELVERLET,2,2,c,d);
 67:   }
 68:   {
 69:     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};
 70:     TSBasicSymplecticRegister(TSBASICSYMPLECTIC3,3,3,c,d);
 71:   }
 72:   {
 73: #define CUBE../../../../..OFTWO 1.2599210498948731647672106
 74:     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};
 75:     TSBasicSymplecticRegister(TSBASICSYMPLECTIC4,4,4,c,d);
 76:   }
 77:   return(0);
 78: }

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

 83:    Not Collective

 85:    Level: advanced

 87: .seealso: TSBasicSymplecticRegister(), TSBasicSymplecticRegisterAll()
 88: @*/
 89: PetscErrorCode TSBasicSymplecticRegisterDestroy(void)
 90: {
 91:   PetscErrorCode            ierr;
 92:   BasicSymplecticSchemeLink link;

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

106: /*@C
107:   TSBasicSymplecticInitializePackage - This function initializes everything in the TSBasicSymplectic package. It is called
108:   from TSInitializePackage().

110:   Level: developer

112: .seealso: PetscInitialize()
113: @*/
114: PetscErrorCode TSBasicSymplecticInitializePackage(void)
115: {

119:   if (TSBasicSymplecticPackageInitialized) return(0);
120:   TSBasicSymplecticPackageInitialized = PETSC_TRUE;
121:   TSBasicSymplecticRegisterAll();
122:   PetscRegisterFinalize(TSBasicSymplecticFinalizePackage);
123:   return(0);
124: }

126: /*@C
127:   TSBasicSymplecticFinalizePackage - This function destroys everything in the TSBasicSymplectic package. It is
128:   called from PetscFinalize().

130:   Level: developer

132: .seealso: PetscFinalize()
133: @*/
134: PetscErrorCode TSBasicSymplecticFinalizePackage(void)
135: {

139:   TSBasicSymplecticPackageInitialized = PETSC_FALSE;
140:   TSBasicSymplecticRegisterDestroy();
141:   return(0);
142: }

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

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

149:    Input Parameters:
150: +  name - identifier for method
151: .  order - approximation order of method
152: .  s - number of stages, this is the dimension of the matrices below
153: .  c - coefficients for updating generalized position (dimension s)
154: -  d - coefficients for updating generalized momentum (dimension s)

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

159:    Level: advanced

161: .seealso: TSBasicSymplectic
162: @*/
163: PetscErrorCode TSBasicSymplecticRegister(TSRosWType name,PetscInt order,PetscInt s,PetscReal c[],PetscReal d[])
164: {
165:   BasicSymplecticSchemeLink link;
166:   BasicSymplecticScheme     scheme;
167:   PetscErrorCode            ierr;


174:   TSBasicSymplecticInitializePackage();
175:   PetscNew(&link);
176:   scheme = &link->sch;
177:   PetscStrallocpy(name,&scheme->name);
178:   scheme->order = order;
179:   scheme->s = s;
180:   PetscMalloc2(s,&scheme->c,s,&scheme->d);
181:   PetscArraycpy(scheme->c,c,s);
182:   PetscArraycpy(scheme->d,d,s);
183:   link->next = BasicSymplecticSchemeList;
184:   BasicSymplecticSchemeList = link;
185:   return(0);
186: }

188: /*
189: The simplified form of the equations are:

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

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

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

198: - Update the velocity of the particle by adding to it its acceleration multiplied by c_1
199: - Update the position of the particle by adding to it its (updated) velocity multiplied by d_1
200: - Update the velocity of the particle by adding to it its acceleration (at the updated position) multiplied by c_2
201: - Update the position of the particle by adding to it its (double-updated) velocity multiplied by d_2

203: */
204: static PetscErrorCode TSStep_BasicSymplectic(TS ts)
205: {
206:   TS_BasicSymplectic    *bsymp = (TS_BasicSymplectic*)ts->data;
207:   BasicSymplecticScheme scheme = bsymp->scheme;
208:   Vec                   solution = ts->vec_sol,update = bsymp->update,q,p,q_update,p_update;
209:   IS                    is_q = bsymp->is_q,is_p = bsymp->is_p;
210:   TS                    subts_q = bsymp->subts_q,subts_p = bsymp->subts_p;
211:   PetscBool             stageok;
212:   PetscReal             next_time_step = ts->time_step;
213:   PetscInt              iter;
214:   PetscErrorCode        ierr;

217:   VecGetSubVector(solution,is_q,&q);
218:   VecGetSubVector(solution,is_p,&p);
219:   VecGetSubVector(update,is_q,&q_update);
220:   VecGetSubVector(update,is_p,&p_update);

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

242:   ts->time_step = next_time_step;
243:   VecRestoreSubVector(solution,is_q,&q);
244:   VecRestoreSubVector(solution,is_p,&p);
245:   VecRestoreSubVector(update,is_q,&q_update);
246:   VecRestoreSubVector(update,is_p,&p_update);
247:   return(0);
248: }

250: static PetscErrorCode DMCoarsenHook_BasicSymplectic(DM fine,DM coarse,void *ctx)
251: {
253:   return(0);
254: }

256: static PetscErrorCode DMRestrictHook_BasicSymplectic(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
257: {
259:   return(0);
260: }

262: static PetscErrorCode DMSubDomainHook_BasicSymplectic(DM dm,DM subdm,void *ctx)
263: {
265:   return(0);
266: }

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

272:   return(0);
273: }

275: static PetscErrorCode TSSetUp_BasicSymplectic(TS ts)
276: {
277:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
278:   DM                 dm;
279:   PetscErrorCode     ierr;

282:   TSRHSSplitGetIS(ts,"position",&bsymp->is_q);
283:   TSRHSSplitGetIS(ts,"momentum",&bsymp->is_p);
284:   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");
285:   TSRHSSplitGetSubTS(ts,"position",&bsymp->subts_q);
286:   TSRHSSplitGetSubTS(ts,"momentum",&bsymp->subts_p);
287:   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");

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

291:   TSGetAdapt(ts,&ts->adapt);
292:   TSAdaptCandidatesClear(ts->adapt); /* make sure to use fixed time stepping */
293:   TSGetDM(ts,&dm);
294:   if (dm) {
295:     DMCoarsenHookAdd(dm,DMCoarsenHook_BasicSymplectic,DMRestrictHook_BasicSymplectic,ts);
296:     DMSubDomainHookAdd(dm,DMSubDomainHook_BasicSymplectic,DMSubDomainRestrictHook_BasicSymplectic,ts);
297:   }
298:   return(0);
299: }

301: static PetscErrorCode TSReset_BasicSymplectic(TS ts)
302: {
303:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
304:   PetscErrorCode     ierr;

307:   VecDestroy(&bsymp->update);
308:   return(0);
309: }

311: static PetscErrorCode TSDestroy_BasicSymplectic(TS ts)
312: {

316:   TSReset_BasicSymplectic(ts);
317:   PetscFree(ts->data);
318:   return(0);
319: }

321: static PetscErrorCode TSSetFromOptions_BasicSymplectic(PetscOptionItems *PetscOptionsObject,TS ts)
322: {
323:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
324:   PetscErrorCode     ierr;

327:   PetscOptionsHead(PetscOptionsObject,"Basic symplectic integrator options");
328:   {
329:     BasicSymplecticSchemeLink link;
330:     PetscInt                  count,choice;
331:     PetscBool                 flg;
332:     const char                **namelist;

334:     for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) ;
335:     PetscMalloc1(count,(char***)&namelist);
336:     for (link=BasicSymplecticSchemeList,count=0; link; link=link->next,count++) namelist[count] = link->sch.name;
337:     PetscOptionsEList("-ts_basicsymplectic_type","Family of basic symplectic integration method","TSBasicSymplecticSetType",(const char*const*)namelist,count,bsymp->scheme->name,&choice,&flg);
338:     if (flg) {TSBasicSymplecticSetType(ts,namelist[choice]);}
339:     PetscFree(namelist);
340:   }
341:   PetscOptionsTail();
342:   return(0);
343: }

345: static PetscErrorCode TSView_BasicSymplectic(TS ts,PetscViewer viewer)
346: {
348:   return(0);
349: }

351: static PetscErrorCode TSInterpolate_BasicSymplectic(TS ts,PetscReal t,Vec X)
352: {
353:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;
354:   Vec                update = bsymp->update;
355:   PetscReal          alpha = (ts->ptime - t)/ts->time_step;
356:   PetscErrorCode     ierr;

359:   VecWAXPY(X,-ts->time_step,update,ts->vec_sol);
360:   VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);
361:   return(0);
362: }

364: static PetscErrorCode TSComputeLinearStability_BasicSymplectic(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
365: {
367:   *yr = 1.0 + xr;
368:   *yi = xi;
369:   return(0);
370: }

372: /*@C
373:   TSBasicSymplecticSetType - Set the type of the basic symplectic method

375:   Logically Collective on TS

377:   Input Parameter:
378: +  ts - timestepping context
379: -  bsymptype - type of the symplectic scheme

381:   Options Database:
382: .  -ts_basicsymplectic_type <scheme>

384:   Notes:
385:   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.

387:   Level: intermediate
388: @*/
389: PetscErrorCode TSBasicSymplecticSetType(TS ts,TSBasicSymplecticType bsymptype)
390: {

395:   PetscTryMethod(ts,"TSBasicSymplecticSetType_C",(TS,TSBasicSymplecticType),(ts,bsymptype));
396:   return(0);
397: }

399: /*@C
400:   TSBasicSymplecticGetType - Get the type of the basic symplectic method

402:   Logically Collective on TS

404:   Input Parameter:
405: +  ts - timestepping context
406: -  bsymptype - type of the basic symplectic scheme

408:   Level: intermediate
409: @*/
410: PetscErrorCode TSBasicSymplecticGetType(TS ts,TSBasicSymplecticType *bsymptype)
411: {

416:   PetscUseMethod(ts,"TSBasicSymplecticGetType_C",(TS,TSBasicSymplecticType*),(ts,bsymptype));
417:   return(0);
418: }

420: static PetscErrorCode TSBasicSymplecticSetType_BasicSymplectic(TS ts,TSBasicSymplecticType bsymptype)
421: {
422:   TS_BasicSymplectic        *bsymp = (TS_BasicSymplectic*)ts->data;
423:   BasicSymplecticSchemeLink link;
424:   PetscBool                 match;
425:   PetscErrorCode            ierr;

428:   if (bsymp->scheme) {
429:     PetscStrcmp(bsymp->scheme->name,bsymptype,&match);
430:     if (match) return(0);
431:   }
432:   for (link = BasicSymplecticSchemeList; link; link=link->next) {
433:     PetscStrcmp(link->sch.name,bsymptype,&match);
434:     if (match) {
435:       bsymp->scheme = &link->sch;
436:       return(0);
437:     }
438:   }
439:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",bsymptype);
440: }

442: static PetscErrorCode  TSBasicSymplecticGetType_BasicSymplectic(TS ts,TSBasicSymplecticType *bsymptype)
443: {
444:   TS_BasicSymplectic *bsymp = (TS_BasicSymplectic*)ts->data;

447:   *bsymptype = bsymp->scheme->name;
448:   return(0);
449: }

451: /*MC
452:   TSBasicSymplectic - ODE solver using basic symplectic integration schemes

454:   These methods are intened for separable Hamiltonian systems

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

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

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

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

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

468:   and solved iteratively with

470: $  q_new = q_old + d_i*h*f(p_old,t_old)
471: $  t_new = t_old + d_i*h
472: $  p_new = p_old + c_i*h*g(p_new,t_new)
473: $  i=0,1,...,n.

475:   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.
476:   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.

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

480:   Level: beginner

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

484: M*/
485: PETSC_EXTERN PetscErrorCode TSCreate_BasicSymplectic(TS ts)
486: {
487:   TS_BasicSymplectic *bsymp;
488:   PetscErrorCode     ierr;

491:   TSBasicSymplecticInitializePackage();
492:   PetscNewLog(ts,&bsymp);
493:   ts->data = (void*)bsymp;

495:   ts->ops->setup           = TSSetUp_BasicSymplectic;
496:   ts->ops->step            = TSStep_BasicSymplectic;
497:   ts->ops->reset           = TSReset_BasicSymplectic;
498:   ts->ops->destroy         = TSDestroy_BasicSymplectic;
499:   ts->ops->setfromoptions  = TSSetFromOptions_BasicSymplectic;
500:   ts->ops->view            = TSView_BasicSymplectic;
501:   ts->ops->interpolate     = TSInterpolate_BasicSymplectic;
502:   ts->ops->linearstability = TSComputeLinearStability_BasicSymplectic;

504:   PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticSetType_C",TSBasicSymplecticSetType_BasicSymplectic);
505:   PetscObjectComposeFunction((PetscObject)ts,"TSBasicSymplecticGetType_C",TSBasicSymplecticGetType_BasicSymplectic);

507:   TSBasicSymplecticSetType(ts,TSBasicSymplecticDefault);
508:   return(0);
509: }