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