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