Actual source code: adaptglee.c
petsc-3.11.4 2019-09-28
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->E) {VecDuplicate(X,&glee->E);}
28: E = glee->E;
29: TSGetTimeError(ts,0,&E);
30: /* this should be called with Y (the solution at the beginning of the step)*/
31: TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);
32: } else {
33: /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */
34: TSGetSolution(ts,&X);
35: if (!glee->Y) {VecDuplicate(X,&glee->Y);}
36: Y = glee->Y;
37: TSEvaluateStep(ts,order-1,Y,NULL);
38: TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
39: }
41: if (enorm < 0) {
42: *accept = PETSC_TRUE;
43: *next_h = h; /* Reuse the old step */
44: *wlte = -1; /* Weighted error was not evaluated */
45: *wltea = -1; /* Weighted absolute error was not evaluated */
46: *wlter = -1; /* Weighted relative error was not evaluated */
47: return(0);
48: }
50: if (enorm > 1. || enorma > 1. || enormr > 1.) {
51: if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
52: if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
53: 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);
54: *accept = PETSC_TRUE;
55: } else if (adapt->always_accept) {
56: 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);
57: *accept = PETSC_TRUE;
58: } else {
59: 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);
60: *accept = PETSC_FALSE;
61: }
62: } else {
63: 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);
64: *accept = PETSC_TRUE;
65: }
67: if (bGTEMethod){
68: /* The optimal new step based on the current global truncation error. */
69: if (enorm > 0) {
70: /* factor based on the absolute tolerance */
71: hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/order);
72: /* factor based on the relative tolerance */
73: hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/order);
74: /* pick the minimum time step among the relative and absolute tolerances */
75: hfac_lte = PetscMin(hfac_ltea,hfac_lter);
76: } else {
77: hfac_lte = safety * PETSC_INFINITY;
78: }
79: h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
80: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
81: } else {
82: /* The optimal new step based purely on local truncation error for this step. */
83: if (enorm > 0) {
84: /* factor based on the absolute tolerance */
85: hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order);
86: /* factor based on the relative tolerance */
87: hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order);
88: /* pick the minimum time step among the relative and absolute tolerances */
89: hfac_lte = PetscMin(hfac_ltea,hfac_lter);
90: } else {
91: hfac_lte = safety * PETSC_INFINITY;
92: }
93: h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
94: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
95: }
96: *wlte = enorm;
97: *wltea = enorma;
98: *wlter = enormr;
99: return(0);
100: }
102: static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
103: {
104: TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data;
108: VecDestroy(&glee->Y);
109: VecDestroy(&glee->E);
110: return(0);
111: }
113: static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt)
114: {
118: TSAdaptReset_GLEE(adapt);
119: PetscFree(adapt->data);
120: return(0);
121: }
123: /*MC
124: TSADAPTGLEE - GLEE adaptive controller for time stepping
126: Level: intermediate
128: .seealso: TS, TSAdapt, TSGetAdapt()
129: M*/
130: PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
131: {
133: TSAdapt_GLEE *glee;
136: PetscNewLog(adapt,&glee);
137: adapt->data = (void*)glee;
138: adapt->ops->choose = TSAdaptChoose_GLEE;
139: adapt->ops->reset = TSAdaptReset_GLEE;
140: adapt->ops->destroy = TSAdaptDestroy_GLEE;
141: return(0);
142: }