1 /* 2 Code for Timestepping with explicit Euler. 3 */ 4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 5 6 typedef struct { 7 Vec update; /* work vector where new solution is formed */ 8 } TS_Euler; 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "TSStep_Euler" 12 static PetscErrorCode TSStep_Euler(TS ts) 13 { 14 TS_Euler *euler = (TS_Euler*)ts->data; 15 Vec sol = ts->vec_sol,update = euler->update; 16 PetscErrorCode ierr; 17 PetscBool accept; 18 19 PetscFunctionBegin; 20 ierr = TSPreStage(ts,ts->ptime);CHKERRQ(ierr); 21 ierr = TSComputeRHSFunction(ts,ts->ptime,sol,update);CHKERRQ(ierr); 22 ierr = VecAXPY(sol,ts->time_step,update);CHKERRQ(ierr); 23 ierr = TSFunctionDomainError(ts,ts->ptime,sol,&accept);CHKERRQ(ierr); 24 if(!accept) { 25 ts->reason = TS_DIVERGED_STEP_REJECTED; 26 PetscFunctionReturn(0); 27 } 28 ierr = TSPostStage(ts,ts->ptime,0,&sol);CHKERRQ(ierr); 29 ts->ptime += ts->time_step; 30 ts->steps++; 31 PetscFunctionReturn(0); 32 } 33 /*------------------------------------------------------------*/ 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "TSSetUp_Euler" 37 static PetscErrorCode TSSetUp_Euler(TS ts) 38 { 39 TS_Euler *euler = (TS_Euler*)ts->data; 40 PetscErrorCode ierr; 41 TSRHSFunction rhsfunction; 42 TSIFunction ifunction; 43 44 PetscFunctionBegin; 45 ierr = TSGetIFunction(ts,NULL,&ifunction,NULL);CHKERRQ(ierr); 46 ierr = TSGetRHSFunction(ts,NULL,&rhsfunction,NULL);CHKERRQ(ierr); 47 if (!rhsfunction || ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must define RHSFunction() and leave IFunction() undefined in order to use -ts_type euler"); 48 ierr = VecDuplicate(ts->vec_sol,&euler->update);CHKERRQ(ierr); 49 PetscFunctionReturn(0); 50 } 51 52 #undef __FUNCT__ 53 #define __FUNCT__ "TSReset_Euler" 54 static PetscErrorCode TSReset_Euler(TS ts) 55 { 56 TS_Euler *euler = (TS_Euler*)ts->data; 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 ierr = VecDestroy(&euler->update);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "TSDestroy_Euler" 66 static PetscErrorCode TSDestroy_Euler(TS ts) 67 { 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 ierr = TSReset_Euler(ts);CHKERRQ(ierr); 72 ierr = PetscFree(ts->data);CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 /*------------------------------------------------------------*/ 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "TSSetFromOptions_Euler" 79 static PetscErrorCode TSSetFromOptions_Euler(PetscOptionItems *PetscOptionsObject,TS ts) 80 { 81 PetscFunctionBegin; 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "TSView_Euler" 87 static PetscErrorCode TSView_Euler(TS ts,PetscViewer viewer) 88 { 89 PetscFunctionBegin; 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "TSInterpolate_Euler" 95 static PetscErrorCode TSInterpolate_Euler(TS ts,PetscReal t,Vec X) 96 { 97 TS_Euler *euler = (TS_Euler*)ts->data; 98 Vec update = euler->update; 99 PetscReal alpha = (ts->ptime - t)/ts->time_step; 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 ierr = VecWAXPY(X,-ts->time_step,update,ts->vec_sol);CHKERRQ(ierr); 104 ierr = VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 #undef __FUNCT__ 109 #define __FUNCT__ "TSComputeLinearStability_Euler" 110 PetscErrorCode TSComputeLinearStability_Euler(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 111 { 112 PetscFunctionBegin; 113 *yr = 1.0 + xr; 114 *yi = xi; 115 PetscFunctionReturn(0); 116 } 117 /* ------------------------------------------------------------ */ 118 119 /*MC 120 TSEULER - ODE solver using the explicit forward Euler method 121 122 Level: beginner 123 124 .seealso: TSCreate(), TS, TSSetType(), TSBEULER 125 126 M*/ 127 #undef __FUNCT__ 128 #define __FUNCT__ "TSCreate_Euler" 129 PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts) 130 { 131 TS_Euler *euler; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 ts->ops->setup = TSSetUp_Euler; 136 ts->ops->step = TSStep_Euler; 137 ts->ops->reset = TSReset_Euler; 138 ts->ops->destroy = TSDestroy_Euler; 139 ts->ops->setfromoptions = TSSetFromOptions_Euler; 140 ts->ops->view = TSView_Euler; 141 ts->ops->interpolate = TSInterpolate_Euler; 142 ts->ops->linearstability = TSComputeLinearStability_Euler; 143 144 /* does not have adaptivity so delete adaptivity object */ 145 ierr = TSAdaptDestroy(&ts->adapt);CHKERRQ(ierr); 146 147 ierr = PetscNewLog(ts,&euler);CHKERRQ(ierr); 148 ts->data = (void*)euler; 149 PetscFunctionReturn(0); 150 } 151