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,0,&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 43 PetscFunctionBegin; 44 ierr = VecDuplicate(ts->vec_sol,&euler->update);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "TSReset_Euler" 50 static PetscErrorCode TSReset_Euler(TS ts) 51 { 52 TS_Euler *euler = (TS_Euler*)ts->data; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 ierr = VecDestroy(&euler->update);CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "TSDestroy_Euler" 62 static PetscErrorCode TSDestroy_Euler(TS ts) 63 { 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 ierr = TSReset_Euler(ts);CHKERRQ(ierr); 68 ierr = PetscFree(ts->data);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 /*------------------------------------------------------------*/ 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "TSSetFromOptions_Euler" 75 static PetscErrorCode TSSetFromOptions_Euler(PetscOptions *PetscOptionsObject,TS ts) 76 { 77 PetscFunctionBegin; 78 PetscFunctionReturn(0); 79 } 80 81 #undef __FUNCT__ 82 #define __FUNCT__ "TSView_Euler" 83 static PetscErrorCode TSView_Euler(TS ts,PetscViewer viewer) 84 { 85 PetscFunctionBegin; 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNCT__ 90 #define __FUNCT__ "TSInterpolate_Euler" 91 static PetscErrorCode TSInterpolate_Euler(TS ts,PetscReal t,Vec X) 92 { 93 PetscReal alpha = (ts->ptime - t)/ts->time_step; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = VecAXPBY(ts->vec_sol,1.0-alpha,alpha,X);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "TSComputeLinearStability_Euler" 103 PetscErrorCode TSComputeLinearStability_Euler(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 104 { 105 PetscFunctionBegin; 106 *yr = 1.0 + xr; 107 *yi = xi; 108 PetscFunctionReturn(0); 109 } 110 /* ------------------------------------------------------------ */ 111 112 /*MC 113 TSEULER - ODE solver using the explicit forward Euler method 114 115 Level: beginner 116 117 .seealso: TSCreate(), TS, TSSetType(), TSBEULER 118 119 M*/ 120 #undef __FUNCT__ 121 #define __FUNCT__ "TSCreate_Euler" 122 PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts) 123 { 124 TS_Euler *euler; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ts->ops->setup = TSSetUp_Euler; 129 ts->ops->step = TSStep_Euler; 130 ts->ops->reset = TSReset_Euler; 131 ts->ops->destroy = TSDestroy_Euler; 132 ts->ops->setfromoptions = TSSetFromOptions_Euler; 133 ts->ops->view = TSView_Euler; 134 ts->ops->interpolate = TSInterpolate_Euler; 135 ts->ops->linearstability = TSComputeLinearStability_Euler; 136 137 ierr = PetscNewLog(ts,&euler);CHKERRQ(ierr); 138 ts->data = (void*)euler; 139 PetscFunctionReturn(0); 140 } 141