Actual source code: adaptglee.c
petsc-3.14.6 2021-03-30
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;
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;
18: *next_sc = 0; /* Reuse the same order scheme */
19: safety = adapt->safety;
20: PetscObjectTypeCompare((PetscObject)ts,TSGLEE,&bGTEMethod);
21: order = adapt->candidates.order[0];
23: if (bGTEMethod){/* the method is of GLEE type */
24: DM dm;
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: TSGetDM(ts,&dm);
32: DMGetGlobalVector(dm,&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: DMRestoreGlobalVector(dm,&E);
40: } else {
41: /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */
42: TSGetSolution(ts,&X);
43: if (!glee->Y) {VecDuplicate(X,&glee->Y);}
44: Y = glee->Y;
45: TSEvaluateStep(ts,order-1,Y,NULL);
46: TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
47: }
49: if (enorm < 0) {
50: *accept = PETSC_TRUE;
51: *next_h = h; /* Reuse the old step */
52: *wlte = -1; /* Weighted error was not evaluated */
53: *wltea = -1; /* Weighted absolute error was not evaluated */
54: *wlter = -1; /* Weighted relative error was not evaluated */
55: return(0);
56: }
58: if (enorm > 1. || enorma > 1. || enormr > 1.) {
59: if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */
60: if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
61: 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);
62: *accept = PETSC_TRUE;
63: } else if (adapt->always_accept) {
64: 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);
65: *accept = PETSC_TRUE;
66: } else {
67: 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);
68: *accept = PETSC_FALSE;
69: }
70: } else {
71: 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);
72: *accept = PETSC_TRUE;
73: }
75: if (bGTEMethod){
76: if (*accept == PETSC_TRUE && adapt->glee_use_local) {
77: /* If step is accepted, then overwrite previous step error with the current error to be used on the next step */
78: /* WARNING: if the adapters are composable, then the accept test will not be reliable*/
79: TSGetTimeError(ts,0,&glee->Y);
80: }
82: /* The optimal new step based on the current global truncation error. */
83: if (enorm > 0) {
84: /* factor based on the absolute tolerance */
85: hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/(order+1));
86: /* factor based on the relative tolerance */
87: hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/(order+1));
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: } else {
96: /* The optimal new step based purely on local truncation error for this step. */
97: if (enorm > 0) {
98: /* factor based on the absolute tolerance */
99: hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order);
100: /* factor based on the relative tolerance */
101: hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order);
102: /* pick the minimum time step among the relative and absolute tolerances */
103: hfac_lte = PetscMin(hfac_ltea,hfac_lter);
104: } else {
105: hfac_lte = safety * PETSC_INFINITY;
106: }
107: h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]);
108: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
109: }
110: *wlte = enorm;
111: *wltea = enorma;
112: *wlter = enormr;
113: return(0);
114: }
116: static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
117: {
118: TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data;
122: VecDestroy(&glee->Y);
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: }