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