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