1 2 #ifndef lint 3 static char vcid[] = "$Id: euler.c,v 1.4 1996/08/08 14:45:37 bsmith Exp curfman $"; 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 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,double *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 Scalar dt = ts->time_step; 32 33 *steps = -ts->steps; 34 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol); CHKERRQ(ierr); 35 36 for ( i=0; i<max_steps; i++ ) { 37 ts->ptime += ts->time_step; 38 ierr = TSComputeRHSFunction(ts,ts->ptime,sol,update); CHKERRQ(ierr); 39 ierr = VecAXPY(&dt,update,sol); CHKERRQ(ierr); 40 ts->steps++; 41 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol); CHKERRQ(ierr); 42 if (ts->ptime > ts->max_time) break; 43 } 44 45 *steps += ts->steps; 46 *time = ts->ptime; 47 return 0; 48 } 49 /*------------------------------------------------------------*/ 50 static int TSDestroy_Euler(PetscObject obj ) 51 { 52 TS ts = (TS) obj; 53 TS_Euler *euler = (TS_Euler*) ts->data; 54 55 VecDestroy(euler->update); 56 PetscFree(euler); 57 return 0; 58 } 59 /*------------------------------------------------------------*/ 60 61 static int TSSetFromOptions_Euler(TS ts) 62 { 63 64 return 0; 65 } 66 67 static int TSPrintHelp_Euler(TS ts) 68 { 69 70 return 0; 71 } 72 73 static int TSView_Euler(PetscObject obj,Viewer viewer) 74 { 75 return 0; 76 } 77 78 /* ------------------------------------------------------------ */ 79 int TSCreate_Euler(TS ts ) 80 { 81 TS_Euler *euler; 82 83 ts->type = TS_EULER; 84 ts->setup = TSSetUp_Euler; 85 ts->step = TSStep_Euler; 86 ts->destroy = TSDestroy_Euler; 87 ts->printhelp = TSPrintHelp_Euler; 88 ts->setfromoptions = TSSetFromOptions_Euler; 89 ts->view = TSView_Euler; 90 91 euler = PetscNew(TS_Euler); CHKPTRQ(euler); 92 ts->data = (void *) euler; 93 94 return 0; 95 } 96 97 98 99 100 101