184df9cb4SJed Brown #include <private/tsimpl.h> /*I "petscts.h" I*/ 284df9cb4SJed Brown 384df9cb4SJed Brown typedef struct { 4*1c3436cfSJed Brown PetscBool always_accept; 5*1c3436cfSJed Brown PetscReal clip[2]; /* admissible decrease/increase factors */ 6*1c3436cfSJed Brown PetscReal safety; /* safety factor relative to target error */ 7*1c3436cfSJed Brown Vec Y; 884df9cb4SJed Brown } TSAdapt_Basic; 984df9cb4SJed Brown 1084df9cb4SJed Brown #undef __FUNCT__ 11*1c3436cfSJed Brown #define __FUNCT__ "TSAdaptChoose_Basic" 12*1c3436cfSJed Brown static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept) 1384df9cb4SJed Brown { 14*1c3436cfSJed Brown TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 15*1c3436cfSJed Brown PetscErrorCode ierr; 16*1c3436cfSJed Brown Vec X,Y; 17*1c3436cfSJed Brown PetscReal enorm,hfac_lte,h_lte; 18*1c3436cfSJed Brown PetscInt order,stepno; 1984df9cb4SJed Brown 2084df9cb4SJed Brown PetscFunctionBegin; 21*1c3436cfSJed Brown ierr = TSGetTimeStepNumber(ts,&stepno);CHKERRQ(ierr); 22*1c3436cfSJed Brown ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 23*1c3436cfSJed Brown if (!basic->Y) {ierr = VecDuplicate(X,&basic->Y);CHKERRQ(ierr);} 24*1c3436cfSJed Brown Y = basic->Y; 25*1c3436cfSJed Brown order = adapt->candidates.order[0]; 26*1c3436cfSJed Brown ierr = TSEvaluateStep(ts,order-1,Y,PETSC_NULL);CHKERRQ(ierr); 27*1c3436cfSJed Brown 28*1c3436cfSJed Brown ierr = TSErrorNormWRMS(ts,Y,&enorm);CHKERRQ(ierr); 29*1c3436cfSJed Brown if (enorm > 1.) { 30*1c3436cfSJed Brown ierr = PetscInfo1(adapt,"Estimated scaled local truncation error %G, step should be rejected\n",enorm);CHKERRQ(ierr); 31*1c3436cfSJed Brown } else { 32*1c3436cfSJed Brown ierr = PetscInfo1(adapt,"Estimated scaled local truncation error %G, step accepted\n",enorm);CHKERRQ(ierr); 33*1c3436cfSJed Brown } 34*1c3436cfSJed Brown 35*1c3436cfSJed Brown /* The optimal new step based purely on local truncation error for this step. */ 36*1c3436cfSJed Brown hfac_lte = basic->safety * PetscRealPart(PetscPowScalar(enorm,-1./order)); 37*1c3436cfSJed Brown h_lte = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]); 38*1c3436cfSJed Brown 3984df9cb4SJed Brown *next_sc = 0; 40*1c3436cfSJed Brown *next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max); 4184df9cb4SJed Brown *accept = PETSC_TRUE; 4284df9cb4SJed Brown PetscFunctionReturn(0); 4384df9cb4SJed Brown } 4484df9cb4SJed Brown 4584df9cb4SJed Brown #undef __FUNCT__ 4684df9cb4SJed Brown #define __FUNCT__ "TSAdaptDestroy_Basic" 4784df9cb4SJed Brown static PetscErrorCode TSAdaptDestroy_Basic(TSAdapt adapt) 4884df9cb4SJed Brown { 49*1c3436cfSJed Brown TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 5084df9cb4SJed Brown PetscErrorCode ierr; 5184df9cb4SJed Brown 5284df9cb4SJed Brown PetscFunctionBegin; 53*1c3436cfSJed Brown ierr = VecDestroy(&basic->Y);CHKERRQ(ierr); 5484df9cb4SJed Brown ierr = PetscFree(adapt->data);CHKERRQ(ierr); 5584df9cb4SJed Brown PetscFunctionReturn(0); 5684df9cb4SJed Brown } 5784df9cb4SJed Brown 58*1c3436cfSJed Brown #undef __FUNCT__ 59*1c3436cfSJed Brown #define __FUNCT__ "TSAdaptSetFromOptions_Basic" 60*1c3436cfSJed Brown static PetscErrorCode TSAdaptSetFromOptions_Basic(TSAdapt adapt) 61*1c3436cfSJed Brown { 62*1c3436cfSJed Brown TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data; 63*1c3436cfSJed Brown PetscErrorCode ierr; 64*1c3436cfSJed Brown PetscInt two; 65*1c3436cfSJed Brown PetscBool set; 66*1c3436cfSJed Brown 67*1c3436cfSJed Brown PetscFunctionBegin; 68*1c3436cfSJed Brown ierr = PetscOptionsHead("Basic adaptive controller options");CHKERRQ(ierr); 69*1c3436cfSJed Brown two = 2; 70*1c3436cfSJed Brown ierr = PetscOptionsRealArray("-ts_adapt_basic_clip","Admissible decrease/increase in step size","",basic->clip,&two,&set);CHKERRQ(ierr); 71*1c3436cfSJed Brown if (set && (two != 2 || basic->clip[0] > basic->clip[1])) 72*1c3436cfSJed Brown SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_basic_clip"); 73*1c3436cfSJed Brown ierr = PetscOptionsReal("-ts_adapt_basic_safety","Safety factor relative to target error","",basic->safety,&basic->safety,PETSC_NULL);CHKERRQ(ierr); 74*1c3436cfSJed 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); 75*1c3436cfSJed Brown if (!basic->always_accept) SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_SUP,"step rejection not implemented yet"); 76*1c3436cfSJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 77*1c3436cfSJed Brown PetscFunctionReturn(0); 78*1c3436cfSJed Brown } 79*1c3436cfSJed Brown 8084df9cb4SJed Brown EXTERN_C_BEGIN 8184df9cb4SJed Brown #undef __FUNCT__ 8284df9cb4SJed Brown #define __FUNCT__ "TSAdaptCreate_Basic" 8384df9cb4SJed Brown /*MC 8484df9cb4SJed Brown TSADAPTBASIC - Basic adaptive controller for time stepping 8584df9cb4SJed Brown 8684df9cb4SJed Brown Level: intermediate 8784df9cb4SJed Brown 8884df9cb4SJed Brown .seealso: TS, TSAdapt, TSSetAdapt() 8984df9cb4SJed Brown M*/ 9084df9cb4SJed Brown PetscErrorCode TSAdaptCreate_Basic(TSAdapt adapt) 9184df9cb4SJed Brown { 9284df9cb4SJed Brown PetscErrorCode ierr; 9384df9cb4SJed Brown TSAdapt_Basic *a; 9484df9cb4SJed Brown 9584df9cb4SJed Brown PetscFunctionBegin; 9684df9cb4SJed Brown ierr = PetscNewLog(adapt,TSAdapt_Basic,&a);CHKERRQ(ierr); 9784df9cb4SJed Brown adapt->data = (void*)a; 98*1c3436cfSJed Brown adapt->ops->choose = TSAdaptChoose_Basic; 99*1c3436cfSJed Brown adapt->ops->setfromoptions = TSAdaptSetFromOptions_Basic; 10084df9cb4SJed Brown adapt->ops->destroy = TSAdaptDestroy_Basic; 101*1c3436cfSJed Brown 102*1c3436cfSJed Brown a->clip[0] = 0.1; 103*1c3436cfSJed Brown a->clip[1] = 10.; 104*1c3436cfSJed Brown a->safety = 0.9; 105*1c3436cfSJed Brown a->always_accept = PETSC_TRUE; /* fix this */ 10684df9cb4SJed Brown PetscFunctionReturn(0); 10784df9cb4SJed Brown } 10884df9cb4SJed Brown EXTERN_C_END 109