1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdm.h> 3 4 typedef struct { 5 Vec Y; 6 } TSAdapt_GLEE; 7 8 static PetscErrorCode TSAdaptChoose_GLEE(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter) 9 { 10 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 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; 16 17 PetscFunctionBegin; 18 *next_sc = 0; /* Reuse the same order scheme */ 19 safety = adapt->safety; 20 ierr = PetscObjectTypeCompare((PetscObject)ts,TSGLEE,&bGTEMethod);CHKERRQ(ierr); 21 order = adapt->candidates.order[0]; 22 23 if (bGTEMethod) {/* the method is of GLEE type */ 24 DM dm; 25 26 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 27 if (!glee->Y && adapt->glee_use_local) { 28 ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);/*create vector to store previous step global error*/ 29 ierr = VecZeroEntries(glee->Y);CHKERRQ(ierr); /*set error to zero on the first step - may not work if error is not zero initially*/ 30 } 31 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 32 ierr = DMGetGlobalVector(dm,&E);CHKERRQ(ierr); 33 ierr = TSGetTimeError(ts,0,&E);CHKERRQ(ierr); 34 35 if (adapt->glee_use_local) {ierr = VecAXPY(E,-1.0,glee->Y);CHKERRQ(ierr);} /* local error = current error - previous step error */ 36 37 /* this should be called with the solution at the beginning of the step too*/ 38 ierr = TSErrorWeightedENorm(ts,E,X,X,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 39 ierr = DMRestoreGlobalVector(dm,&E);CHKERRQ(ierr); 40 } else { 41 /* the method is NOT of GLEE type; use the stantard basic augmented by separate atol and rtol */ 42 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 43 if (!glee->Y) {ierr = VecDuplicate(X,&glee->Y);CHKERRQ(ierr);} 44 Y = glee->Y; 45 ierr = TSEvaluateStep(ts,order-1,Y,NULL);CHKERRQ(ierr); 46 ierr = TSErrorWeightedNorm(ts,X,Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 47 } 48 49 if (enorm < 0) { 50 *accept = PETSC_TRUE; 51 *next_h = h; /* Reuse the old step */ 52 *wlte = -1; /* Weighted error was not evaluated */ 53 *wltea = -1; /* Weighted absolute error was not evaluated */ 54 *wlter = -1; /* Weighted relative error was not evaluated */ 55 PetscFunctionReturn(0); 56 } 57 58 if (enorm > 1. || enorma > 1. || enormr > 1.) { 59 if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */ 60 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 61 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); 62 *accept = PETSC_TRUE; 63 } else if (adapt->always_accept) { 64 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); 65 *accept = PETSC_TRUE; 66 } else { 67 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); 68 *accept = PETSC_FALSE; 69 } 70 } else { 71 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); 72 *accept = PETSC_TRUE; 73 } 74 75 if (bGTEMethod) { 76 if (*accept == PETSC_TRUE && adapt->glee_use_local) { 77 /* If step is accepted, then overwrite previous step error with the current error to be used on the next step */ 78 /* WARNING: if the adapters are composable, then the accept test will not be reliable*/ 79 ierr = TSGetTimeError(ts,0,&glee->Y);CHKERRQ(ierr); 80 } 81 82 /* The optimal new step based on the current global truncation error. */ 83 if (enorm > 0) { 84 /* factor based on the absolute tolerance */ 85 hfac_ltea = safety * PetscPowReal(1./enorma,((PetscReal)1)/(order+1)); 86 /* factor based on the relative tolerance */ 87 hfac_lter = safety * PetscPowReal(1./enormr,((PetscReal)1)/(order+1)); 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 } else { 96 /* The optimal new step based purely on local truncation error for this step. */ 97 if (enorm > 0) { 98 /* factor based on the absolute tolerance */ 99 hfac_ltea = safety * PetscPowReal(enorma,((PetscReal)-1)/order); 100 /* factor based on the relative tolerance */ 101 hfac_lter = safety * PetscPowReal(enormr,((PetscReal)-1)/order); 102 /* pick the minimum time step among the relative and absolute tolerances */ 103 hfac_lte = PetscMin(hfac_ltea,hfac_lter); 104 } else { 105 hfac_lte = safety * PETSC_INFINITY; 106 } 107 h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]); 108 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 109 } 110 *wlte = enorm; 111 *wltea = enorma; 112 *wlter = enormr; 113 PetscFunctionReturn(0); 114 } 115 116 static PetscErrorCode TSAdaptReset_GLEE(TSAdapt adapt) 117 { 118 TSAdapt_GLEE *glee = (TSAdapt_GLEE*)adapt->data; 119 PetscErrorCode ierr; 120 121 PetscFunctionBegin; 122 ierr = VecDestroy(&glee->Y);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode TSAdaptDestroy_GLEE(TSAdapt adapt) 127 { 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 ierr = TSAdaptReset_GLEE(adapt);CHKERRQ(ierr); 132 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 /*MC 137 TSADAPTGLEE - GLEE adaptive controller for time stepping 138 139 Level: intermediate 140 141 .seealso: TS, TSAdapt, TSGetAdapt() 142 M*/ 143 PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt adapt) 144 { 145 PetscErrorCode ierr; 146 TSAdapt_GLEE *glee; 147 148 PetscFunctionBegin; 149 ierr = PetscNewLog(adapt,&glee);CHKERRQ(ierr); 150 adapt->data = (void*)glee; 151 adapt->ops->choose = TSAdaptChoose_GLEE; 152 adapt->ops->reset = TSAdaptReset_GLEE; 153 adapt->ops->destroy = TSAdaptDestroy_GLEE; 154 PetscFunctionReturn(0); 155 } 156