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 static PetscErrorCode TSStep_Euler(TS ts) 11 { 12 TS_Euler *euler = (TS_Euler*)ts->data; 13 Vec solution = ts->vec_sol,update = euler->update; 14 PetscBool stageok,accept = PETSC_TRUE; 15 PetscReal next_time_step = ts->time_step; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 ierr = TSPreStage(ts,ts->ptime);CHKERRQ(ierr); 20 ierr = TSComputeRHSFunction(ts,ts->ptime,solution,update);CHKERRQ(ierr); 21 ierr = VecAYPX(update,ts->time_step,solution);CHKERRQ(ierr); 22 ierr = TSPostStage(ts,ts->ptime,0,&solution);CHKERRQ(ierr); 23 ierr = TSAdaptCheckStage(ts->adapt,ts,ts->ptime,solution,&stageok);CHKERRQ(ierr); 24 if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; PetscFunctionReturn(0);} 25 ierr = TSFunctionDomainError(ts,ts->ptime+ts->time_step,update,&stageok);CHKERRQ(ierr); 26 if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; PetscFunctionReturn(0);} 27 28 ierr = TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 29 if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; PetscFunctionReturn(0);} 30 ierr = VecCopy(update,solution);CHKERRQ(ierr); 31 32 ts->ptime += ts->time_step; 33 ts->time_step = next_time_step; 34 PetscFunctionReturn(0); 35 } 36 /*------------------------------------------------------------*/ 37 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 51 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 52 ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode TSReset_Euler(TS ts) 57 { 58 TS_Euler *euler = (TS_Euler*)ts->data; 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = VecDestroy(&euler->update);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 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 static PetscErrorCode TSSetFromOptions_Euler(PetscOptionItems *PetscOptionsObject,TS ts) 78 { 79 PetscFunctionBegin; 80 PetscFunctionReturn(0); 81 } 82 83 static PetscErrorCode TSView_Euler(TS ts,PetscViewer viewer) 84 { 85 PetscFunctionBegin; 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode TSInterpolate_Euler(TS ts,PetscReal t,Vec X) 90 { 91 TS_Euler *euler = (TS_Euler*)ts->data; 92 Vec update = euler->update; 93 PetscReal alpha = (ts->ptime - t)/ts->time_step; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = VecWAXPY(X,-ts->time_step,update,ts->vec_sol);CHKERRQ(ierr); 98 ierr = VecAXPBY(X,1.0-alpha,alpha,ts->vec_sol);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode TSComputeLinearStability_Euler(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 103 { 104 PetscFunctionBegin; 105 *yr = 1.0 + xr; 106 *yi = xi; 107 PetscFunctionReturn(0); 108 } 109 /* ------------------------------------------------------------ */ 110 111 /*MC 112 TSEULER - ODE solver using the explicit forward Euler method 113 114 Level: beginner 115 116 .seealso: TSCreate(), TS, TSSetType(), TSBEULER 117 118 M*/ 119 PETSC_EXTERN PetscErrorCode TSCreate_Euler(TS ts) 120 { 121 TS_Euler *euler; 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 ierr = PetscNewLog(ts,&euler);CHKERRQ(ierr); 126 ts->data = (void*)euler; 127 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 PetscFunctionReturn(0); 137 } 138