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 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 TSAdapt_GLEE *glee = (TSAdapt_GLEE *)adapt->data; 10 Vec X, Y, E; 11 PetscReal enorm, enorma, enormr, hfac_lte, hfac_ltea, hfac_lter, h_lte, safety; 12 PetscInt order; 13 PetscBool bGTEMethod; 14 15 PetscFunctionBegin; 16 *next_sc = 0; /* Reuse the same order scheme */ 17 safety = adapt->safety; 18 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSGLEE, &bGTEMethod)); 19 order = adapt->candidates.order[0]; 20 21 if (bGTEMethod) { /* the method is of GLEE type */ 22 DM dm; 23 24 PetscCall(TSGetSolution(ts, &X)); 25 if (!glee->Y && adapt->glee_use_local) { 26 PetscCall(VecDuplicate(X, &glee->Y)); /*create vector to store previous step global error*/ 27 PetscCall(VecZeroEntries(glee->Y)); /*set error to zero on the first step - may not work if error is not zero initially*/ 28 } 29 PetscCall(TSGetDM(ts, &dm)); 30 PetscCall(DMGetGlobalVector(dm, &E)); 31 PetscCall(TSGetTimeError(ts, 0, &E)); 32 33 if (adapt->glee_use_local) PetscCall(VecAXPY(E, -1.0, glee->Y)); /* local error = current error - previous step error */ 34 35 /* this should be called with the solution at the beginning of the step too*/ 36 PetscCall(TSErrorWeightedENorm(ts, E, X, X, adapt->wnormtype, &enorm, &enorma, &enormr)); 37 PetscCall(DMRestoreGlobalVector(dm, &E)); 38 } else { 39 /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */ 40 PetscCall(TSGetSolution(ts, &X)); 41 if (!glee->Y) PetscCall(VecDuplicate(X, &glee->Y)); 42 Y = glee->Y; 43 PetscCall(TSEvaluateStep(ts, order - 1, Y, NULL)); 44 PetscCall(TSErrorWeightedNorm(ts, X, Y, adapt->wnormtype, &enorm, &enorma, &enormr)); 45 } 46 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 PetscFunctionReturn(0); 54 } 55 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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 } 72 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 PetscCall(TSGetTimeError(ts, 0, &glee->Y)); 78 } 79 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 PetscFunctionReturn(0); 112 } 113 114 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt) { 115 TSAdapt_GLEE *glee = (TSAdapt_GLEE *)adapt->data; 116 117 PetscFunctionBegin; 118 PetscCall(VecDestroy(&glee->Y)); 119 PetscFunctionReturn(0); 120 } 121 122 static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt) { 123 PetscFunctionBegin; 124 PetscCall(TSAdaptReset_GLEE(adapt)); 125 PetscCall(PetscFree(adapt->data)); 126 PetscFunctionReturn(0); 127 } 128 129 /*MC 130 TSADAPTGLEE - GLEE adaptive controller for time stepping 131 132 Level: intermediate 133 134 .seealso: `TS`, `TSAdapt`, `TSGetAdapt()` 135 M*/ 136 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt) { 137 TSAdapt_GLEE *glee; 138 139 PetscFunctionBegin; 140 PetscCall(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 PetscFunctionReturn(0); 146 } 147