1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 3 typedef struct { 4 Vec Y; 5 } TSAdapt_Basic; 6 7 static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter) 8 { 9 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 10 PetscInt order = PETSC_DECIDE; 11 PetscReal enorm = -1; 12 PetscReal enorma,enormr; 13 PetscReal safety = adapt->safety; 14 PetscReal hfac_lte,h_lte; 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 *next_sc = 0; /* Reuse the same order scheme */ 19 *wltea = -1; /* Weighted absolute local truncation error is not used */ 20 *wlter = -1; /* Weighted relative local truncation error is not used */ 21 22 if (ts->ops->evaluatewlte) { 23 ierr = TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);CHKERRQ(ierr); 24 if (enorm >= 0 && order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order); 25 } else if (ts->ops->evaluatestep) { 26 if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered"); 27 if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n); 28 if (!basic->Y) {ierr = VecDuplicate(ts->vec_sol,&basic->Y);CHKERRQ(ierr);} 29 order = adapt->candidates.order[0]; 30 ierr = TSEvaluateStep(ts,order-1,basic->Y,NULL);CHKERRQ(ierr); 31 ierr = TSErrorWeightedNorm(ts,ts->vec_sol,basic->Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 32 } 33 34 if (enorm < 0) { 35 *accept = PETSC_TRUE; 36 *next_h = h; /* Reuse the old step */ 37 *wlte = -1; /* Weighted local truncation error was not evaluated */ 38 PetscFunctionReturn(0); 39 } 40 41 /* Determine whether the step is accepted of rejected */ 42 if (enorm > 1) { 43 if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */ 44 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 45 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);CHKERRQ(ierr); 46 *accept = PETSC_TRUE; 47 } else if (adapt->always_accept) { 48 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); 49 *accept = PETSC_TRUE; 50 } else { 51 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr); 52 *accept = PETSC_FALSE; 53 } 54 } else { 55 ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr); 56 *accept = PETSC_TRUE; 57 } 58 59 /* The optimal new step based purely on local truncation error for this step. */ 60 if (enorm > 0) 61 hfac_lte = safety * PetscPowReal(enorm,((PetscReal)-1)/order); 62 else 63 hfac_lte = safety * PETSC_INFINITY; 64 h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]); 65 66 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 67 *wlte = enorm; 68 PetscFunctionReturn(0); 69 } 70 71 static PetscErrorCode TSAdaptReset_Basic(TSAdapt adapt) 72 { 73 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 74 PetscErrorCode ierr; 75 76 PetscFunctionBegin; 77 ierr = VecDestroy(&basic->Y);CHKERRQ(ierr); 78 PetscFunctionReturn(0); 79 } 80 81 static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt) 82 { 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 ierr = TSAdaptReset_Basic(adapt);CHKERRQ(ierr); 87 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 91 /*MC 92 TSADAPTBASIC - Basic adaptive controller for time stepping 93 94 Level: intermediate 95 96 .seealso: TS, TSAdapt, TSSetAdapt() 97 M*/ 98 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 99 { 100 PetscErrorCode ierr; 101 TSAdapt_Basic *basic; 102 103 PetscFunctionBegin; 104 ierr= PetscNewLog(adapt,&basic);CHKERRQ(ierr); 105 adapt->data = (void*)basic; 106 adapt->ops->choose = TSAdaptChoose_Basic; 107 adapt->ops->reset = TSAdaptReset_Basic; 108 adapt->ops->destroy = TSAdaptDestroy_Basic; 109 PetscFunctionReturn(0); 110 } 111