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