Actual source code: traj.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  2:  #include <petsc/private/tsimpl.h>

  4: PetscFunctionList TSTrajectoryList              = NULL;
  5: PetscBool         TSTrajectoryRegisterAllCalled = PETSC_FALSE;
  6: PetscClassId      TSTRAJECTORY_CLASSID;
  7: PetscLogEvent     TSTrajectory_Set, TSTrajectory_Get;

  9: /*@C
 10:   TSTrajectoryRegister - Adds a way of storing trajectories to the TS package

 12:   Not Collective

 14:   Input Parameters:
 15: + name        - the name of a new user-defined creation routine
 16: - create_func - the creation routine itself

 18:   Notes:
 19:   TSTrajectoryRegister() may be called multiple times to add several user-defined tses.

 21:   Level: developer

 23: .keywords: TS, trajectory, timestep, register

 25: .seealso: TSTrajectoryRegisterAll()
 26: @*/
 27: PetscErrorCode TSTrajectoryRegister(const char sname[],PetscErrorCode (*function)(TSTrajectory,TS))
 28: {

 32:   PetscFunctionListAdd(&TSTrajectoryList,sname,function);
 33:   return(0);
 34: }

 36: PetscErrorCode TSTrajectorySet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
 37: {

 41:   if (!tj) return(0);
 42:   PetscLogEventBegin(TSTrajectory_Set,tj,ts,0,0);
 43:   (*tj->ops->set)(tj,ts,stepnum,time,X);
 44:   PetscLogEventEnd(TSTrajectory_Set,tj,ts,0,0);
 45:   return(0);
 46: }

 48: PetscErrorCode TSTrajectoryGet(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *time)
 49: {

 53:   if (!tj) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TS solver did not save trajectory");
 54:   if (stepnum < 0) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_PLIB,"Requesting negative step number");
 55:   PetscLogEventBegin(TSTrajectory_Get,tj,ts,0,0);
 56:   (*tj->ops->get)(tj,ts,stepnum,time);
 57:   PetscLogEventEnd(TSTrajectory_Get,tj,ts,0,0);
 58:   return(0);
 59: }

 61: /*@C
 62:     TSTrajectoryView - Prints information about the trajectory object

 64:     Collective on TSTrajectory

 66:     Input Parameters:
 67: +   tj - the TSTrajectory context obtained from TSTrajectoryCreate()
 68: -   viewer - visualization context

 70:     Options Database Key:
 71: .   -ts_trajectory_view - calls TSTrajectoryView() at end of TSAdjointStep()

 73:     Notes:
 74:     The available visualization contexts include
 75: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 76: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 77:          output where only the first processor opens
 78:          the file.  All other processors send their
 79:          data to the first processor to print.

 81:     The user can open an alternative visualization context with
 82:     PetscViewerASCIIOpen() - output to a specified file.

 84:     Level: developer

 86: .keywords: TS, trajectory, timestep, view

 88: .seealso: PetscViewerASCIIOpen()
 89: @*/
 90: PetscErrorCode  TSTrajectoryView(TSTrajectory tj,PetscViewer viewer)
 91: {
 93:   PetscBool      iascii;

 97:   if (!viewer) {
 98:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj),&viewer);
 99:   }

103:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
104:   if (iascii) {
105:     PetscObjectPrintClassNamePrefixType((PetscObject)tj,viewer);
106:     PetscViewerASCIIPrintf(viewer,"  total number of recomputations for adjoint calculation = %D\n",tj->recomps);
107:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint reads = %D\n",tj->diskreads);
108:     PetscViewerASCIIPrintf(viewer,"  disk checkpoint writes = %D\n",tj->diskwrites);
109:     if (tj->ops->view) {
110:       PetscViewerASCIIPushTab(viewer);
111:       (*tj->ops->view)(tj,viewer);
112:       PetscViewerASCIIPopTab(viewer);
113:     }
114:   }
115:   return(0);
116: }

