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