Actual source code: trajbasic.c

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

      4: typedef struct {
      5:   PetscViewer viewer;
      6: } TSTrajectory_Basic;

      8: static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
      9: {
     10:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
     11:   char               filename[PETSC_MAX_PATH_LEN];
     12:   PetscInt           ns,i;
     13:   PetscErrorCode     ierr;

     16:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
     17:   PetscViewerFileSetName(tjbasic->viewer,filename); /* this triggers PetscViewer to be set up again */
     18:   PetscViewerSetUp(tjbasic->viewer);
     19:   VecView(X,tjbasic->viewer);
     20:   PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL);
     21:   if (stepnum && !tj->solution_only) {
     22:     Vec       *Y;
     23:     PetscReal tprev;

     25:     TSGetStages(ts,&ns,&Y);
     26:     for (i=0; i<ns; i++) {
     27:       VecView(Y[i],tjbasic->viewer);
     28:     }
     29:     TSGetPrevTime(ts,&tprev);
     30:     PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL);
     31:   }
     32:   /* Tangent linear sensitivities needed by second-order adjoint */
     33:   if (ts->forward_solve) {
     34:     Mat A,*S;

     36:     TSForwardGetSensitivities(ts,NULL,&A);
     37:     MatView(A,tjbasic->viewer);
     38:     if (stepnum) {
     39:       TSForwardGetStages(ts,&ns,&S);
     40:       for (i=0; i<ns; i++) {
     41:         MatView(S[i],tjbasic->viewer);
     42:       }
     43:     }
     44:   }
     45:   return(0);
     46: }

     48: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
     49: {

     53:   PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");
     54:   PetscOptionsTail();
     55:   return(0);
     56: }

     58: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
     59: {
     60:   PetscViewer    viewer;
     61:   char           filename[PETSC_MAX_PATH_LEN];
     63:   Vec            Sol;
     64:   PetscInt       ns,i;

     67:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
     68:   PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);
     69:   TSGetSolution(ts,&Sol);
     70:   PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
     71:   VecLoad(Sol,viewer);
     72:   PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
     73:   if (stepnum && !tj->solution_only) {
     74:     Vec       *Y;
     75:     PetscReal timepre;

     77:     TSGetStages(ts,&ns,&Y);
     78:     for (i=0; i<ns; i++) {
     79:       VecLoad(Y[i],viewer);
     80:     }
     81:     PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
     82:     if (tj->adjoint_solve_mode) {
     83:       TSSetTimeStep(ts,-(*t)+timepre);
     84:     }
     85:   }
     86:   /* Tangent linear sensitivities needed by second-order adjoint */
     87:   if (ts->forward_solve) {
     88:     Mat A,*S;

     90:     TSForwardGetSensitivities(ts,NULL,&A);
     91:     MatLoad(A,viewer);
     92:     if (stepnum) {
     93:       TSForwardGetStages(ts,&ns,&S);
     94:       for (i=0; i<ns; i++) {
     95:         MatLoad(S[i],viewer);
     96:       }
     97:     }
     98:   }
     99:   PetscViewerDestroy(&viewer);
    100:   return(0);
    101: }

    103: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
    104: {
    105:   MPI_Comm       comm;
    106:   PetscMPIInt    rank;

    110:   PetscObjectGetComm((PetscObject)tj,&comm);
    111:   MPI_Comm_rank(comm,&rank);
    112:   if (!rank) {
    113:     char      *dir = tj->dirname;
    114:     PetscBool flg;

    116:     if (!dir) {
    117:       char dtempname[16] = "TS-data-XXXXXX";
    118:       PetscMkdtemp(dtempname);
    119:       PetscStrallocpy(dtempname,&tj->dirname);
    120:     } else {
    121:       PetscTestDirectory(dir,'w',&flg);
    122:       if (!flg) {
    123:         PetscTestFile(dir,'r',&flg);
    124:         if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
    125:         PetscMkdir(dir);
    126:       } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname);
    127:     }
    128:   }
    129:   PetscBarrier((PetscObject)tj);
    130:   return(0);
    131: }

    133: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
    134: {
    135:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
    136:   PetscErrorCode     ierr;

    139:   PetscViewerDestroy(&tjbasic->viewer);
    140:   PetscFree(tjbasic);
    141:   return(0);
    142: }

    144: /*MC
    145:       TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file

    147:       Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.

    149:       This version saves the solutions at all the stages

    151:       $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format

    153:   Level: intermediate

    155: .seealso:  TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()

    157: M*/
    158: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
    159: {
    160:   TSTrajectory_Basic *tjbasic;
    161:   PetscErrorCode     ierr;

    164:   PetscNew(&tjbasic);

    166:   PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
    167:   PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
    168:   PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
    169:   PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
    170:   tj->data = tjbasic;

    172:   tj->ops->set            = TSTrajectorySet_Basic;
    173:   tj->ops->get            = TSTrajectoryGet_Basic;
    174:   tj->ops->setup          = TSTrajectorySetUp_Basic;
    175:   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
    176:   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
    177:   return(0);
    178: }