118: /*@C
119:    TSTrajectorySetVariableNames - Sets the name of each component in the solution vector so that it may be saved with the trajectory

121:    Collective on TSTrajectory

123:    Input Parameters:
124: +  tr - the trajectory context
125: -  names - the names of the components, final string must be NULL

127:    Level: intermediate

129:    Note: Fortran interface is not possible because of the string array argument

131: .keywords: TS, TSTrajectory, vector, monitor, view

133: .seealso: TSTrajectory, TSGetTrajectory()
134: @*/
135: PetscErrorCode  TSTrajectorySetVariableNames(TSTrajectory ctx,const char * const *names)
136: {
137:   PetscErrorCode    ierr;

140:   PetscStrArrayDestroy(&ctx->names);
141:   PetscStrArrayallocpy(names,&ctx->names);
142:   return(0);
143: }

145: /*@C
146:    TSTrjactorySetTransform - Solution vector will be transformed by provided function before being saved to disk

148:    Collective on TSLGCtx

150:    Input Parameters:
151: +  tj - the TSTrajectory context
152: .  transform - the transform function
153: .  destroy - function to destroy the optional context
154: -  ctx - optional context used by transform function

156:    Level: intermediate

158: .keywords: TSTrajectory,  vector, monitor, view

160: .seealso:  TSTrajectorySetVariableNames(), TSTrajectory, TSMonitorLGSetTransform()
161: @*/
162: PetscErrorCode  TSTrajectorySetTransform(TSTrajectory tj,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
163: {
165:   tj->transform        = transform;
166:   tj->transformdestroy = destroy;
167:   tj->transformctx     = tctx;
168:   return(0);
169: }


172: /*@
173:   TSTrajectoryCreate - This function creates an empty trajectory object used to store the time dependent solution of an ODE/DAE

175:   Collective on MPI_Comm

177:   Input Parameter:
178: . comm - the communicator

180:   Output Parameter:
181: . tj   - the trajectory object

183:   Level: developer

185:   Notes: Usually one does not call this routine, it is called automatically when one calls TSSetSaveTrajectory().

187: .keywords: TS, trajectory, create

189: .seealso: TSTrajectorySetUp(), TSTrajectoryDestroy(), TSTrajectorySetType(), TSTrajectorySetVariableNames(), TSGetTrajectory(), TSTrajectorySetKeepFiles()
190: @*/
191: PetscErrorCode  TSTrajectoryCreate(MPI_Comm comm,TSTrajectory *tj)
192: {
193:   TSTrajectory   t;

198:   *tj = NULL;
199:   TSInitializePackage();

201:   PetscHeaderCreate(t,TSTRAJECTORY_CLASSID,"TSTrajectory","Time stepping","TS",comm,TSTrajectoryDestroy,TSTrajectoryView);
202:   t->setupcalled = PETSC_FALSE;
203:   t->keepfiles   = PETSC_TRUE;
204:   *tj  = t;
205:   TSTrajectorySetDirname(t,"SA-data");
206:   TSTrajectorySetFiletemplate(t,"SA-%06D.bin");
207:   return(0);
208: }

210: /*@C
211:   TSTrajectorySetType - Sets the storage method to be used as in a trajectory

213:   Collective on TS

215:   Input Parameters:
216: + tj   - the TSTrajectory context
217: . ts   - the TS context
218: - type - a known method

220:   Options Database Command:
221: . -ts_trajectory_type <type> - Sets the method; use -help for a list of available methods (for instance, basic)

223:    Level: developer

225: .keywords: TS, trajectory, timestep, set, type

227: .seealso: TS, TSTrajectoryCreate(), TSTrajectorySetFromOptions(), TSTrajectoryDestroy()

229: @*/
230: PetscErrorCode  TSTrajectorySetType(TSTrajectory tj,TS ts,const TSTrajectoryType type)
231: {
232:   PetscErrorCode (*r)(TSTrajectory,TS);
233:   PetscBool      match;

238:   PetscObjectTypeCompare((PetscObject)tj,type,&match);
239:   if (match) return(0);

241:   PetscFunctionListFind(TSTrajectoryList,type,&r);
242:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSTrajectory type: %s",type);
243:   if (tj->ops->destroy) {
244:     (*(tj)->ops->destroy)(tj);

246:     tj->ops->destroy = NULL;
247:   }
248:   PetscMemzero(tj->ops,sizeof(*tj->ops));

250:   PetscObjectChangeTypeName((PetscObject)tj,type);
251:   (*r)(tj,ts);
252:   return(0);
253: }

255: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory,TS);
256: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory,TS);
257: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory,TS);
258: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory,TS);

