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 if (adapt->timestepjustincreased){ 65 hfac_lte = PetscMin(hfac_lte,1.0); 66 adapt->timestepjustincreased--; 67 } 68 h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]); 69 70 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 71 *wlte = enorm; 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode TSAdaptReset_Basic(TSAdapt adapt) 76 { 77 TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = VecDestroy(&basic->Y);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt) 86 { 87 PetscErrorCode ierr; 88 89 PetscFunctionBegin; 90 ierr = TSAdaptReset_Basic(adapt);CHKERRQ(ierr); 91 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 /*MC 96 TSADAPTBASIC - Basic adaptive controller for time stepping 97 98 Level: intermediate 99 100 .seealso: TS, TSAdapt, TSGetAdapt() 101 M*/ 102 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 103 { 104 PetscErrorCode ierr; 105 TSAdapt_Basic *basic; 106 107 PetscFunctionBegin; 108 ierr= PetscNewLog(adapt,&basic);CHKERRQ(ierr); 109 adapt->data = (void*)basic; 110 adapt->ops->choose = TSAdaptChoose_Basic; 111 adapt->ops->reset = TSAdaptReset_Basic; 112 adapt->ops->destroy = TSAdaptDestroy_Basic; 113 PetscFunctionReturn(0); 114 } 115