Actual source code: adaptglee.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdm.h>
4: typedef struct {
5: Vec Y;
6: } TSAdapt_GLEE;
8: static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter)
9: {
10: TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data;
11: Vec X,Y,E;
12: PetscReal enorm,enorma,enormr,hfac_lte,hfac_ltea,hfac_lter,h_lte,safety;
13: PetscInt order;
14: PetscBool bGTEMethod;
16: *next_sc = 0; /* Reuse the same order scheme */
17: safety = adapt->safety;
18: PetscObjectTypeCompare((PetscObject)ts,TSGLEE,&bGTEMethod);
19: order = adapt->candidates.order[0];
21: if (bGTEMethod) {/* the method is of GLEE type */
22: DM dm;
24: TSGetSolution(ts,&X);
25: if (!glee->Y && adapt->glee_use_local) {
26: VecDuplicate(X,&glee->Y);/*create vector to store previous step global error*/
27: VecZeroEntries(glee->Y); /*set error to zero on the first step - may not work if error is not zero initially*/
28: }
29: TSGetDM(ts,&dm);
30: DMGetGlobalVector(dm,&E);
31: TSGetTimeError(ts,0,&E);
33: if (adapt->glee_use_local) VecAXPY(E,-1.0,glee->Y); /* local error = current error - previous step error */
35: /* this should be called with the solution at the beginning of the step too*/
36: TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);
37: DMRestoreGlobalVector(dm,&E);
38: } else {
39: /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */
40: TSGetSolution(ts,&X);
41: if (!glee->Y) VecDuplicate(X,&glee->Y);
42: Y = glee->Y;
43: TSEvaluateStep(ts,order-1,Y,NULL);
44: TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
45: }
47: if (enorm < 0) {
48: *accept = PETSC_TRUE;
49: *next_h = h; /* Reuse the old step */
50: *wlte = -1; /* Weighted error was not evaluated */
51: *wltea = -1; /* Weighted absolute error was not evaluated */
52: *wlter = -1; /* Weighted relative error was not evaluated */
53: return 0;
54: }
56: if (enorm > 1. || enorma > 1. || enormr > 1.) {
57: if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
58: if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
59: PetscInfo(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);
60: *accept = PETSC_TRUE;
61: } else if (adapt->always_accept) {
62: PetscInfo(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);
63: *accept = PETSC_TRUE;
64: } else {
65: PetscInfo(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);
66: *accept = PETSC_FALSE;
67: }
68: } else {
69: PetscInfo(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);
70: *accept = PETSC_TRUE;
71: }
73: if (bGTEMethod) {
74: if (*accept == PETSC_TRUE && adapt->glee_use_local) {
75: /* If step is accepted, then overwrite previous step error with the current error to be used on the next step */
76: /* WARNING: if the adapters are composable, then the accept test will not be reliable*/
77: TSGetTimeError(ts,0,&glee->Y);
78: }
80: /* The optimal new step based on the current global truncation error. */
81: if (enorm > 0) {
82: /* factor based on the absolute tolerance */
83: hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/(order+1));
84: /* factor based on the relative tolerance */
85: hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/(order+1));
86: /* pick the minimum time step among the relative and absolute tolerances */
87: hfac_lte = PetscMin(hfac_ltea,hfac_lter);
88: } else {
89: hfac_lte = safety * PETSC_INFINITY;
90: }
91: h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
92: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
93: } else {
94: /* The optimal new step based purely on local truncation error for this step. */
95: if (enorm > 0) {
96: /* factor based on the absolute tolerance */
97: hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order);
98: /* factor based on the relative tolerance */
99: hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order);
100: /* pick the minimum time step among the relative and absolute tolerances */
101: hfac_lte = PetscMin(hfac_ltea,hfac_lter);
102: } else {
103: hfac_lte = safety * PETSC_INFINITY;
104: }
105: h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
106: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
107: }
108: *wlte = enorm;
109: *wltea = enorma;
110: *wlter = enormr;
111: return 0;
112: }
114: static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
115: {
116: TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data;
118: VecDestroy(&glee->Y);
119: return 0;
120: }
122: static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt)
123: {
124: TSAdaptReset_GLEE(adapt);
125: PetscFree(adapt->data);
126: return 0;
127: }
129: /*MC
130: TSADAPTGLEE - GLEE adaptive controller for time stepping
132: Level: intermediate
134: .seealso: TS, TSAdapt, TSGetAdapt()
135: M*/
136: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
137: {
138: TSAdapt_GLEE *glee;
140: PetscNewLog(adapt,&glee);
141: adapt->data = (void*)glee;
142: adapt->ops->choose = TSAdaptChoose_GLEE;
143: adapt->ops->reset = TSAdaptReset_GLEE;
144: adapt->ops->destroy = TSAdaptDestroy_GLEE;
145: return 0;
146: }