1 #include <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 Vec Y; 8 } TSAdapt_Basic; 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "TSAdaptChoose_Basic" 12 static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept) 13 { 14 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 15 PetscErrorCode ierr; 16 Vec X,Y; 17 PetscReal enorm,hfac_lte,h_lte; 18 PetscInt order,stepno; 19 20 PetscFunctionBegin; 21 ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr); 22 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 23 if (!basic->Y) {ierr = VecDuplicate(X,&basic->Y);CHKERRQ(ierr);} 24 Y = basic->Y; 25 order = adapt->candidates.order[0]; 26 ierr = TSEvaluateStep(ts,order-1,Y,PETSC_NULL);CHKERRQ(ierr); 27 28 ierr = TSErrorNormWRMS(ts,Y,&enorm);CHKERRQ(ierr); 29 if (enorm > 1.) { 30 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 31 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, accepting because step size %G is at minimum\n",enorm,h);CHKERRQ(ierr); 32 *accept = PETSC_TRUE; 33 } else { 34 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, rejecting step of size %G\n",enorm,h);CHKERRQ(ierr); 35 *accept = PETSC_FALSE; 36 } 37 } else { 38 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, accepting step of size %G\n",enorm,h);CHKERRQ(ierr); 39 *accept = PETSC_TRUE; 40 } 41 42 /* The optimal new step based purely on local truncation error for this step. */ 43 hfac_lte = basic->safety * PetscRealPart(PetscPowScalar((PetscScalar)enorm,(PetscReal)(-1./order))); 44 h_lte = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]); 45 46 *next_sc = 0; 47 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 48 PetscFunctionReturn(0); 49 } 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "TSAdaptDestroy_Basic" 53 static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt) 54 { 55 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = VecDestroy(&basic->Y);CHKERRQ(ierr); 60 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "TSAdaptSetFromOptions_Basic" 66 static PetscErrorCode TSAdaptSetFromOptions_Basic(TSAdapt adapt) 67 { 68 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 69 PetscErrorCode ierr; 70 PetscInt two; 71 PetscBool set; 72 73 PetscFunctionBegin; 74 ierr = PetscOptionsHead("Basic adaptive controller options");CHKERRQ(ierr); 75 two = 2; 76 ierr = PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase in step size","",basic->clip,&two,&set);CHKERRQ(ierr); 77 if (set && (two != 2 || basic->clip[0] > basic->clip[1])) 78 SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip"); 79 ierr = PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,PETSC_NULL);CHKERRQ(ierr); 80 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,PETSC_NULL);CHKERRQ(ierr); 81 if (!basic->always_accept) SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_SUP,"step rejection not implemented yet"); 82 ierr = PetscOptionsTail();CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 EXTERN_C_BEGIN 87 #undef __FUNCT__ 88 #define __FUNCT__ "TSAdaptCreate_Basic" 89 /*MC 90 TSADAPTBASIC - Basic adaptive controller for time stepping 91 92 Level: intermediate 93 94 .seealso: TS, TSAdapt, TSSetAdapt() 95 M*/ 96 PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 97 { 98 PetscErrorCode ierr; 99 TSAdapt_Basic *a; 100 101 PetscFunctionBegin; 102 ierr = PetscNewLog(adapt,TSAdapt_Basic,&a);CHKERRQ(ierr); 103 adapt->data = (void*)a; 104 adapt->ops->choose = TSAdaptChoose_Basic; 105 adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic; 106 adapt->ops->destroy = TSAdaptDestroy_Basic; 107 108 a->clip[0] = 0.1; 109 a->clip[1] = 10.; 110 a->safety = 0.9; 111 a->always_accept = PETSC_TRUE; /* fix this */ 112 PetscFunctionReturn(0); 113 } 114 EXTERN_C_END 115