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 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 *next_sc = 0; /* Reuse the same order scheme */ 17 *wltea = -1; /* Weighted absolute local truncation error is not used */ 18 *wlter = -1; /* Weighted relative local truncation error is not used */ 19 20 if (ts->ops->evaluatewlte) { 21 ierr = TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);CHKERRQ(ierr); 22 if (enorm >= 0 && order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order); 23 } else if (ts->ops->evaluatestep) { 24 if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered"); 25 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); 26 order = adapt->candidates.order[0]; 27 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 28 ierr = DMGetGlobalVector(dm,&Y);CHKERRQ(ierr); 29 ierr = TSEvaluateStep(ts,order-1,Y,NULL);CHKERRQ(ierr); 30 ierr = TSErrorWeightedNorm(ts,ts->vec_sol,Y,adapt->wnormtype,&enorm,&enorma,&enormr);CHKERRQ(ierr); 31 ierr = DMRestoreGlobalVector(dm,&Y);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->timestepjustdecreased) { 65 hfac_lte = PetscMin(hfac_lte,1.0); 66 adapt->timestepjustdecreased--; 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 /*MC 76 TSADAPTBASIC - Basic adaptive controller for time stepping 77 78 Level: intermediate 79 80 .seealso: TS, TSAdapt, TSGetAdapt() 81 M*/ 82 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 83 { 84 PetscFunctionBegin; 85 adapt->ops->choose = TSAdaptChoose_Basic; 86 PetscFunctionReturn(0); 87 } 88