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