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