Actual source code: adaptbasic.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdm.h>
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;
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 */
20: if (ts->ops->evaluatewlte) {
21: TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);
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: TSGetDM(ts,&dm);
28: DMGetGlobalVector(dm,&Y);
29: TSEvaluateStep(ts,order-1,Y,NULL);
30: TSErrorWeightedNorm(ts,ts->vec_sol,Y,adapt->wnormtype,&enorm,&enorma,&enormr);
31: DMRestoreGlobalVector(dm,&Y);
32: }
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: return(0);
39: }
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: PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);
46: *accept = PETSC_TRUE;
47: } else if (adapt->always_accept) {
48: PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);
49: *accept = PETSC_TRUE;
50: } else {
51: PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);
52: *accept = PETSC_FALSE;
53: }
54: } else {
55: PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);
56: *accept = PETSC_TRUE;
57: }
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]);
70: *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
71: *wlte = enorm;
72: return(0);
73: }
75: /*MC
76: TSADAPTBASIC - Basic adaptive controller for time stepping
78: Level: intermediate
80: .seealso: TS, TSAdapt, TSGetAdapt()
81: M*/
82: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
83: {
85: adapt->ops->choose = TSAdaptChoose_Basic;
86: return(0);
87: }