260: /*@C
261:   TSTrajectoryRegisterAll - Registers all of the trajectory storage schecmes in the TS package.

263:   Not Collective

265:   Level: developer

267: .keywords: TS, trajectory, register, all

269: .seealso: TSTrajectoryRegister()
270: @*/
271: PetscErrorCode  TSTrajectoryRegisterAll(void)
272: {

276:   if (TSTrajectoryRegisterAllCalled) return(0);
277:   TSTrajectoryRegisterAllCalled = PETSC_TRUE;

279:   TSTrajectoryRegister(TSTRAJECTORYBASIC,TSTrajectoryCreate_Basic);
280:   TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE,TSTrajectoryCreate_Singlefile);
281:   TSTrajectoryRegister(TSTRAJECTORYMEMORY,TSTrajectoryCreate_Memory);
282:   TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION,TSTrajectoryCreate_Visualization);
283:   return(0);
284: }

286: /*@
287:    TSTrajectoryDestroy - Destroys a trajectory context

289:    Collective on TSTrajectory

291:    Input Parameter:
292: .  tj - the TSTrajectory context obtained from TSTrajectoryCreate()

294:    Level: developer

296: .keywords: TS, trajectory, timestep, destroy

298: .seealso: TSTrajectoryCreate(), TSTrajectorySetUp()
299: @*/
300: PetscErrorCode  TSTrajectoryDestroy(TSTrajectory *tj)
301: {

305:   if (!*tj) return(0);
307:   if (--((PetscObject)(*tj))->refct > 0) {*tj = 0; return(0);}

309:   if ((*tj)->transformdestroy) {(*(*tj)->transformdestroy)((*tj)->transformctx);}
310:   if ((*tj)->ops->destroy) {(*(*tj)->ops->destroy)((*tj));}
311:   PetscViewerDestroy(&(*tj)->monitor);
312:   PetscStrArrayDestroy(&(*tj)->names);
313:   PetscFree((*tj)->dirname);
314:   PetscFree((*tj)->filetemplate);
315:   PetscFree((*tj)->dirfiletemplate);
316:   PetscHeaderDestroy(tj);
317:   return(0);
318: }

320: /*
321:   TSTrajectorySetTypeFromOptions_Private - Sets the type of ts from user options.

323:   Collective on TSTrajectory

325:   Input Parameter:
326: + tj - the TSTrajectory context
327: - ts - the TS context

329:   Options Database Keys:
330: . -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION

332:   Level: developer

334: .keywords: TS, trajectory, set, options, type

336: .seealso: TSTrajectorySetFromOptions(), TSTrajectorySetType()
337: */
338: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,TSTrajectory tj,TS ts)
339: {
340:   PetscBool      opt;
341:   const char     *defaultType;
342:   char           typeName[256];
343:   PetscBool      flg;

347:   if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
348:   else defaultType = TSTRAJECTORYBASIC;

350:   TSTrajectoryRegisterAll();
351:   PetscOptionsFList("-ts_trajectory_type","TSTrajectory method","TSTrajectorySetType",TSTrajectoryList,defaultType,typeName,256,&opt);
352:   if (opt) {
353:     PetscStrcmp(typeName,TSTRAJECTORYMEMORY,&flg);
354:     TSTrajectorySetType(tj,ts,typeName);
355:   } else {
356:     TSTrajectorySetType(tj,ts,defaultType);
357:   }
358:   return(0);
359: }

