1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 typedef struct { 4 Vec E,Y; 5 } TSAdapt_GLEE; 6 7 static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter) 8 { 9 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 10 TSType time_scheme; /* Type of time-integration scheme */ 11 PetscErrorCode ierr; 12 Vec X,Y,E; 13 PetscReal enorm,enorma,enormr,hfac_lte,hfac_ltea,hfac_lter,h_lte,safety; 14 PetscInt order,stepno; 15 PetscBool bGTEMethod=PETSC_FALSE; 16 17 PetscFunctionBegin; 18 19 *next_sc = 0; /* Reuse the same order scheme */ 20 safety = adapt->safety; 21 ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr); 22 ierr = TSGetType(ts,&time_scheme);CHKERRQ(ierr); 23 if (!strcmp(time_scheme,TSGLEE)) bGTEMethod=PETSC_TRUE; 24 order = adapt->candidates.order[0]; 25 26 if (bGTEMethod){/* the method is of GLEE type */ 27 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 28 if (!glee->E) {ierr = VecDuplicate(X,&glee->E);CHKERRQ(ierr);} 29 E = glee->E; 30 ierr = TSGetTimeError(ts,0,&E);CHKERRQ(ierr); 31 /* this should be called with Y (the solution at the beginning of the step)*/ 32 ierr = TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 33 } else { 34 /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */ 35 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 36 if (!glee->Y) {ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);} 37 Y = glee->Y; 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 42 if (enorm < 0) { 43 *accept = PETSC_TRUE; 44 *next_h = h; /* Reuse the old step */ 45 *wlte = -1; /* Weighted error was not evaluated */ 46 *wltea = -1; /* Weighted absolute error was not evaluated */ 47 *wlter = -1; /* Weighted relative error was not evaluated */ 48 PetscFunctionReturn(0); 49 } 50 51 if (enorm > 1. || enorma > 1. || enormr > 1.) { 52 if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */ 53 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 54 ierr = PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], accepting because step size %g is at minimum\n",(double)enorm,(double)enorma,(double)enormr,(double)h);CHKERRQ(ierr); 55 *accept = PETSC_TRUE; 56 } else if (adapt->always_accept) { 57 ierr = PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], accepting step of size %g because always_accept is set\n",(double)enorm,(double)enorma,(double)enormr,(double)h);CHKERRQ(ierr); 58 *accept = PETSC_TRUE; 59 } else { 60 ierr = PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative]] [%g, %g, %g], rejecting step of size %g\n",(double)enorm,(double)enorma,(double)enormr,(double)h);CHKERRQ(ierr); 61 *accept = PETSC_FALSE; 62 } 63 } else { 64 ierr = PetscInfo4(adapt,"Estimated scaled truncation error [combined, absolute, relative] [%g, %g, %g], accepting step of size %g\n",(double)enorm,(double)enorma,(double)enormr,(double)h);CHKERRQ(ierr); 65 *accept = PETSC_TRUE; 66 } 67 68 if (bGTEMethod){ 69 /* The optimal new step based on the current global truncation error. */ 70 if (enorm > 0) { 71 /* factor based on the absolute tolerance */ 72 hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/order); 73 /* factor based on the relative tolerance */ 74 hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/order); 75 /* pick the minimum time step among the relative and absolute tolerances */ 76 hfac_lte = PetscMin(hfac_ltea,hfac_lter); 77 } else { 78 hfac_lte = safety * PETSC_INFINITY; 79 } 80 h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]); 81 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 82 } else { 83 /* The optimal new step based purely on local truncation error for this step. */ 84 if (enorm > 0) { 85 /* factor based on the absolute tolerance */ 86 hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order); 87 /* factor based on the relative tolerance */ 88 hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order); 89 /* pick the minimum time step among the relative and absolute tolerances */ 90 hfac_lte = PetscMin(hfac_ltea,hfac_lter); 91 } else { 92 hfac_lte = safety * PETSC_INFINITY; 93 } 94 h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]); 95 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 96 } 97 *wlte = enorm; 98 *wltea = enorma; 99 *wlter = enormr; 100 PetscFunctionReturn(0); 101 } 102 103 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt) 104 { 105 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 106 PetscErrorCode ierr; 107 108 PetscFunctionBegin; 109 ierr = VecDestroy(&glee->Y);CHKERRQ(ierr); 110 ierr = VecDestroy(&glee->E);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 114 static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt) 115 { 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 ierr = TSAdaptReset_GLEE(adapt);CHKERRQ(ierr); 120 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 /*MC 125 TSADAPTGLEE - GLEE adaptive controller for time stepping 126 127 Level: intermediate 128 129 .seealso: TS, TSAdapt, TSSetAdapt() 130 M*/ 131 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt) 132 { 133 PetscErrorCode ierr; 134 TSAdapt_GLEE *glee; 135 136 PetscFunctionBegin; 137 ierr = PetscNewLog(adapt,&glee);CHKERRQ(ierr); 138 adapt->data = (void*)glee; 139 adapt->ops->choose = TSAdaptChoose_GLEE; 140 adapt->ops->reset = TSAdaptReset_GLEE; 141 adapt->ops->destroy = TSAdaptDestroy_GLEE; 142 PetscFunctionReturn(0); 143 } 144