Actual source code: adaptglee.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  1:  #include <petsc/private/tsimpl.h>

  3: typedef struct {
  4:   Vec E,Y;
  5: } TSAdapt_GLEE;

  7: static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter)
  8: {
  9:   TSAdapt_GLEE  *glee = (TSAdapt_GLEE*)adapt->data;
 10:   TSType         time_scheme;      /* Type of time-integration scheme        */
 12:   Vec            X,Y,E;
 13:   PetscReal      enorm,enorma,enormr,hfac_lte,hfac_ltea,hfac_lter,h_lte,safety;
 14:   PetscInt       order;
 15:   PetscBool      bGTEMethod=PETSC_FALSE;


 19:   *next_sc = 0; /* Reuse the same order scheme */
 20:   safety = adapt->safety;
 21:   TSGetType(ts,&time_scheme);
 22:   if (!strcmp(time_scheme,TSGLEE)) bGTEMethod=PETSC_TRUE;
 23:   order = adapt->candidates.order[0];

 25:   if (bGTEMethod){/* the method is of GLEE type */
 26:     TSGetSolution(ts,&X);
 27:     if (!glee->Y && adapt->glee_use_local) {
 28:       VecDuplicate(X,&glee->Y);/*create vector to store previous step global error*/
 29:       VecZeroEntries(glee->Y); /*set error to zero on the first step - may not work if error is not zero initially*/
 30:     }
 31:     if (!glee->E) {VecDuplicate(X,&glee->E);}
 32:     E    = glee->E;
 33:     TSGetTimeError(ts,0,&E);

 35:     if (adapt->glee_use_local) {VecAXPY(E,-1.0,glee->Y);} /* local error = current error - previous step error */

 37:     /* this should be called with the solution at the beginning of the step too*/
 38:     TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);
 39:   } else {
 40:     /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */
 41:     TSGetSolution(ts,&X);
 42:     if (!glee->Y) {VecDuplicate(X,&glee->Y);}
 43:     Y     = glee->Y;
 44:     TSEvaluateStep(ts,order-1,Y,NULL);
 45:     TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
 46:   }

 48:   if (enorm < 0) {
 49:     *accept  = PETSC_TRUE;
 50:     *next_h  = h;            /* Reuse the old step */
 51:     *wlte    = -1;           /* Weighted error was not evaluated */
 52:     *wltea   = -1;           /* Weighted absolute error was not evaluated */
 53:     *wlter   = -1;           /* Weighted relative error was not evaluated */
 54:     return(0);
 55:   }

 57:   if (enorm > 1. || enorma > 1. || enormr > 1.) {
 58:     if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
 59:     if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
 60:       PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], accepting because step size %g is at minimum\n",(double)enorm,(double)enorma,(double)enormr,(double)h);
 61:       *accept = PETSC_TRUE;
 62:     } else if (adapt->always_accept) {
 63:       PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], accepting step of size %g because always_accept is set\n",(double)enorm,(double)enorma,(double)enormr,(double)h);
 64:       *accept = PETSC_TRUE;
 65:     } else {
 66:       PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], rejecting step of size %g\n",(double)enorm,(double)enorma,(double)enormr,(double)h);
 67:       *accept = PETSC_FALSE;
 68:     }
 69:   } else {
 70:     PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative] [%g, %g, %g], accepting step of size %g\n",(double)enorm,(double)enorma,(double)enormr,(double)h);
 71:     *accept = PETSC_TRUE;
 72:   }

 74:   if (bGTEMethod){
 75:     if (*accept == PETSC_TRUE && adapt->glee_use_local) {
 76:       /* If step is accepted, then overwrite previous step error with the current error to be used on the next step */
 77:       /* WARNING: if the adapters are composable, then the accept test will not be reliable*/
 78:       TSGetTimeError(ts,0,&(glee->Y));
 79:     }

 81:     /* The optimal new step based on the current global truncation error. */
 82:     if (enorm > 0) {
 83:       /* factor based on the absolute tolerance */
 84:       hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/(order+1));
 85:       /* factor based on the relative tolerance */
 86:       hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/(order+1));
 87:       /* pick the minimum time step among the relative and absolute tolerances */
 88:       hfac_lte  = PetscMin(hfac_ltea,hfac_lter);
 89:     } else {
 90:       hfac_lte = safety * PETSC_INFINITY;
 91:     }
 92:     h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
 93:     *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
 94:   } else {
 95:     /* The optimal new step based purely on local truncation error for this step. */
 96:     if (enorm > 0) {
 97:       /* factor based on the absolute tolerance */
 98:       hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order);
 99:       /* factor based on the relative tolerance */
100:       hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order);
101:       /* pick the minimum time step among the relative and absolute tolerances */
102:       hfac_lte  = PetscMin(hfac_ltea,hfac_lter);
103:     } else {
104:       hfac_lte = safety * PETSC_INFINITY;
105:     }
106:     h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
107:     *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
108:   }
109:   *wlte   = enorm;
110:   *wltea  = enorma;
111:   *wlter  = enormr;
112:   return(0);
113: }

115: static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
116: {
117:   TSAdapt_GLEE  *glee = (TSAdapt_GLEE*)adapt->data;

121:   VecDestroy(&glee->Y);
122:   VecDestroy(&glee->E);
123:   return(0);
124: }

126: static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt)
127: {

131:   TSAdaptReset_GLEE(adapt);
132:   PetscFree(adapt->data);
133:   return(0);
134: }

136: /*MC
137:    TSADAPTGLEE - GLEE adaptive controller for time stepping

139:    Level: intermediate

141: .seealso: TS, TSAdapt, TSGetAdapt()
142: M*/
143: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
144: {
146:   TSAdapt_GLEE  *glee;

149:   PetscNewLog(adapt,&glee);
150:   adapt->data         = (void*)glee;
151:   adapt->ops->choose  = TSAdaptChoose_GLEE;
152:   adapt->ops->reset   = TSAdaptReset_GLEE;
153:   adapt->ops->destroy = TSAdaptDestroy_GLEE;
154:   return(0);
155: }