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