1*b45d2f2cSJed Brown #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 284df9cb4SJed Brown 384df9cb4SJed Brown typedef struct { 41c3436cfSJed Brown PetscBool always_accept; 51c3436cfSJed Brown PetscReal clip[2]; /* admissible decrease/increase factors */ 61c3436cfSJed Brown PetscReal safety; /* safety factor relative to target error */ 77e1ba4dcSJed Brown PetscReal reject_safety; /* extra safety factor if the last step was rejected */ 81c3436cfSJed Brown Vec Y; 984df9cb4SJed Brown } TSAdapt_Basic; 1084df9cb4SJed Brown 1184df9cb4SJed Brown #undef __FUNCT__ 121c3436cfSJed Brown #define __FUNCT__ "TSAdaptChoose_Basic" 130b99f514SJed Brown static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte) 1484df9cb4SJed Brown { 151c3436cfSJed Brown TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 161c3436cfSJed Brown PetscErrorCode ierr; 171c3436cfSJed Brown Vec X,Y; 187e1ba4dcSJed Brown PetscReal enorm,hfac_lte,h_lte,safety; 191c3436cfSJed Brown PetscInt order,stepno; 2084df9cb4SJed Brown 2184df9cb4SJed Brown PetscFunctionBegin; 221c3436cfSJed Brown ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr); 231c3436cfSJed Brown ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 241c3436cfSJed Brown if (!basic->Y) {ierr = VecDuplicate(X,&basic->Y);CHKERRQ(ierr);} 251c3436cfSJed Brown Y = basic->Y; 261c3436cfSJed Brown order = adapt->candidates.order[0]; 271c3436cfSJed Brown ierr = TSEvaluateStep(ts,order-1,Y,PETSC_NULL);CHKERRQ(ierr); 281c3436cfSJed Brown 297e1ba4dcSJed Brown safety = basic->safety; 301c3436cfSJed Brown ierr = TSErrorNormWRMS(ts,Y,&enorm);CHKERRQ(ierr); 311c3436cfSJed Brown if (enorm > 1.) { 327e1ba4dcSJed Brown if (!*accept) safety *= basic->reject_safety; /* The last attempt also failed, shorten more aggressively */ 33fd94acc0SJed Brown if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) { 34fd94acc0SJed Brown ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, accepting because step size %G is at minimum\n",enorm,h);CHKERRQ(ierr); 35fd94acc0SJed Brown *accept = PETSC_TRUE; 367e1ba4dcSJed Brown } else if (basic->always_accept) { 377e1ba4dcSJed Brown ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, accepting step of size %G because always_accept is set\n",enorm,h);CHKERRQ(ierr); 387e1ba4dcSJed Brown *accept = PETSC_TRUE; 391c3436cfSJed Brown } else { 40fd94acc0SJed Brown ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, rejecting step of size %G\n",enorm,h);CHKERRQ(ierr); 41fd94acc0SJed Brown *accept = PETSC_FALSE; 42fd94acc0SJed Brown } 43fd94acc0SJed Brown } else { 44fd94acc0SJed Brown ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %G, accepting step of size %G\n",enorm,h);CHKERRQ(ierr); 45fd94acc0SJed Brown *accept = PETSC_TRUE; 461c3436cfSJed Brown } 471c3436cfSJed Brown 481c3436cfSJed Brown /* The optimal new step based purely on local truncation error for this step. */ 497e1ba4dcSJed Brown hfac_lte = safety * PetscRealPart(PetscPowScalar((PetscScalar)enorm,(PetscReal)(-1./order))); 501c3436cfSJed Brown h_lte = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]); 511c3436cfSJed Brown 5284df9cb4SJed Brown *next_sc = 0; 531c3436cfSJed Brown *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 540b99f514SJed Brown *wlte = enorm; 5584df9cb4SJed Brown PetscFunctionReturn(0); 5684df9cb4SJed Brown } 5784df9cb4SJed Brown 5884df9cb4SJed Brown #undef __FUNCT__ 5984df9cb4SJed Brown #define __FUNCT__ "TSAdaptDestroy_Basic" 6084df9cb4SJed Brown static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt) 6184df9cb4SJed Brown { 621c3436cfSJed Brown TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 6384df9cb4SJed Brown PetscErrorCode ierr; 6484df9cb4SJed Brown 6584df9cb4SJed Brown PetscFunctionBegin; 661c3436cfSJed Brown ierr = VecDestroy(&basic->Y);CHKERRQ(ierr); 6784df9cb4SJed Brown ierr = PetscFree(adapt->data);CHKERRQ(ierr); 6884df9cb4SJed Brown PetscFunctionReturn(0); 6984df9cb4SJed Brown } 7084df9cb4SJed Brown 711c3436cfSJed Brown #undef __FUNCT__ 721c3436cfSJed Brown #define __FUNCT__ "TSAdaptSetFromOptions_Basic" 731c3436cfSJed Brown static PetscErrorCode TSAdaptSetFromOptions_Basic(TSAdapt adapt) 741c3436cfSJed Brown { 751c3436cfSJed Brown TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 761c3436cfSJed Brown PetscErrorCode ierr; 771c3436cfSJed Brown PetscInt two; 781c3436cfSJed Brown PetscBool set; 791c3436cfSJed Brown 801c3436cfSJed Brown PetscFunctionBegin; 811c3436cfSJed Brown ierr = PetscOptionsHead("Basic adaptive controller options");CHKERRQ(ierr); 821c3436cfSJed Brown two = 2; 831c3436cfSJed Brown ierr = PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase in step size","",basic->clip,&two,&set);CHKERRQ(ierr); 841c3436cfSJed Brown if (set && (two != 2 || basic->clip[0] > basic->clip[1])) 851c3436cfSJed Brown SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip"); 861c3436cfSJed Brown ierr = PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,PETSC_NULL);CHKERRQ(ierr); 877e1ba4dcSJed Brown ierr = PetscOptionsReal("-ts_adapt_basic_reject_safety","Extra safety factor to apply if the last step was rejected","",basic->reject_safety,&basic->reject_safety,PETSC_NULL);CHKERRQ(ierr); 881c3436cfSJed Brown ierr = PetscOptionsBool("-ts_adapt_basic_always_accept","Always accept the step regardless of whether local truncation error meets goal","",basic->always_accept,&basic->always_accept,PETSC_NULL);CHKERRQ(ierr); 891c3436cfSJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 901c3436cfSJed Brown PetscFunctionReturn(0); 911c3436cfSJed Brown } 921c3436cfSJed Brown 9384df9cb4SJed Brown EXTERN_C_BEGIN 9484df9cb4SJed Brown #undef __FUNCT__ 9584df9cb4SJed Brown #define __FUNCT__ "TSAdaptCreate_Basic" 9684df9cb4SJed Brown /*MC 9784df9cb4SJed Brown TSADAPTBASIC - Basic adaptive controller for time stepping 9884df9cb4SJed Brown 9984df9cb4SJed Brown Level: intermediate 10084df9cb4SJed Brown 10184df9cb4SJed Brown .seealso: TS, TSAdapt, TSSetAdapt() 10284df9cb4SJed Brown M*/ 10384df9cb4SJed Brown PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 10484df9cb4SJed Brown { 10584df9cb4SJed Brown PetscErrorCode ierr; 10684df9cb4SJed Brown TSAdapt_Basic *a; 10784df9cb4SJed Brown 10884df9cb4SJed Brown PetscFunctionBegin; 10984df9cb4SJed Brown ierr = PetscNewLog(adapt,TSAdapt_Basic,&a);CHKERRQ(ierr); 11084df9cb4SJed Brown adapt->data = (void*)a; 1111c3436cfSJed Brown adapt->ops->choose = TSAdaptChoose_Basic; 1121c3436cfSJed Brown adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic; 11384df9cb4SJed Brown adapt->ops->destroy = TSAdaptDestroy_Basic; 1141c3436cfSJed Brown 1151c3436cfSJed Brown a->clip[0] = 0.1; 1161c3436cfSJed Brown a->clip[1] = 10.; 1171c3436cfSJed Brown a->safety = 0.9; 1187e1ba4dcSJed Brown a->reject_safety = 0.5; 1197e1ba4dcSJed Brown a->always_accept = PETSC_FALSE; 12084df9cb4SJed Brown PetscFunctionReturn(0); 12184df9cb4SJed Brown } 12284df9cb4SJed Brown EXTERN_C_END 123