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