Actual source code: adaptglee.c
petsc-3.12.5 2020-03-29
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: }