1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: euler.c,v 1.9 1997/06/05 12:56:16 bsmith Exp balay $"; 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 ierr = VecDuplicate(ts->vec_sol,&euler->update); CHKERRQ(ierr); 24 return 0; 25 } 26 27 #undef __FUNC__ 28 #define __FUNC__ "TSStep_Euler" 29 static int TSStep_Euler(TS ts,int *steps,double *time) 30 { 31 TS_Euler *euler = (TS_Euler*) ts->data; 32 Vec sol = ts->vec_sol,update = euler->update; 33 int ierr,i,max_steps = ts->max_steps; 34 Scalar dt = ts->time_step; 35 36 *steps = -ts->steps; 37 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol); CHKERRQ(ierr); 38 39 for ( i=0; i<max_steps; i++ ) { 40 ts->ptime += ts->time_step; 41 ierr = TSComputeRHSFunction(ts,ts->ptime,sol,update); CHKERRQ(ierr); 42 ierr = VecAXPY(&dt,update,sol); CHKERRQ(ierr); 43 ts->steps++; 44 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol); CHKERRQ(ierr); 45 if (ts->ptime > ts->max_time) break; 46 } 47 48 *steps += ts->steps; 49 *time = ts->ptime; 50 return 0; 51 } 52 /*------------------------------------------------------------*/ 53 #undef __FUNC__ 54 #define __FUNC__ "TSDestroy_Euler" 55 static int TSDestroy_Euler(PetscObject obj ) 56 { 57 TS ts = (TS) obj; 58 TS_Euler *euler = (TS_Euler*) ts->data; 59 60 VecDestroy(euler->update); 61 PetscFree(euler); 62 return 0; 63 } 64 /*------------------------------------------------------------*/ 65 66 #undef __FUNC__ 67 #define __FUNC__ "TSSetFromOptions_Euler" 68 static int TSSetFromOptions_Euler(TS ts) 69 { 70 71 return 0; 72 } 73 74 #undef __FUNC__ 75 #define __FUNC__ "TSPrintHelp_Euler" 76 static int TSPrintHelp_Euler(TS ts,char *p) 77 { 78 79 return 0; 80 } 81 82 #undef __FUNC__ 83 #define __FUNC__ "TSView_Euler" 84 static int TSView_Euler(PetscObject obj,Viewer viewer) 85 { 86 return 0; 87 } 88 89 /* ------------------------------------------------------------ */ 90 #undef __FUNC__ 91 #define __FUNC__ "TSCreate_Euler" 92 int TSCreate_Euler(TS ts ) 93 { 94 TS_Euler *euler; 95 96 ts->type = TS_EULER; 97 ts->setup = TSSetUp_Euler; 98 ts->step = TSStep_Euler; 99 ts->destroy = TSDestroy_Euler; 100 ts->printhelp = TSPrintHelp_Euler; 101 ts->setfromoptions = TSSetFromOptions_Euler; 102 ts->view = TSView_Euler; 103 104 euler = PetscNew(TS_Euler); CHKPTRQ(euler); 105 PLogObjectMemory(ts,sizeof(TS_Euler)); 106 ts->data = (void *) euler; 107 108 return 0; 109 } 110 111 112 113 114 115