Actual source code: adaptbasic.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/tsimpl.h>
3: typedef struct {
4: Vec Y;
5: } TSAdapt_Basic;
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;
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 */
22: if (ts->ops->evaluatewlte) {
23: TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);
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) {VecDuplicate(ts->vec_sol,&basic->Y);}
29: order = adapt->candidates.order[0];
30: TSEvaluateStep(ts,order-1,basic->Y,NULL);
31: TSErrorWeightedNorm(ts,ts->vec_sol,basic->Y,adapt->wnormtype,&enorm,&enorma,&enormr);
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: static PetscErrorCode TSAdaptReset_Basic(TSAdapt adapt)
76: {
77: TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
81: VecDestroy(&basic->Y);
82: return(0);
83: }
85: static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt)
86: {
90: TSAdaptReset_Basic(adapt);
91: PetscFree(adapt->data);
92: return(0);
93: }
95: /*MC
96: TSADAPTBASIC - Basic adaptive controller for time stepping
98: Level: intermediate
100: .seealso: TS, TSAdapt, TSGetAdapt()
101: M*/
102: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
103: {
105: TSAdapt_Basic *basic;
108: ierr= PetscNewLog(adapt,&basic);
109: adapt->data = (void*)basic;
110: adapt->ops->choose = TSAdaptChoose_Basic;
111: adapt->ops->reset = TSAdaptReset_Basic;
112: adapt->ops->destroy = TSAdaptDestroy_Basic;
113: return(0);
114: }