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