Actual source code: adapthist.c

  1: #include <petsc/private/tshistoryimpl.h>

  3: typedef struct {
  4:   TSHistory hist;
  5:   PetscBool bw;
  6: } TSAdapt_History;

  8: static PetscErrorCode TSAdaptChoose_History(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter)
  9: {
 10:   PetscErrorCode  ierr;
 11:   PetscInt        step;
 12:   TSAdapt_History *thadapt = (TSAdapt_History*)adapt->data;

 15:   if (!thadapt->hist) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ORDER,"Need to call TSAdaptHistorySetHistory() first");
 16:   TSGetStepNumber(ts,&step);
 17:   TSHistoryGetTimeStep(thadapt->hist,thadapt->bw,step+1,next_h);
 18:   *accept  = PETSC_TRUE;
 19:   *next_sc = 0;
 20:   *wlte    = -1;
 21:   *wltea   = -1;
 22:   *wlter   = -1;
 23:   return(0);
 24: }

 26: static PetscErrorCode TSAdaptReset_History(TSAdapt adapt)
 27: {
 28:   TSAdapt_History *thadapt = (TSAdapt_History*)adapt->data;
 29:   PetscErrorCode  ierr;

 32:   TSHistoryDestroy(&thadapt->hist);
 33:   return(0);
 34: }

 36: static PetscErrorCode TSAdaptDestroy_History(TSAdapt adapt)
 37: {

 41:   TSAdaptReset_History(adapt);
 42:   PetscFree(adapt->data);
 43:   return(0);
 44: }

 46: /* this is not public, as TSHistory is not a public object */
 47: PetscErrorCode TSAdaptHistorySetTSHistory(TSAdapt adapt, TSHistory hist, PetscBool backward)
 48: {
 49:   PetscReal      *hist_t;
 50:   PetscInt       n;
 51:   PetscBool      flg;

 57:   PetscObjectTypeCompare((PetscObject)adapt,TSADAPTHISTORY,&flg);
 58:   if (!flg) return(0);
 59:   TSHistoryGetHistory(hist,&n,(const PetscReal**)&hist_t,NULL,NULL);
 60:   TSAdaptHistorySetHistory(adapt,n,hist_t,backward);
 61:   return(0);
 62: }

 64: /*@
 65:    TSAdaptHistoryGetStep - Gets time and time step for a given step number in the history

 67:    Logically Collective on TSAdapt

 69:    Input Parameters:
 70: +  adapt    - the TSAdapt context
 71: -  step     - the step number

 73:    Output Parameters:
 74: +  t  - the time corresponding to the requested step (can be NULL)
 75: -  dt - the time step to be taken at the requested step (can be NULL)

 77:    Notes: The time history is internally copied, and the user can free the hist array. The user still needs to specify the termination of the solve via TSSetMaxSteps().

 79:    Level: advanced

 81: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptHistorySetTrajectory(), TSADAPTHISTORY
 82: @*/
 83: PetscErrorCode TSAdaptHistoryGetStep(TSAdapt adapt, PetscInt step, PetscReal *t, PetscReal *dt)
 84: {
 85:   TSAdapt_History *thadapt;
 86:   PetscBool       flg;
 87:   PetscErrorCode  ierr;

 92:   PetscObjectTypeCompare((PetscObject)adapt,TSADAPTHISTORY,&flg);
 93:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Not for type %s",((PetscObject)adapt)->type_name);
 94:   thadapt = (TSAdapt_History*)adapt->data;
 95:   TSHistoryGetTimeStep(thadapt->hist,thadapt->bw,step,dt);
 96:   TSHistoryGetTime(thadapt->hist,thadapt->bw,step,t);
 97:   return(0);
 98: }

100: /*@
101:    TSAdaptHistorySetHistory - Sets a time history in the adaptor

103:    Logically Collective on TSAdapt

105:    Input Parameters:
106: +  adapt    - the TSAdapt context
107: .  n        - size of the time history
108: .  hist     - the time history
109: -  backward - if the time history has to be followed backward

111:    Notes: The time history is internally copied, and the user can free the hist array. The user still needs to specify the termination of the solve via TSSetMaxSteps().

113:    Level: advanced

115: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptHistorySetTrajectory(), TSADAPTHISTORY
116: @*/
117: PetscErrorCode TSAdaptHistorySetHistory(TSAdapt adapt, PetscInt n, PetscReal hist[], PetscBool backward)
118: {
119:   TSAdapt_History *thadapt;
120:   PetscBool       flg;
121:   PetscErrorCode  ierr;

128:   PetscObjectTypeCompare((PetscObject)adapt,TSADAPTHISTORY,&flg);
129:   if (!flg) return(0);
130:   thadapt = (TSAdapt_History*)adapt->data;
131:   TSHistoryDestroy(&thadapt->hist);
132:   TSHistoryCreate(PetscObjectComm((PetscObject)adapt),&thadapt->hist);
133:   TSHistorySetHistory(thadapt->hist,n,hist,NULL,PETSC_FALSE);
134:   thadapt->bw = backward;
135:   return(0);
136: }

138: /*@
139:    TSAdaptHistorySetTrajectory - Sets a time history in the adaptor from a given TSTrajectory

141:    Logically Collective on TSAdapt

143:    Input Parameters:
144: +  adapt    - the TSAdapt context
145: .  tj       - the TSTrajectory context
146: -  backward - if the time history has to be followed backward

148:    Notes: The time history is internally copied, and the user can destroy the TSTrajectory if not needed. The user still needs to specify the termination of the solve via TSSetMaxSteps().

150:    Level: advanced

152: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptHistorySetHistory(), TSADAPTHISTORY
153: @*/
154: PetscErrorCode TSAdaptHistorySetTrajectory(TSAdapt adapt, TSTrajectory tj, PetscBool backward)
155: {
156:   PetscBool       flg;
157:   PetscErrorCode  ierr;

163:   PetscObjectTypeCompare((PetscObject)adapt,TSADAPTHISTORY,&flg);
164:   if (!flg) return(0);
165:   TSAdaptHistorySetTSHistory(adapt,tj->tsh,backward);
166:   return(0);
167: }


170: /*MC
171:    TSADAPTHISTORY - Time stepping controller that follows a given time history, used for Tangent Linear Model simulations

173:    Level: developer

175: .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptHistorySetHistory()
176: M*/
177: PETSC_EXTERN PetscErrorCode TSAdaptCreate_History(TSAdapt adapt)
178: {
179:   PetscErrorCode     ierr;
180:   TSAdapt_History *thadapt;

183:   PetscNew(&thadapt);
184:   adapt->matchstepfac[0] = PETSC_SMALL; /* prevent from accumulation errors */
185:   adapt->matchstepfac[1] = 0.0;         /* we will always match the final step, prevent TSAdaptChoose to mess with it */
186:   adapt->data            = thadapt;

188:   adapt->ops->choose  = TSAdaptChoose_History;
189:   adapt->ops->reset   = TSAdaptReset_History;
190:   adapt->ops->destroy = TSAdaptDestroy_History;
191:   return(0);
192: }