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 PetscCheckFalse(enorm >= 0 && order < 1,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 PetscCheckFalse(adapt->candidates.n < 1,PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered"); 90 PetscCheckFalse(!adapt->candidates.inuse_set,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 PetscCheckFalse(!tab,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 PetscCheckFalse(set && !n,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 PetscCheckFalse(set && !n,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 PetscCheckFalse(set && !n,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: 315 . * - http://dx.doi.org/10.1145/641876.641877 316 317 Notes: 318 Valid filter names are 319 + "basic" - similar to TSADAPTBASIC but with different criteria for step rejections. 320 . "PI30", "PI42", "PI33", "PI34" - PI controllers. 321 . "PC11", "PC47", "PC36" - predictive controllers. 322 . "H0211", "H211b", "H211PI" - digital filters with orders dynamics=2, adaptivity=1, filter=1. 323 . "H0312", "H312b", "H312PID" - digital filters with orders dynamics=3, adaptivity=1, filter=2. 324 - "H0321", "H321" - digital filters with orders dynamics=3, adaptivity=2, filter=1. 325 326 Options Database: 327 . -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers 328 329 .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID() 330 @*/ 331 PetscErrorCode TSAdaptDSPSetFilter(TSAdapt adapt,const char *name) 332 { 333 PetscErrorCode ierr; 334 PetscFunctionBegin; 335 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 336 PetscValidCharPointer(name,2); 337 ierr = PetscTryMethod(adapt,"TSAdaptDSPSetFilter_C",(TSAdapt,const char*),(adapt,name));CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 /*@ 342 TSAdaptDSPSetPID - Set the PID controller parameters 343 344 Input Parameters: 345 + adapt - adaptive controller context 346 . kkI - Integral parameter 347 . kkP - Proportional parameter 348 - kkD - Derivative parameter 349 350 Level: intermediate 351 352 References: 353 . * - http://dx.doi.org/10.1016/j.cam.2005.03.008 354 355 Options Database: 356 . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters 357 358 .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetFilter() 359 @*/ 360 PetscErrorCode TSAdaptDSPSetPID(TSAdapt adapt,PetscReal kkI,PetscReal kkP,PetscReal kkD) 361 { 362 PetscErrorCode ierr; 363 PetscFunctionBegin; 364 PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); 365 PetscValidLogicalCollectiveReal(adapt,kkI,2); 366 PetscValidLogicalCollectiveReal(adapt,kkP,3); 367 PetscValidLogicalCollectiveReal(adapt,kkD,4); 368 ierr = PetscTryMethod(adapt,"TSAdaptDSPSetPID_C",(TSAdapt,PetscReal,PetscReal,PetscReal),(adapt,kkI,kkP,kkD));CHKERRQ(ierr); 369 PetscFunctionReturn(0); 370 } 371 372 /*MC 373 TSADAPTDSP - Adaptive controller for time-stepping based on digital signal processing (DSP) 374 375 Level: intermediate 376 377 References: 378 + * - http://dx.doi.org/10.1145/641876.641877 379 - * - http://dx.doi.org/10.1016/j.cam.2005.03.008 380 381 Options Database: 382 + -ts_adapt_dsp_filter <name> - Sets predefined controller by name; use -help for a list of available controllers 383 . -ts_adapt_dsp_pid <kkI,kkP,kkD> - Sets PID controller parameters 384 . -ts_adapt_dsp_kbeta <b1,b2,b2> - Sets general filter parameters 385 - -ts_adapt_dsp_alpha <a2,a3> - Sets general filter parameters 386 387 .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID(), TSAdaptDSPSetFilter() 388 M*/ 389 PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt adapt) 390 { 391 TSAdapt_DSP *dsp; 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 ierr = PetscNewLog(adapt,&dsp);CHKERRQ(ierr); 396 adapt->reject_safety = 1.0; /* unused */ 397 398 adapt->data = (void*)dsp; 399 adapt->ops->choose = TSAdaptChoose_DSP; 400 adapt->ops->setfromoptions = TSAdaptSetFromOptions_DSP; 401 adapt->ops->destroy = TSAdaptDestroy_DSP; 402 adapt->ops->view = TSAdaptView_DSP; 403 404 ierr = PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetFilter_C",TSAdaptDSPSetFilter_DSP);CHKERRQ(ierr); 405 ierr = PetscObjectComposeFunction((PetscObject)adapt,"TSAdaptDSPSetPID_C",TSAdaptDSPSetPID_DSP);CHKERRQ(ierr); 406 407 ierr = TSAdaptDSPSetFilter_DSP(adapt,"PI42");CHKERRQ(ierr); 408 ierr = TSAdaptRestart_DSP(adapt);CHKERRQ(ierr); 409 PetscFunctionReturn(0); 410 } 411