1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: euler.c,v 1.12 1998/03/06 00:17:32 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Code for Timestepping with explicit Euler. 6 */ 7 #include <math.h> 8 #include "src/ts/tsimpl.h" /*I "ts.h" I*/ 9 #include "pinclude/pviewer.h" 10 11 12 typedef struct { 13 Vec update; /* work vector where F(t[i],u[i]) is stored */ 14 } TS_Euler; 15 16 #undef __FUNC__ 17 #define __FUNC__ "TSSetUp_Euler" 18 static int TSSetUp_Euler(TS ts) 19 { 20 TS_Euler *euler = (TS_Euler*) ts->data; 21 int ierr; 22 23 PetscFunctionBegin; 24 ierr = VecDuplicate(ts->vec_sol,&euler->update); CHKERRQ(ierr); 25 PetscFunctionReturn(0); 26 } 27 28 #undef __FUNC__ 29 #define __FUNC__ "TSStep_Euler" 30 static int TSStep_Euler(TS ts,int *steps,double *time) 31 { 32 TS_Euler *euler = (TS_Euler*) ts->data; 33 Vec sol = ts->vec_sol,update = euler->update; 34 int ierr,i,max_steps = ts->max_steps; 35 Scalar dt = ts->time_step; 36 37 PetscFunctionBegin; 38 *steps = -ts->steps; 39 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol); CHKERRQ(ierr); 40 41 for ( i=0; i<max_steps; i++ ) { 42 ts->ptime += ts->time_step; 43 ierr = TSComputeRHSFunction(ts,ts->ptime,sol,update); CHKERRQ(ierr); 44 ierr = VecAXPY(&dt,update,sol); CHKERRQ(ierr); 45 ts->steps++; 46 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol); CHKERRQ(ierr); 47 if (ts->ptime > ts->max_time) break; 48 } 49 50 *steps += ts->steps; 51 *time = ts->ptime; 52 PetscFunctionReturn(0); 53 } 54 /*------------------------------------------------------------*/ 55 #undef __FUNC__ 56 #define __FUNC__ "TSDestroy_Euler" 57 static int TSDestroy_Euler(TS ts) 58 { 59 TS_Euler *euler = (TS_Euler*) ts->data; 60 61 PetscFunctionBegin; 62 VecDestroy(euler->update); 63 PetscFree(euler); 64 PetscFunctionReturn(0); 65 } 66 /*------------------------------------------------------------*/ 67 68 #undef __FUNC__ 69 #define __FUNC__ "TSSetFromOptions_Euler" 70 static int TSSetFromOptions_Euler(TS ts) 71 { 72 PetscFunctionBegin; 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNC__ 77 #define __FUNC__ "TSPrintHelp_Euler" 78 static int TSPrintHelp_Euler(TS ts,char *p) 79 { 80 PetscFunctionBegin; 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNC__ 85 #define __FUNC__ "TSView_Euler" 86 static int TSView_Euler(TS ts,Viewer viewer) 87 { 88 PetscFunctionBegin; 89 PetscFunctionReturn(0); 90 } 91 92 /* ------------------------------------------------------------ */ 93 #undef __FUNC__ 94 #define __FUNC__ "TSCreate_Euler" 95 int TSCreate_Euler(TS ts ) 96 { 97 TS_Euler *euler; 98 99 PetscFunctionBegin; 100 ts->setup = TSSetUp_Euler; 101 ts->step = TSStep_Euler; 102 ts->destroy = TSDestroy_Euler; 103 ts->printhelp = TSPrintHelp_Euler; 104 ts->setfromoptions = TSSetFromOptions_Euler; 105 ts->view = TSView_Euler; 106 107 euler = PetscNew(TS_Euler); CHKPTRQ(euler); 108 PLogObjectMemory(ts,sizeof(TS_Euler)); 109 ts->data = (void *) euler; 110 111 PetscFunctionReturn(0); 112 } 113 114 115 116 117 118