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