1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdm.h> 3 4 static const char *citation[] = { 5 "@article{Soderlind2003,\n" 6 " author = {S\"{o}derlind, Gustaf},\n" 7 " title = {Digital Filters in Adaptive Time-stepping},\n" 8 " journal = {ACM Transactions on Mathematical Software},\n" 9 " volume = {29},\n" 10 " number = {1},\n" 11 " pages = {1--26},\n" 12 " year = {2003},\n" 13 " issn = {0098-3500},\n" 14 " doi = {http://dx.doi.org/10.1145/641876.641877},\n" 15 "}\n", 16 "@article{Soderlind2006,\n" 17 " author = {Gustaf S\"{o}derlind and Lina Wang},\n" 18 " title = {Adaptive time-stepping and computational stability},\n" 19 " journal = {Journal of Computational and Applied Mathematics},\n" 20 " volume = {185},\n" 21 " number = {2},\n" 22 " pages = {225--243},\n" 23 " year = {2006},\n" 24 " issn = {0377-0427},\n" 25 " doi = {http://dx.doi.org/10.1016/j.cam.2005.03.008},\n" 26 "}\n", 27 }; 28 static PetscBool cited[] = {PETSC_FALSE, PETSC_FALSE}; 29 30 typedef struct { 31 PetscReal kBeta[3]; /* filter parameters */ 32 PetscReal Alpha[2]; /* filter parameters */ 33 PetscReal cerror[3]; /* control error (controller input) history */ 34 PetscReal hratio[3]; /* stepsize ratio (controller output) history */ 35 PetscBool rollback; 36 } TSAdapt_DSP; 37 38 static PetscReal Limiter(PetscReal value, PetscReal kappa) { 39 return 1 + kappa * PetscAtanReal((value - 1) / kappa); 40 } 41 42 static PetscErrorCode TSAdaptRestart_DSP(TSAdapt adapt) { 43 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 44 PetscFunctionBegin; 45 dsp->cerror[0] = dsp->hratio[0] = 1.0; 46 dsp->cerror[1] = dsp->hratio[1] = 1.0; 47 dsp->cerror[2] = dsp->hratio[2] = 1.0; 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode TSAdaptRollBack_DSP(TSAdapt adapt) { 52 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 53 PetscFunctionBegin; 54 dsp->cerror[0] = dsp->cerror[1]; 55 dsp->cerror[1] = dsp->cerror[2]; 56 dsp->cerror[2] = 1.0; 57 dsp->hratio[0] = dsp->hratio[1]; 58 dsp->hratio[1] = dsp->hratio[2]; 59 dsp->hratio[2] = 1.0; 60 PetscFunctionReturn(0); 61 } 62 63 static PetscErrorCode TSAdaptChoose_DSP(TSAdapt adapt, TS ts, PetscReal h, PetscInt *next_sc, PetscReal *next_h, PetscBool *accept, PetscReal *wlte, PetscReal *wltea, PetscReal *wlter) { 64 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 65 PetscInt order = PETSC_DECIDE; 66 PetscReal enorm = -1; 67 PetscReal enorma, enormr; 68 PetscReal safety = adapt->safety * (PetscReal)0.9; 69 PetscReal hnew, hfac = PETSC_INFINITY; 70 PetscReal hmin = adapt->dt_min * (1 + PETSC_SQRT_MACHINE_EPSILON); 71 72 PetscFunctionBegin; 73 *next_sc = 0; /* Reuse the same order scheme */ 74 *wltea = -1; /* Weighted absolute local truncation error is not used */ 75 *wlter = -1; /* Weighted relative local truncation error is not used */ 76 77 if (ts->ops->evaluatewlte) { 78 PetscCall(TSEvaluateWLTE(ts, adapt->wnormtype, &order, &enorm)); 79 PetscCheck(enorm < 0 || order >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_OUTOFRANGE, "Computed error order %" PetscInt_FMT " must be positive", order); 80 } else if (ts->ops->evaluatestep) { 81 DM dm; 82 Vec Y; 83 84 PetscCheck(adapt->candidates.n >= 1, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONGSTATE, "No candidate has been registered"); 85 PetscCheck(adapt->candidates.inuse_set, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONGSTATE, "The current in-use scheme is not among the %" PetscInt_FMT " candidates", adapt->candidates.n); 86 order = adapt->candidates.order[0]; 87 PetscCall(TSGetDM(ts, &dm)); 88 PetscCall(DMGetGlobalVector(dm, &Y)); 89 PetscCall(TSEvaluateStep(ts, order - 1, Y, NULL)); 90 PetscCall(TSErrorWeightedNorm(ts, ts->vec_sol, Y, adapt->wnormtype, &enorm, &enorma, &enormr)); 91 PetscCall(DMRestoreGlobalVector(dm, &Y)); 92 } 93 if (enorm < 0) { 94 PetscCall(TSAdaptRestart_DSP(adapt)); 95 *accept = PETSC_TRUE; /* Accept the step */ 96 *next_h = h; /* Reuse the old step size */ 97 *wlte = -1; /* Weighted local truncation error was not evaluated */ 98 PetscFunctionReturn(0); 99 } 100 101 PetscCall(PetscCitationsRegister(citation[0], &cited[0])); 102 PetscCall(PetscCitationsRegister(citation[1], &cited[1])); 103 104 /* Update history after rollback */ 105 if (!ts->steprollback) dsp->rollback = PETSC_FALSE; 106 else if (!dsp->rollback) { 107 dsp->rollback = PETSC_TRUE; 108 PetscCall(TSAdaptRollBack_DSP(adapt)); 109 } 110 /* Reset history after restart */ 111 if (ts->steprestart) PetscCall(TSAdaptRestart_DSP(adapt)); 112 113 { 114 PetscReal k = (PetscReal)order; 115 PetscReal b1 = dsp->kBeta[0]; 116 PetscReal b2 = dsp->kBeta[1]; 117 PetscReal b3 = dsp->kBeta[2]; 118 PetscReal a2 = dsp->Alpha[0]; 119 PetscReal a3 = dsp->Alpha[1]; 120 121 PetscReal ctr0; 122 PetscReal ctr1 = dsp->cerror[0]; 123 PetscReal ctr2 = dsp->cerror[1]; 124 PetscReal rho0; 125 PetscReal rho1 = dsp->hratio[0]; 126 PetscReal rho2 = dsp->hratio[1]; 127 128 /* Compute the step size ratio */ 129 enorm = PetscMax(enorm, PETSC_SMALL); 130 ctr0 = PetscPowReal(1 / enorm, 1 / k); 131 rho0 = PetscPowReal(ctr0, b1); 132 rho0 *= PetscPowReal(ctr1, b2); 133 rho0 *= PetscPowReal(ctr2, b3); 134 rho0 *= PetscPowReal(rho1, -a2); 135 rho0 *= PetscPowReal(rho2, -a3); 136 rho0 = Limiter(rho0, 1); 137 138 /* Determine whether the step is accepted or rejected */ 139 if (rho0 >= safety) *accept = PETSC_TRUE; 140 else if (adapt->always_accept) *accept = PETSC_TRUE; 141 else if (h < hmin) *accept = PETSC_TRUE; 142 else *accept = PETSC_FALSE; 143 144 /* Update history after accept */ 145 if (*accept) { 146 dsp->cerror[2] = dsp->cerror[1]; 147 dsp->cerror[1] = dsp->cerror[0]; 148 dsp->cerror[0] = ctr0; 149 dsp->hratio[2] = dsp->hratio[1]; 150 dsp->hratio[1] = dsp->hratio[0]; 151 dsp->hratio[0] = rho0; 152 dsp->rollback = PETSC_FALSE; 153 } 154 155 hfac = rho0; 156 } 157 158 hnew = h * PetscClipInterval(hfac, adapt->clip[0], adapt->clip[1]); 159 *next_h = PetscClipInterval(hnew, adapt->dt_min, adapt->dt_max); 160 *wlte = enorm; 161 PetscFunctionReturn(0); 162 } 163 164 static PetscErrorCode TSAdaptDestroy_DSP(TSAdapt adapt) { 165 PetscFunctionBegin; 166 PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", NULL)); 167 PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", NULL)); 168 PetscCall(PetscFree(adapt->data)); 169 PetscFunctionReturn(0); 170 } 171 172 static PetscErrorCode TSAdaptView_DSP(TSAdapt adapt, PetscViewer viewer) { 173 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 174 PetscBool iascii; 175 176 PetscFunctionBegin; 177 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 178 if (iascii) { 179 double a2 = (double)dsp->Alpha[0], a3 = (double)dsp->Alpha[1]; 180 double b1 = (double)dsp->kBeta[0], b2 = (double)dsp->kBeta[1], b3 = (double)dsp->kBeta[2]; 181 PetscCall(PetscViewerASCIIPrintf(viewer, "filter parameters kBeta=[%g,%g,%g] Alpha=[%g,%g]\n", b1, b2, b3, a2, a3)); 182 } 183 PetscFunctionReturn(0); 184 } 185 186 struct FilterTab { 187 const char *name; 188 PetscReal scale; 189 PetscReal kBeta[3]; 190 PetscReal Alpha[2]; 191 }; 192 193 static struct FilterTab filterlist[] = { 194 {"basic", 1, {1, 0, 0}, {0, 0} }, 195 196 {"PI30", 3, {1, 0, 0}, {0, 0} }, 197 {"PI42", 5, {3, -1, 0}, {0, 0} }, 198 {"PI33", 3, {2, -1, 0}, {0, 0} }, 199 {"PI34", 10, {7, -4, 0}, {0, 0} }, 200 201 {"PC11", 1, {2, -1, 0}, {-1, 0} }, 202 {"PC47", 10, {11, -7, 0}, {-10, 0} }, 203 {"PC36", 10, {9, -6, 0}, {-10, 0} }, 204 205 {"H0211", 2, {1, 1, 0}, {1, 0} }, 206 {"H211b", 4, {1, 1, 0}, {1, 0} }, 207 {"H211PI", 6, {1, 1, 0}, {0, 0} }, 208 209 {"H0312", 4, {1, 2, 1}, {3, 1} }, 210 {"H312b", 8, {1, 2, 1}, {3, 1} }, 211 {"H312PID", 18, {1, 2, 1}, {0, 0} }, 212 213 {"H0321", 4, {5, 2, -3}, {-1, -3} }, 214 {"H321", 18, {6, 1, -5}, {-15, -3}}, 215 }; 216 217 static PetscErrorCode TSAdaptDSPSetFilter_DSP(TSAdapt adapt, const char *name) { 218 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 219 PetscInt i, count = (PetscInt)(sizeof(filterlist) / sizeof(filterlist[0])); 220 struct FilterTab *tab = NULL; 221 PetscBool match; 222 223 PetscFunctionBegin; 224 for (i = 0; i < count; i++) { 225 PetscCall(PetscStrcasecmp(name, filterlist[i].name, &match)); 226 if (match) { 227 tab = &filterlist[i]; 228 break; 229 } 230 } 231 PetscCheck(tab, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_UNKNOWN_TYPE, "Filter name %s not found", name); 232 dsp->kBeta[0] = tab->kBeta[0] / tab->scale; 233 dsp->kBeta[1] = tab->kBeta[1] / tab->scale; 234 dsp->kBeta[2] = tab->kBeta[2] / tab->scale; 235 dsp->Alpha[0] = tab->Alpha[0] / tab->scale; 236 dsp->Alpha[1] = tab->Alpha[1] / tab->scale; 237 PetscFunctionReturn(0); 238 } 239 240 static PetscErrorCode TSAdaptDSPSetPID_DSP(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD) { 241 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 242 243 PetscFunctionBegin; 244 dsp->kBeta[0] = kkI + kkP + kkD; 245 dsp->kBeta[1] = -(kkP + 2 * kkD); 246 dsp->kBeta[2] = kkD; 247 dsp->Alpha[0] = 0; 248 dsp->Alpha[1] = 0; 249 PetscFunctionReturn(0); 250 } 251 252 static PetscErrorCode TSAdaptSetFromOptions_DSP(TSAdapt adapt, PetscOptionItems *PetscOptionsObject) { 253 TSAdapt_DSP *dsp = (TSAdapt_DSP *)adapt->data; 254 const char *names[sizeof(filterlist) / sizeof(filterlist[0])]; 255 PetscInt count = (PetscInt)(sizeof(filterlist) / sizeof(filterlist[0])); 256 PetscInt index = 2; /* PI42 */ 257 PetscReal pid[3] = {1, 0, 0}; 258 PetscInt i, n; 259 PetscBool set; 260 261 PetscFunctionBegin; 262 for (i = 0; i < count; i++) names[i] = filterlist[i].name; 263 PetscOptionsHeadBegin(PetscOptionsObject, "DSP adaptive controller options"); 264 265 PetscCall(PetscOptionsEList("-ts_adapt_dsp_filter", "Filter name", "TSAdaptDSPSetFilter", names, count, names[index], &index, &set)); 266 if (set) PetscCall(TSAdaptDSPSetFilter(adapt, names[index])); 267 268 PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_pid", "PID parameters <kkI,kkP,kkD>", "TSAdaptDSPSetPID", pid, (n = 3, &n), &set)); 269 PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for PID parameters"); 270 if (set) PetscCall(TSAdaptDSPSetPID(adapt, pid[0], pid[1], pid[2])); 271 272 PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_kbeta", "Filter parameters", "", dsp->kBeta, (n = 3, &n), &set)); 273 PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for parameter kbeta"); 274 if (set) 275 for (i = n; i < 3; i++) dsp->kBeta[i] = 0; 276 277 PetscCall(PetscOptionsRealArray("-ts_adapt_dsp_alpha", "Filter parameters", "", dsp->Alpha, (n = 2, &n), &set)); 278 PetscCheck(!set || n, PetscObjectComm((PetscObject)adapt), PETSC_ERR_ARG_WRONG, "Must provide at least one value for parameter alpha"); 279 if (set) 280 for (i = n; i < 2; i++) dsp->Alpha[i] = 0; 281 282 PetscOptionsHeadEnd(); 283 PetscFunctionReturn(0); 284 } 285 286 /*@C 287 TSAdaptDSPSetFilter - Sets internal parameters corresponding to the named filter 288 289 Collective on TSAdapt 290 291 Input Parameters: 292 + adapt - adaptive controller context 293 - name - filter name 294 295 Level: intermediate 296 297 References: 298 . * - http://dx.doi.org/10.1145/641876.641877 299 300 Notes: 301 Valid filter names are 302 + "basic" - similar to TSADAPTBASIC but with different criteria for step rejections. 303 . "PI30", "PI42", "PI33", "PI34" - PI controllers. 304 . "PC11", "PC47", "PC36" - predictive controllers. 305 . "H0211", "H211b", "H211PI" - digital filters with orders dynamics=2, adaptivity=1, filter=1. 306 . "H0312", "H312b", "H312PID" - digital filters with orders dynamics=3, adaptivity=1, filter=2. 307 - "H0321", "H321" - digital filters with orders dynamics=3, adaptivity=2, filter=1. 308 309 Options Database: 310 . -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers 311 312 .seealso: `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()` 313 @*/ 314 PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt, const char *name) { 315 PetscFunctionBegin; 316 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 317 PetscValidCharPointer(name, 2); 318 PetscTryMethod(adapt, "TSAdaptDSPSetFilter_C", (TSAdapt, const char *), (adapt, name)); 319 PetscFunctionReturn(0); 320 } 321 322 /*@ 323 TSAdaptDSPSetPID - Set the PID controller parameters 324 325 Input Parameters: 326 + adapt - adaptive controller context 327 . kkI - Integral parameter 328 . kkP - Proportional parameter 329 - kkD - Derivative parameter 330 331 Level: intermediate 332 333 References: 334 . * - http://dx.doi.org/10.1016/j.cam.2005.03.008 335 336 Options Database: 337 . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters 338 339 .seealso: `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetFilter()` 340 @*/ 341 PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt, PetscReal kkI, PetscReal kkP, PetscReal kkD) { 342 PetscFunctionBegin; 343 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 344 PetscValidLogicalCollectiveReal(adapt, kkI, 2); 345 PetscValidLogicalCollectiveReal(adapt, kkP, 3); 346 PetscValidLogicalCollectiveReal(adapt, kkD, 4); 347 PetscTryMethod(adapt, "TSAdaptDSPSetPID_C", (TSAdapt, PetscReal, PetscReal, PetscReal), (adapt, kkI, kkP, kkD)); 348 PetscFunctionReturn(0); 349 } 350 351 /*MC 352 TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP) 353 354 Level: intermediate 355 356 References: 357 + * - http://dx.doi.org/10.1145/641876.641877 358 - * - http://dx.doi.org/10.1016/j.cam.2005.03.008 359 360 Options Database: 361 + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers 362 . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters 363 . -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters 364 - -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters 365 366 .seealso: `TS`, `TSAdapt`, `TSGetAdapt()`, `TSAdaptDSPSetPID()`, `TSAdaptDSPSetFilter()` 367 M*/ 368 PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt) { 369 TSAdapt_DSP *dsp; 370 371 PetscFunctionBegin; 372 PetscCall(PetscNew(&dsp)); 373 adapt->reject_safety = 1.0; /* unused */ 374 375 adapt->data = (void *)dsp; 376 adapt->ops->choose = TSAdaptChoose_DSP; 377 adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP; 378 adapt->ops->destroy = TSAdaptDestroy_DSP; 379 adapt->ops->view = TSAdaptView_DSP; 380 381 PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetFilter_C", TSAdaptDSPSetFilter_DSP)); 382 PetscCall(PetscObjectComposeFunction((PetscObject)adapt, "TSAdaptDSPSetPID_C", TSAdaptDSPSetPID_DSP)); 383 384 PetscCall(TSAdaptDSPSetFilter_DSP(adapt, "PI42")); 385 PetscCall(TSAdaptRestart_DSP(adapt)); 386 PetscFunctionReturn(0); 387 } 388