1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 typedef struct { 4 PetscReal clip[2]; /* admissible decrease/increase factors */ 5 PetscReal safety; /* safety factor relative to target error */ 6 PetscReal reject_safety; /* extra safety factor if the last step was rejected */ 7 Vec X; 8 Vec Y; 9 Vec E; 10 } TSAdapt_GLEE; 11 12 static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter) 13 { 14 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 15 TSType time_scheme; /* Type of time-integration scheme */ 16 PetscErrorCode ierr; 17 Vec X,Y,E; 18 PetscReal enorm,enorma,enormr,hfac_lte,hfac_ltea,hfac_lter,h_lte,safety; 19 PetscInt order,stepno; 20 PetscBool bGTEMethod=PETSC_FALSE; 21 22 PetscFunctionBegin; 23 24 *next_sc = 0; /* Reuse the same order scheme */ 25 safety = glee->safety; 26 ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr); 27 ierr = TSGetType(ts,&time_scheme);CHKERRQ(ierr); 28 if (!strcmp(time_scheme,TSGLEE)) bGTEMethod=PETSC_TRUE; 29 order = adapt->candidates.order[0]; 30 31 if (bGTEMethod){/* the method is of GLEE type */ 32 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 33 if (!glee->E) {ierr = VecDuplicate(X,&glee->E);CHKERRQ(ierr);} 34 E = glee->E; 35 ierr = TSGetTimeError(ts,0,&E);CHKERRQ(ierr); 36 /* this should be called with Y (the solution at the beginning of the step)*/ 37 ierr = TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 38 } else { 39 /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */ 40 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 41 if (!glee->Y) {ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);} 42 Y = glee->Y; 43 ierr = TSEvaluateStep(ts,order-1,Y,NULL);CHKERRQ(ierr); 44 ierr = TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 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 *= glee->reject_safety; /* The last attempt also failed, shorten more aggressively */ 58 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 59 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); 60 *accept = PETSC_TRUE; 61 } else if (adapt->always_accept) { 62 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); 63 *accept = PETSC_TRUE; 64 } else { 65 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); 66 *accept = PETSC_FALSE; 67 } 68 } else { 69 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); 70 *accept = PETSC_TRUE; 71 } 72 73 if (bGTEMethod){ 74 /* The optimal new step based on the current global truncation error. */ 75 if (enorm > 0) { 76 /* factor based on the absolute tolerance */ 77 hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/order); 78 /* factor based on the relative tolerance */ 79 hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/order); 80 /* pick the minimum time step among the relative and absolute tolerances */ 81 hfac_lte = PetscMin(hfac_ltea,hfac_lter); 82 } else { 83 hfac_lte = safety * PETSC_INFINITY; 84 } 85 h_lte = h * PetscClipInterval(hfac_lte,glee->clip[0],glee->clip[1]); 86 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 87 } else { 88 /* The optimal new step based purely on local truncation error for this step. */ 89 if (enorm > 0) { 90 /* factor based on the absolute tolerance */ 91 hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order); 92 /* factor based on the relative tolerance */ 93 hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order); 94 /* pick the minimum time step among the relative and absolute tolerances */ 95 hfac_lte = PetscMin(hfac_ltea,hfac_lter); 96 } else { 97 hfac_lte = safety * PETSC_INFINITY; 98 } 99 h_lte = h * PetscClipInterval(hfac_lte,glee->clip[0],glee->clip[1]); 100 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 101 } 102 *wlte = enorm; 103 *wltea = enorma; 104 *wlter = enormr; 105 PetscFunctionReturn(0); 106 } 107 108 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt) 109 { 110 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = VecDestroy(&glee->Y);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt) 119 { 120 PetscErrorCode ierr; 121 122 PetscFunctionBegin; 123 ierr = TSAdaptReset_GLEE(adapt);CHKERRQ(ierr); 124 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 static PetscErrorCode TSAdaptSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TSAdapt adapt) 129 { 130 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 131 PetscErrorCode ierr; 132 PetscInt two; 133 PetscBool set; 134 135 PetscFunctionBegin; 136 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE adaptive controller options");CHKERRQ(ierr); 137 two = 2; 138 ierr = PetscOptionsRealArray("-ts_adapt_glee_clip","Admissible decrease/increase in step size","",glee->clip,&two,&set);CHKERRQ(ierr); 139 if (set && (two != 2 || glee->clip[0] > glee->clip[1])) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_glee_clip"); 140 ierr = PetscOptionsReal("-ts_adapt_glee_safety","Safety factor relative to target error","",glee->safety,&glee->safety,NULL);CHKERRQ(ierr); 141 ierr = PetscOptionsReal("-ts_adapt_glee_reject_safety","Extra safety factor to apply if the last step was rejected","",glee->reject_safety,&glee->reject_safety,NULL);CHKERRQ(ierr); 142 ierr = PetscOptionsTail();CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 146 static PetscErrorCode TSAdaptView_GLEE(TSAdapt adapt,PetscViewer viewer) 147 { 148 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 149 PetscErrorCode ierr; 150 PetscBool iascii; 151 152 PetscFunctionBegin; 153 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 154 if (iascii) { 155 ierr = PetscViewerASCIIPrintf(viewer," GLEE: clip fastest decrease %g, fastest increase %g\n",(double)glee->clip[0],(double)glee->clip[1]);CHKERRQ(ierr); 156 ierr = PetscViewerASCIIPrintf(viewer," GLEE: safety factor %g, extra factor after step rejection %g\n",(double)glee->safety,(double)glee->reject_safety);CHKERRQ(ierr); 157 } 158 PetscFunctionReturn(0); 159 } 160 161 /*MC 162 TSADAPTGLEE - GLEE adaptive controller for time stepping 163 164 Level: intermediate 165 166 .seealso: TS, TSAdapt, TSSetAdapt() 167 M*/ 168 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt) 169 { 170 PetscErrorCode ierr; 171 TSAdapt_GLEE *a; 172 173 PetscFunctionBegin; 174 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr); 175 adapt->data = (void*)a; 176 adapt->ops->choose = TSAdaptChoose_GLEE; 177 adapt->ops->setfromoptions = TSAdaptSetFromOptions_GLEE; 178 adapt->ops->destroy = TSAdaptDestroy_GLEE; 179 adapt->ops->view = TSAdaptView_GLEE; 180 181 a->clip[0] = 0.1; 182 a->clip[1] = 10.; 183 a->safety = 0.9; 184 a->reject_safety = 0.5; 185 PetscFunctionReturn(0); 186 } 187