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