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