1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2 #include <petscdm.h>
3
TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt * next_sc,PetscReal * next_h,PetscBool * accept,PetscReal * wlte,PetscReal * wltea,PetscReal * wlter)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*/
TSAdaptCreate_Basic(TSAdapt adapt)82 PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt)
83 {
84 PetscFunctionBegin;
85 adapt->ops->choose = TSAdaptChoose_Basic;
86 PetscFunctionReturn(0);
87 }
88