1 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 4 typedef struct { 5 PetscViewer viewer; 6 } TSTrajectory_Basic; 7 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; 14 15 PetscFunctionBegin; 16 ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);CHKERRQ(ierr); 17 ierr = PetscViewerFileSetName(tjbasic->viewer,filename);CHKERRQ(ierr); /* this triggers PetscViewer to be set up again */ 18 ierr = PetscViewerSetUp(tjbasic->viewer);CHKERRQ(ierr); 19 ierr = VecView(X,tjbasic->viewer);CHKERRQ(ierr); 20 ierr = PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL);CHKERRQ(ierr); 21 if (stepnum && !tj->solution_only) { 22 Vec *Y; 23 PetscReal tprev; 24 25 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 26 for (i=0; i<ns; i++) { 27 ierr = VecView(Y[i],tjbasic->viewer);CHKERRQ(ierr); 28 } 29 ierr = TSGetPrevTime(ts,&tprev);CHKERRQ(ierr); 30 ierr = PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL);CHKERRQ(ierr); 31 } 32 /* Tangent linear sensitivities needed by second-order adjoint */ 33 if (ts->forward_solve) { 34 Mat A,*S; 35 36 ierr = TSForwardGetSensitivities(ts,NULL,&A);CHKERRQ(ierr); 37 ierr = MatView(A,tjbasic->viewer);CHKERRQ(ierr); 38 if (stepnum) { 39 ierr = TSForwardGetStages(ts,&ns,&S);CHKERRQ(ierr); 40 for (i=0; i<ns; i++) { 41 ierr = MatView(S[i],tjbasic->viewer);CHKERRQ(ierr); 42 } 43 } 44 } 45 PetscFunctionReturn(0); 46 } 47 48 static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj) 49 { 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 ierr = PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");CHKERRQ(ierr); 54 ierr = PetscOptionsTail();CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 59 { 60 PetscViewer viewer; 61 char filename[PETSC_MAX_PATH_LEN]; 62 PetscErrorCode ierr; 63 Vec Sol; 64 PetscInt ns,i; 65 66 PetscFunctionBegin; 67 ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);CHKERRQ(ierr); 68 ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 69 ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr); 70 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr); 71 ierr = VecLoad(Sol,viewer);CHKERRQ(ierr); 72 ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr); 73 if (stepnum && !tj->solution_only) { 74 Vec *Y; 75 PetscReal timepre; 76 77 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 78 for (i=0; i<ns; i++) { 79 ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr); 80 } 81 ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr); 82 if (tj->adjoint_solve_mode) { 83 ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr); 84 } 85 } 86 /* Tangent linear sensitivities needed by second-order adjoint */ 87 if (ts->forward_solve) { 88 Mat A,*S; 89 90 ierr = TSForwardGetSensitivities(ts,NULL,&A);CHKERRQ(ierr); 91 ierr = MatLoad(A,viewer);CHKERRQ(ierr); 92 if (stepnum) { 93 ierr = TSForwardGetStages(ts,&ns,&S);CHKERRQ(ierr); 94 for (i=0; i<ns; i++) { 95 ierr = MatLoad(S[i],viewer);CHKERRQ(ierr); 96 } 97 } 98 } 99 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts) 104 { 105 MPI_Comm comm; 106 PetscMPIInt rank; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = PetscObjectGetComm((PetscObject)tj,&comm);CHKERRQ(ierr); 111 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 112 if (!rank) { 113 char *dir = tj->dirname; 114 PetscBool flg; 115 116 if (!dir) { 117 char dtempname[16] = "TS-data-XXXXXX"; 118 ierr = PetscMkdtemp(dtempname);CHKERRQ(ierr); 119 ierr = PetscStrallocpy(dtempname,&tj->dirname);CHKERRQ(ierr); 120 } else { 121 ierr = PetscTestDirectory(dir,'w',&flg);CHKERRQ(ierr); 122 if (!flg) { 123 ierr = PetscTestFile(dir,'r',&flg);CHKERRQ(ierr); 124 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir); 125 ierr = PetscMkdir(dir);CHKERRQ(ierr); 126 } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname); 127 } 128 } 129 ierr = PetscBarrier((PetscObject)tj);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj) 134 { 135 TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 ierr = PetscViewerDestroy(&tjbasic->viewer);CHKERRQ(ierr); 140 ierr = PetscFree(tjbasic);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 /*MC 145 TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file 146 147 Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed. 148 149 This version saves the solutions at all the stages 150 151 $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format 152 153 Level: intermediate 154 155 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile() 156 157 M*/ 158 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts) 159 { 160 TSTrajectory_Basic *tjbasic; 161 PetscErrorCode ierr; 162 163 PetscFunctionBegin; 164 ierr = PetscNew(&tjbasic);CHKERRQ(ierr); 165 166 ierr = PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);CHKERRQ(ierr); 167 ierr = PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 168 ierr = PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr); 169 ierr = PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 170 tj->data = tjbasic; 171 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 PetscFunctionReturn(0); 178 } 179