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