1 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscsys.h> 4 5 typedef struct _StackElement { 6 PetscInt stepnum; 7 Vec X; 8 PetscReal time; 9 Vec *Y; 10 PetscReal timeprev; 11 } *StackElement; 12 13 typedef struct _Stack { 14 PetscInt top; /* The top of the stack */ 15 PetscInt maxelements; /* The maximum stack size */ 16 PetscInt numY; 17 MPI_Comm comm; 18 StackElement *stack; /* The storage */ 19 } Stack; 20 21 static PetscErrorCode StackCreate(MPI_Comm,Stack *,PetscInt,PetscInt); 22 static PetscErrorCode StackDestroy(Stack*); 23 static PetscErrorCode StackPush(Stack*,StackElement); 24 static PetscErrorCode StackPop(Stack*,StackElement*); 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "StackCreate" 28 static PetscErrorCode StackCreate(MPI_Comm comm,Stack *s,PetscInt size,PetscInt ny) 29 { 30 PetscErrorCode ierr; 31 32 PetscFunctionBeginUser; 33 s->top = -1; 34 s->maxelements = size; 35 s->comm = comm; 36 s->numY = ny; 37 38 ierr = PetscMalloc1(s->maxelements*sizeof(StackElement),&s->stack);CHKERRQ(ierr); 39 ierr = PetscMemzero(s->stack,s->maxelements*sizeof(StackElement));CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "StackDestroy" 45 static PetscErrorCode StackDestroy(Stack *s) 46 { 47 PetscInt i; 48 PetscErrorCode ierr; 49 50 PetscFunctionBeginUser; 51 if (s->top>-1) { 52 for (i=0;i<=s->top;i++) { 53 ierr = VecDestroy(&s->stack[i]->X);CHKERRQ(ierr); 54 ierr = VecDestroyVecs(s->numY,&s->stack[i]->Y);CHKERRQ(ierr); 55 ierr = PetscFree(s->stack[i]);CHKERRQ(ierr); 56 } 57 } 58 ierr = PetscFree(s->stack);CHKERRQ(ierr); 59 ierr = PetscFree(s);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 #undef __FUNCT__ 64 #define __FUNCT__ "StackPush" 65 static PetscErrorCode StackPush(Stack *s,StackElement e) 66 { 67 PetscFunctionBeginUser; 68 if (s->top+1 >= s->maxelements) SETERRQ1(s->comm,PETSC_ERR_MEMC,"Maximum stack size (%D) exceeded",s->maxelements); 69 s->stack[++s->top] = e; 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "StackPop" 75 static PetscErrorCode StackPop(Stack *s,StackElement *e) 76 { 77 PetscFunctionBeginUser; 78 if (s->top == -1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Emptry stack"); 79 *e = s->stack[s->top--]; 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "TSTrajectorySet_Memory" 85 PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 86 { 87 PetscInt ns,i; 88 Vec *Y; 89 PetscReal timeprev; 90 StackElement e; 91 Stack *s = (Stack*)tj->data; 92 PetscErrorCode ierr; 93 94 PetscFunctionBeginUser; 95 ierr = PetscCalloc1(1,&e); 96 e->stepnum = stepnum; 97 ierr = VecDuplicate(X,&e->X);CHKERRQ(ierr); 98 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 99 e->time = time; 100 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 101 ierr = VecDuplicateVecs(Y[0],ns,&e->Y);CHKERRQ(ierr); 102 for (i=0;i<ns;i++) { 103 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 104 } 105 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 106 e->timeprev = timeprev; 107 ierr = StackPush(s,e);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "TSTrajectoryGet_Memory" 113 PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt step,PetscReal *t) 114 { 115 Vec Sol,*Y; 116 PetscInt nr,i; 117 StackElement e; 118 Stack *s = (Stack*)tj->data; 119 PetscErrorCode ierr; 120 121 PetscFunctionBeginUser; 122 ierr = StackPop(s,&e);CHKERRQ(ierr); 123 ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr); 124 ierr = VecCopy(e->X,Sol);CHKERRQ(ierr); 125 126 ierr = TSGetStages(ts,&nr,&Y);CHKERRQ(ierr); 127 for (i=0;i<nr ;i++) { 128 ierr = VecCopy(e->Y[i],Y[i]);CHKERRQ(ierr); 129 } 130 *t = e->time; 131 132 ierr = TSSetTimeStep(ts,-(*t)+e->timeprev);CHKERRQ(ierr); 133 134 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 135 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 136 ierr = PetscFree(e);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "TSTrajectoryDestroy_Memory" 142 PETSC_EXTERN PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 143 { 144 Stack *s = (Stack*)tj->data; 145 PetscErrorCode ierr; 146 147 PetscFunctionBeginUser; 148 ierr = StackDestroy(s);CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 /*MC 153 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 154 155 Level: intermediate 156 157 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 158 159 M*/ 160 #undef __FUNCT__ 161 #define __FUNCT__ "TSTrajectoryCreate_Memory" 162 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 163 { 164 PetscInt nr; 165 Stack *s; 166 PetscErrorCode ierr; 167 168 PetscFunctionBegin; 169 tj->ops->set = TSTrajectorySet_Memory; 170 tj->ops->get = TSTrajectoryGet_Memory; 171 tj->ops->destroy = TSTrajectoryDestroy_Memory; 172 173 ierr = PetscCalloc1(1,&s); 174 ierr = TSGetStages(ts,&nr,PETSC_IGNORE);CHKERRQ(ierr); 175 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,ts->max_steps,nr);CHKERRQ(ierr); 176 177 tj->data = s; 178 PetscFunctionReturn(0); 179 } 180