1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdm.h> 3 4 static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte,PetscReal *wltea,PetscReal *wlter) 5 { 6 Vec Y; 7 DM dm; 8 PetscInt order = PETSC_DECIDE; 9 PetscReal enorm = -1; 10 PetscReal enorma,enormr; 11 PetscReal safety = adapt->safety; 12 PetscReal hfac_lte,h_lte; 13 14 PetscFunctionBegin; 15 *next_sc = 0; /* Reuse the same order scheme */ 16 *wltea = -1; /* Weighted absolute local truncation error is not used */ 17 *wlter = -1; /* Weighted relative local truncation error is not used */ 18 19 if (ts->ops->evaluatewlte) { 20 PetscCall(TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm)); 21 PetscCheckFalse(enorm >= 0 && order < 1,PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order); 22 } else if (ts->ops->evaluatestep) { 23 PetscCheckFalse(adapt->candidates.n < 1,PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered"); 24 PetscCheck(adapt->candidates.inuse_set,PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n); 25 order = adapt->candidates.order[0]; 26 PetscCall(TSGetDM(ts,&dm)); 27 PetscCall(DMGetGlobalVector(dm,&Y)); 28 PetscCall(TSEvaluateStep(ts,order-1,Y,NULL)); 29 PetscCall(TSErrorWeightedNorm(ts,ts->vec_sol,Y,adapt->wnormtype,&enorm,&enorma,&enormr)); 30 PetscCall(DMRestoreGlobalVector(dm,&Y)); 31 } 32 33 if (enorm < 0) { 34 *accept = PETSC_TRUE; 35 *next_h = h; /* Reuse the old step */ 36 *wlte = -1; /* Weighted local truncation error was not evaluated */ 37 PetscFunctionReturn(0); 38 } 39 40 /* Determine whether the step is accepted of rejected */ 41 if (enorm > 1) { 42 if (!*accept) safety *= adapt->reject_safety; /* The last attempt also failed, shorten more aggressively */ 43 if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 44 PetscCall(PetscInfo(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h)); 45 *accept = PETSC_TRUE; 46 } else if (adapt->always_accept) { 47 PetscCall(PetscInfo(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h)); 48 *accept = PETSC_TRUE; 49 } else { 50 PetscCall(PetscInfo(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h)); 51 *accept = PETSC_FALSE; 52 } 53 } else { 54 PetscCall(PetscInfo(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h)); 55 *accept = PETSC_TRUE; 56 } 57 58 /* The optimal new step based purely on local truncation error for this step. */ 59 if (enorm > 0) 60 hfac_lte = safety * PetscPowReal(enorm,((PetscReal)-1)/order); 61 else 62 hfac_lte = safety * PETSC_INFINITY; 63 if (adapt->timestepjustdecreased) { 64 hfac_lte = PetscMin(hfac_lte,1.0); 65 adapt->timestepjustdecreased--; 66 } 67 h_lte = h * PetscClipInterval(hfac_lte,adapt->clip[0],adapt->clip[1]); 68 69 *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 70 *wlte = enorm; 71 PetscFunctionReturn(0); 72 } 73 74 /*MC 75 TSADAPTBASIC - Basic adaptive controller for time stepping 76 77 Level: intermediate 78 79 .seealso: TS, TSAdapt, TSGetAdapt() 80 M*/ 81 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 82 { 83 PetscFunctionBegin; 84 adapt->ops->choose = TSAdaptChoose_Basic; 85 PetscFunctionReturn(0); 86 } 87