"Fossies" - the Fresh Open Source Software Archive

Member "petsc-3.11.3/src/ts/adapt/impls/dsp/adaptdsp.c" (12 Sep 2018, 14646 Bytes) of package /linux/misc/petsc-3.11.3.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "adaptdsp.c" see the Fossies "Dox" file reference documentation.

    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[0];
  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 controlers.
  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>
  333 
  334 .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetFilter()
  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>
  361 
  362 .seealso: TS, TSAdapt, TSGetAdapt(), TSAdaptDSPSetPID()
  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>
  385 .   -ts_adapt_dsp_pid <kkI,kkP,kkD>
  386 .   -ts_adapt_dsp_kbeta <b1,b2,b2>
  387 -   -ts_adapt_dsp_alpha <a2,a3>
  388 
  389 .seealso: TS, TSAdapt, TSGetAdapt()
  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 }