1 2 #ifndef lint 3 static char vcid[] = "$Id: euler.c,v 1.1 1996/01/06 16:31:06 bsmith Exp bsmith $"; 4 #endif 5 /* 6 Code for Time Stepping with explicit Euler. 7 */ 8 #include <math.h> 9 #include "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 static int TSSetUp_Euler(TS ts) 18 { 19 TS_Euler *euler = (TS_Euler*) ts->data; 20 int ierr; 21 22 ierr = VecDuplicate(ts->vec_sol,&euler->update); CHKERRQ(ierr); 23 return 0; 24 } 25 26 static int TSStep_Euler(TS ts,int *steps,Scalar *time) 27 { 28 TS_Euler *euler = (TS_Euler*) ts->data; 29 Vec sol = ts->vec_sol,update = euler->update; 30 int ierr,i,max_steps = ts->max_steps; 31 32 *steps = -ts->steps; 33 ierr = (*ts->monitor)(ts,ts->steps,ts->ptime,sol,ts->monP); CHKERRQ(ierr); 34 35 for ( i=0; i<max_steps; i++ ) { 36 ierr = TSComputeRHSFunction(ts,ts->ptime,sol,update); CHKERRQ(ierr); 37 ierr = VecAXPY(&ts->time_step,update,sol); CHKERRQ(ierr); 38 ts->ptime += ts->time_step; 39 ts->steps++; 40 ierr = (*ts->monitor)(ts,ts->steps,ts->ptime,sol,ts->monP); CHKERRQ(ierr); 41 if (ts->ptime > ts->max_time) break; 42 } 43 44 *steps += ts->steps; 45 *time = ts->ptime; 46 return 0; 47 } 48 /*------------------------------------------------------------*/ 49 static int TSDestroy_Euler(PetscObject obj ) 50 { 51 TS ts = (TS) obj; 52 TS_Euler *euler = (TS_Euler*) ts->data; 53 54 VecDestroy(euler->update); 55 PetscFree(euler); 56 return 0; 57 } 58 /*------------------------------------------------------------*/ 59 60 static int TSSetFromOptions_Euler(TS ts) 61 { 62 63 return 0; 64 } 65 66 static int TSPrintHelp_Euler(TS ts) 67 { 68 69 return 0; 70 } 71 72 static int TSView_Euler(PetscObject obj,Viewer viewer) 73 { 74 return 0; 75 } 76 77 /* ------------------------------------------------------------ */ 78 int TSCreate_Euler(TS ts ) 79 { 80 TS_Euler *euler; 81 82 ts->type = 0; 83 ts->setup = TSSetUp_Euler; 84 ts->step = TSStep_Euler; 85 ts->destroy = TSDestroy_Euler; 86 ts->printhelp = TSPrintHelp_Euler; 87 ts->setfromoptions = TSSetFromOptions_Euler; 88 ts->view = TSView_Euler; 89 90 euler = PetscNew(TS_Euler); CHKPTRQ(euler); 91 ts->data = (void *) euler; 92 93 return 0; 94 } 95