xref: /petsc/src/ts/adapt/impls/basic/adaptbasic.c (revision 1996d6cf3e56d9a238de3dc30a74d1e7e96eca44)
1 #include <private/tsimpl.h> /*I "petscts.h" I*/
2 
3 typedef struct {
4   PetscBool always_accept;
5   PetscReal clip[2];            /* admissible decrease/increase factors */
6   PetscReal safety;             /* safety factor relative to target error */
7   Vec       Y;
8 } TSAdapt_Basic;
9 
10 #undef __FUNCT__
11 #define __FUNCT__ "TSAdaptChoose_Basic"
12 static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
13 {
14   TSAdapt_Basic     *basic = (TSAdapt_Basic*)adapt->data;
15   PetscErrorCode    ierr;
16   Vec               X,Y;
17   PetscReal         enorm,hfac_lte,h_lte;
18   PetscInt          order,stepno;
19 
20   PetscFunctionBegin;
21   ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr);
22   ierr = TSGetSolution(ts,&X);CHKERRQ(ierr);
23   if (!basic->Y) {ierr = VecDuplicate(X,&basic->Y);CHKERRQ(ierr);}
24   Y = basic->Y;
25   order = adapt->candidates.order[0];
26   ierr = TSEvaluateStep(ts,order-1,Y,PETSC_NULL);CHKERRQ(ierr);
27 
28   ierr = TSErrorNormWRMS(ts,Y,&enorm);CHKERRQ(ierr);
29   if (enorm > 1.) {
30     ierr = PetscInfo1(adapt,"Estimated scaled local truncation error %G, step should be rejected\n",enorm);CHKERRQ(ierr);
31   } else {
32     ierr = PetscInfo1(adapt,"Estimated scaled local truncation error %G, step accepted\n",enorm);CHKERRQ(ierr);
33   }
34 
35   /* The optimal new step based purely on local truncation error for this step. */
36   hfac_lte = basic->safety * PetscRealPart(PetscPowScalar((PetscScalar)enorm,-1./order));
37   h_lte = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]);
38 
39   *next_sc = 0;
40   *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
41   *accept = PETSC_TRUE;
42   PetscFunctionReturn(0);
43 }
44 
45 #undef __FUNCT__
46 #define __FUNCT__ "TSAdaptDestroy_Basic"
47 static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt)
48 {
49   TSAdapt_Basic  *basic = (TSAdapt_Basic*)adapt->data;
50   PetscErrorCode ierr;
51 
52   PetscFunctionBegin;
53   ierr = VecDestroy(&basic->Y);CHKERRQ(ierr);
54   ierr = PetscFree(adapt->data);CHKERRQ(ierr);
55   PetscFunctionReturn(0);
56 }
57 
58 #undef __FUNCT__
59 #define __FUNCT__ "TSAdaptSetFromOptions_Basic"
60 static PetscErrorCode TSAdaptSetFromOptions_Basic(TSAdapt adapt)
61 {
62   TSAdapt_Basic  *basic = (TSAdapt_Basic*)adapt->data;
63   PetscErrorCode ierr;
64   PetscInt       two;
65   PetscBool      set;
66 
67   PetscFunctionBegin;
68   ierr = PetscOptionsHead("Basic adaptive controller options");CHKERRQ(ierr);
69   two = 2;
70   ierr = PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase in step size","",basic->clip,&two,&set);CHKERRQ(ierr);
71   if (set && (two != 2 || basic->clip[0] > basic->clip[1]))
72     SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip");
73   ierr = PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,PETSC_NULL);CHKERRQ(ierr);
74   ierr = PetscOptionsBool("-ts_adapt_basic_always_accept","Always accept the step regardless of whether local truncation error meets goal","",basic->always_accept,&basic->always_accept,PETSC_NULL);CHKERRQ(ierr);
75   if (!basic->always_accept) SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_SUP,"step rejection not implemented yet");
76   ierr = PetscOptionsTail();CHKERRQ(ierr);
77   PetscFunctionReturn(0);
78 }
79 
80 EXTERN_C_BEGIN
81 #undef __FUNCT__
82 #define __FUNCT__ "TSAdaptCreate_Basic"
83 /*MC
84    TSADAPTBASIC - Basic adaptive controller for time stepping
85 
86    Level: intermediate
87 
88 .seealso: TS, TSAdapt, TSSetAdapt()
89 M*/
90 PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
91 {
92   PetscErrorCode ierr;
93   TSAdapt_Basic *a;
94 
95   PetscFunctionBegin;
96   ierr = PetscNewLog(adapt,TSAdapt_Basic,&a);CHKERRQ(ierr);
97   adapt->data = (void*)a;
98   adapt->ops->choose         = TSAdaptChoose_Basic;
99   adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic;
100   adapt->ops->destroy        = TSAdaptDestroy_Basic;
101 
102   a->clip[0] = 0.1;
103   a->clip[1] = 10.;
104   a->safety  = 0.9;
105   a->always_accept = PETSC_TRUE; /* fix this */
106   PetscFunctionReturn(0);
107 }
108 EXTERN_C_END
109