1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2 #include <petscdm.h>
3
4 typedef struct {
5 Vec Y;
6 } TSAdapt_GLEE;
7
TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept,PetscReal * wlte,PetscReal * wltea,PetscReal * wlter)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 PetscErrorCode ierr;
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;
16
17 PetscFunctionBegin;
18 *next_sc = 0; /* Reuse the same order scheme */
19 safety = adapt->safety;
20 ierr = PetscObjectTypeCompare((PetscObject)ts,TSGLEE,&bGTEMethod);CHKERRQ(ierr);
21 order = adapt->candidates.order[0];
22
23 if (bGTEMethod){/* the method is of GLEE type */
24 DM dm;
25
26 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr);
27 if (!glee->Y && adapt->glee_use_local) {
28 ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);/*create vector to store previous step global error*/
29 ierr = VecZeroEntries(glee->Y);CHKERRQ(ierr); /*set error to zero on the first step - may not work if error is not zero initially*/
30 }
31 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
32 ierr = DMGetGlobalVector(dm,&E);CHKERRQ(ierr);
33 ierr = TSGetTimeError(ts,0,&E);CHKERRQ(ierr);
34
35 if (adapt->glee_use_local) {ierr = VecAXPY(E,-1.0,glee->Y);CHKERRQ(ierr);} /* local error = current error - previous step error */
36
37 /* this should be called with the solution at the beginning of the step too*/
38 ierr = TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr);
39 ierr = DMRestoreGlobalVector(dm,&E);CHKERRQ(ierr);
40 } else {
41 /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */
42 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr);
43 if (!glee->Y) {ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);}
44 Y = glee->Y;
45 ierr = TSEvaluateStep(ts,order-1,Y,NULL);CHKERRQ(ierr);
46 ierr = TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr);
47 }
48
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 PetscFunctionReturn(0);
56 }
57
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 ierr = 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);CHKERRQ(ierr);
62 *accept = PETSC_TRUE;
63 } else if (adapt->always_accept) {
64 ierr = 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);CHKERRQ(ierr);
65 *accept = PETSC_TRUE;
66 } else {
67 ierr = 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);CHKERRQ(ierr);
68 *accept = PETSC_FALSE;
69 }
70 } else {
71 ierr = 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);CHKERRQ(ierr);
72 *accept = PETSC_TRUE;
73 }
74
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 ierr = TSGetTimeError(ts,0,&glee->Y);CHKERRQ(ierr);
80 }
81
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 PetscFunctionReturn(0);
114 }
115
TSAdaptReset_GLEE(TSAdapt adapt)116 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt)
117 {
118 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data;
119 PetscErrorCode ierr;
120
121 PetscFunctionBegin;
122 ierr = VecDestroy(&glee->Y);CHKERRQ(ierr);
123 PetscFunctionReturn(0);
124 }
125
TSAdaptDestroy_GLEE(TSAdapt adapt)126 static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt)
127 {
128 PetscErrorCode ierr;
129
130 PetscFunctionBegin;
131 ierr = TSAdaptReset_GLEE(adapt);CHKERRQ(ierr);
132 ierr = PetscFree(adapt->data);CHKERRQ(ierr);
133 PetscFunctionReturn(0);
134 }
135
136 /*MC
137 TSADAPTGLEE - GLEE adaptive controller for time stepping
138
139 Level: intermediate
140
141 .seealso: TS, TSAdapt, TSGetAdapt()
142 M*/
TSAdaptCreate_GLEE(TSAdapt adapt)143 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt)
144 {
145 PetscErrorCode ierr;
146 TSAdapt_GLEE *glee;
147
148 PetscFunctionBegin;
149 ierr = PetscNewLog(adapt,&glee);CHKERRQ(ierr);
150 adapt->data = (void*)glee;
151 adapt->ops->choose = TSAdaptChoose_GLEE;
152 adapt->ops->reset = TSAdaptReset_GLEE;
153 adapt->ops->destroy = TSAdaptDestroy_GLEE;
154 PetscFunctionReturn(0);
155 }
156