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