361: /*@
362:    TSTrajectorySetMonitor - Monitor the schedules generated by the checkpointing controller

364:    Collective on TSTrajectory

366:    Input Arguments:
367: +  tj - the TSTrajectory context
368: -  flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable

370:    Options Database Keys:
371: .  -ts_trajectory_monitor - print TSTrajectory information

373:    Level: developer

375: .keywords: TS, trajectory, set, monitor

377: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp()
378: @*/
379: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj,PetscBool flg)
380: {

386:   if (flg) {
387:     if (!tj->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)tj),"stdout",&tj->monitor);}
388:   } else {
389:     PetscViewerDestroy(&tj->monitor);
390:   }
391:   return(0);
392: }

394: /*@
395:    TSTrajectorySetKeepFiles - Keep the files generated by the TSTrajectory

397:    Collective on TSTrajectory

399:    Input Arguments:
400: +  tj - the TSTrajectory context
401: -  flg - PETSC_TRUE to save, PETSC_FALSE to disable

403:    Options Database Keys:
404: .  -ts_trajectory_keep_files - have it keep the files

406:    Notes: By default the TSTrajectory used for adjoint computations, TSTRAJECTORYBASIC, removes the files it generates at the end of the run. This causes the files to be kept.

408:    Level: advanced

410: .keywords: TS, trajectory, set, monitor

412: .seealso: TSTrajectoryCreate(), TSTrajectoryDestroy(), TSTrajectorySetUp(), TSTrajectorySetMonitor()
413: @*/
414: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj,PetscBool flg)
415: {
419:   tj->keepfiles = flg;
420:   return(0);
421: }

423: /*@C
424:    TSTrajectorySetDirname - Specify the name of the directory where disk checkpoints are stored.

426:    Collective on TSTrajectory

428:    Input Arguments:
429: +  tj      - the TSTrajectory context
430: -  dirname - the directory name

432:    Options Database Keys:
433: .  -ts_trajectory_dirname - set the directory name

435:    Notes: The final location of the files is determined by dirname/filetemplate where filetemplate was provided by TSTrajectorySetFiletemplate()

437:    Level: developer

439: .keywords: TS, trajectory, set

441: .seealso: TSTrajectorySetFiletemplate(),TSTrajectorySetUp()
442: @*/
443: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj,const char dirname[])
444: {
448:   if (tj->dirfiletemplate) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set directoryname after it TSTrajectory has been setup");
449:   PetscFree(tj->dirname);
450:   PetscStrallocpy(dirname,&tj->dirname);
451:   return(0);
452: }

454: /*@C
455:    TSTrajectorySetFiletemplate - Specify the name template for the files storing checkpoints.

457:    Collective on TSTrajectory

459:    Input Arguments:
460: +  tj      - the TSTrajectory context
461: -  filetemplate - the template

463:    Options Database Keys:
464: .  -ts_trajectory_file_template - set the file name template

466:    Notes: The name template should be of the form, for example filename-%06D.bin It should not begin with a leading /

468:    The final location of the files is determined by dirname/filetemplate where dirname was provided by TSTrajectorySetDirname(). The %06D is replaced by the 
469:    timestep counter

471:    Level: developer

473: .keywords: TS, trajectory, set

475: .seealso: TSTrajectorySetDirname(),TSTrajectorySetUp()
476: @*/
477: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj,const char filetemplate[])
478: {
480:   const char     *ptr,*ptr2;

484:   if (tj->dirfiletemplate) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_ARG_WRONGSTATE,"Cannot set filetemplate after TSTrajectory has been setup");

486:   if (!filetemplate[0]) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
487:   /* Do some cursory validation of the input. */
488:   PetscStrstr(filetemplate,"%",(char**)&ptr);
489:   if (!ptr) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"-ts_trajectory_file_template requires a file name template, e.g. filename-%%06D.bin");
490:   for (ptr++; ptr && *ptr; ptr++) {
491:     PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
492:     if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_USER,"Invalid file template argument to -ts_trajectory_file_template, should look like filename-%%06D.bin");
493:     if (ptr2) break;
494:   }
495:   PetscFree(tj->filetemplate);
496:   PetscStrallocpy(filetemplate,&tj->filetemplate);
497:   return(0);
498: }

