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