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