500: /*@
501:    TSTrajectorySetFromOptions - Sets various TSTrajectory parameters from user options.

503:    Collective on TSTrajectory

505:    Input Parameter:
506: +  tj - the TSTrajectory context obtained from TSTrajectoryCreate()
507: -  ts - the TS context

509:    Options Database Keys:
510: +  -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
511: .  -ts_trajectory_keep_files <true,false> - keep the files generated by the code after the program ends. This is true by default for TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION
512: -  -ts_trajectory_monitor - print TSTrajectory information

514:    Level: developer

516:    Notes: This is not normally called directly by users

518: .keywords: TS, trajectory, timestep, set, options, database

520: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
521: @*/
522: PetscErrorCode  TSTrajectorySetFromOptions(TSTrajectory tj,TS ts)
523: {
524:   PetscBool      set,flg;
525:   char           dirname[PETSC_MAX_PATH_LEN],filetemplate[PETSC_MAX_PATH_LEN];

531:   PetscObjectOptionsBegin((PetscObject)tj);
532:   TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject,tj,ts);
533:   PetscOptionsBool("-ts_trajectory_monitor","Print checkpointing schedules","TSTrajectorySetMonitor",tj->monitor ? PETSC_TRUE:PETSC_FALSE,&flg,&set);
534:   if (set) {TSTrajectorySetMonitor(tj,flg);}

536:   PetscOptionsBool("-ts_trajectory_keep_files","Keep any trajectory files generated during the run","TSTrajectorySetKeepFiles",tj->keepfiles,&flg,&set);
537:   if (set) {TSTrajectorySetKeepFiles(tj,flg);}

539:   PetscOptionsString("-ts_trajectory_dirname","Directory name for TSTrajectory file","TSTrajectorySetDirname",0,dirname,PETSC_MAX_PATH_LEN-14,&set);
540:   if (set) {
541:     TSTrajectorySetDirname(tj,dirname);
542:   }

544:   PetscOptionsString("-ts_trajectory_file_template","Template for TSTrajectory file name, use filename-%06D.bin","TSTrajectorySetFiletemplate",0,filetemplate,PETSC_MAX_PATH_LEN,&set);
545:   if (set) {
546:     TSTrajectorySetFiletemplate(tj,filetemplate);
547:   }

549:   /* Handle specific TSTrajectory options */
550:   if (tj->ops->setfromoptions) {
551:     (*tj->ops->setfromoptions)(PetscOptionsObject,tj);
552:   }
553:   PetscOptionsEnd();
554:   return(0);
555: }

557: /*@
558:    TSTrajectorySetUp - Sets up the internal data structures, e.g. stacks, for the later use
559:    of a TS trajectory.

561:    Collective on TS

563:    Input Parameter:
564: +  ts - the TS context obtained from TSCreate()
565: -  tj - the TS trajectory context

567:    Level: developer

569: .keywords: TS, trajectory, setup

571: .seealso: TSSetSaveTrajectory(), TSTrajectoryCreate(), TSTrajectoryDestroy()
572: @*/
573: PetscErrorCode  TSTrajectorySetUp(TSTrajectory tj,TS ts)
574: {
576:   size_t         s1,s2;

579:   if (!tj) return(0);
582:   if (tj->setupcalled) return(0);

584:   if (!((PetscObject)tj)->type_name) {
585:     TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
586:   }
587:   if (tj->ops->setup) {
588:     (*tj->ops->setup)(tj,ts);
589:   }

591:   tj->setupcalled = PETSC_TRUE;

593:   /* Set the counters to zero */
594:   tj->recomps    = 0;
595:   tj->diskreads  = 0;
596:   tj->diskwrites = 0;
597:   PetscStrlen(tj->dirname,&s1);
598:   PetscStrlen(tj->filetemplate,&s2);
599:   PetscMalloc((s1 + s2 + 10)*sizeof(char),&tj->dirfiletemplate);
600:   PetscSNPrintf(tj->dirfiletemplate,s1+s2+10,"%s/%s",tj->dirname,tj->filetemplate);
601:   return(0);
602: }