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 } TSAdapt_GLEE; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "TSAdaptChoose_GLEE" 14 static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte) 15 { 16 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 17 TSType time_scheme; /* Type of time-integration scheme */ 18 PetscErrorCode ierr; 19 Vec X,Y,Z; 20 PetscReal enorm,enorma,enormr,hfac_lte,h_lte,safety; 21 PetscInt order,order_add,stepno; 22 23 PetscFunctionBegin; 24 ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr); 25 26 safety = glee->safety; 27 ierr = TSGetType(ts,&time_scheme);CHKERRQ(ierr); 28 if (!strcmp(time_scheme,TSGLEE)){ 29 /* the method is of GLEE type */ 30 order_add=1; /* typically same order estimates */ 31 if (!glee->X) { 32 ierr = TSGetSolution(ts,&Z);CHKERRQ(ierr); 33 ierr = VecDuplicate(Z,&glee->X);CHKERRQ(ierr); 34 } 35 X = glee->X; 36 if (!glee->Y) {ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);} 37 Y = glee->Y; 38 ierr = TSGetTimeError(ts,-1,&X);CHKERRQ(ierr); 39 ierr = TSGetTimeError(ts, 0,&Y);CHKERRQ(ierr); 40 ierr = TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 41 } else { 42 /* the method is NOT of GLEE type */ 43 order_add=0; /* typically lower order estimates */ 44 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 45 if (!glee->Y) {ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);} 46 Y = glee->Y; 47 order = adapt->candidates.order[0]; 48 ierr = TSEvaluateStep(ts,order-1,Y,NULL);CHKERRQ(ierr); 49 ierr = TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 50 } 51 if (enorm > 1.) { 52 if (!*accept) safety *= glee->reject_safety; /* The last attempt also failed, shorten more aggressively */ 53 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 54 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);CHKERRQ(ierr); 55 *accept = PETSC_TRUE; 56 } else if (glee->always_accept) { 57 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);CHKERRQ(ierr); 58 *accept = PETSC_TRUE; 59 } else { 60 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr); 61 *accept = PETSC_FALSE; 62 } 63 } else { 64 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr); 65 *accept = PETSC_TRUE; 66 } 67 68 /* The optimal new step based purely on local truncation error for this step. */ 69 if (enorm == 0.0) { 70 hfac_lte = safety * PETSC_INFINITY; 71 } else { 72 hfac_lte = safety * PetscPowReal(enorm,-1./(order+order_add)); 73 } 74 h_lte = h * PetscClipInterval(hfac_lte,glee->clip[0],glee->clip[1]); 75 76 *next_sc = 0; 77 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 78 *wlte = enorm; 79 PetscFunctionReturn(0); 80 } 81 82 #undef __FUNCT__ 83 #define __FUNCT__ "TSAdaptReset_GLEE" 84 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt) 85 { 86 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 87 PetscErrorCode ierr; 88 89 PetscFunctionBegin; 90 ierr = VecDestroy(&glee->Y);CHKERRQ(ierr); 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "TSAdaptDestroy_GLEE" 96 static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt) 97 { 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 ierr = TSAdaptReset_GLEE(adapt);CHKERRQ(ierr); 102 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "TSAdaptSetFromOptions_GLEE" 108 static PetscErrorCode TSAdaptSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TSAdapt adapt) 109 { 110 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 111 PetscErrorCode ierr; 112 PetscInt two; 113 PetscBool set; 114 115 PetscFunctionBegin; 116 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE adaptive controller options");CHKERRQ(ierr); 117 two = 2; 118 ierr = PetscOptionsRealArray("-ts_adapt_glee_clip","Admissible decrease/increase in step size","",glee->clip,&two,&set);CHKERRQ(ierr); 119 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"); 120 ierr = PetscOptionsReal("-ts_adapt_glee_safety","Safety factor relative to target error","",glee->safety,&glee->safety,NULL);CHKERRQ(ierr); 121 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); 122 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); 123 ierr = PetscOptionsTail();CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "TSAdaptView_GLEE" 129 static PetscErrorCode TSAdaptView_GLEE(TSAdapt adapt,PetscViewer viewer) 130 { 131 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 132 PetscErrorCode ierr; 133 PetscBool iascii; 134 135 PetscFunctionBegin; 136 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 137 if (iascii) { 138 if (glee->always_accept) {ierr = PetscViewerASCIIPrintf(viewer," GLEE: always accepting steps\n");CHKERRQ(ierr);} 139 ierr = PetscViewerASCIIPrintf(viewer," GLEE: clip fastest decrease %g, fastest increase %g\n",(double)glee->clip[0],(double)glee->clip[1]);CHKERRQ(ierr); 140 ierr = PetscViewerASCIIPrintf(viewer," GLEE: safety factor %g, extra factor after step rejection %g\n",(double)glee->safety,(double)glee->reject_safety);CHKERRQ(ierr); 141 } 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "TSAdaptCreate_GLEE" 147 /*MC 148 TSADAPTGLEE - GLEE adaptive controller for time stepping 149 150 Level: intermediate 151 152 .seealso: TS, TSAdapt, TSSetAdapt() 153 M*/ 154 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt) 155 { 156 PetscErrorCode ierr; 157 TSAdapt_GLEE *a; 158 159 PetscFunctionBegin; 160 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr); 161 adapt->data = (void*)a; 162 adapt->ops->choose = TSAdaptChoose_GLEE; 163 adapt->ops->setfromoptions = TSAdaptSetFromOptions_GLEE; 164 adapt->ops->destroy = TSAdaptDestroy_GLEE; 165 adapt->ops->view = TSAdaptView_GLEE; 166 167 a->clip[0] = 0.1; 168 a->clip[1] = 10.; 169 a->safety = 0.9; 170 a->reject_safety = 0.5; 171 a->always_accept = PETSC_FALSE; 172 PetscFunctionReturn(0); 173 } 174