1*8b1af7b3SBarry Smith 24b33d51aSBarry Smith #ifndef lint 3*8b1af7b3SBarry Smith static char vcid[] = "$Id: euler.c,v 1.1 1996/01/06 16:31:06 bsmith Exp bsmith $"; 44b33d51aSBarry Smith #endif 54b33d51aSBarry Smith /* 6*8b1af7b3SBarry Smith Code for Time Stepping with explicit Euler. 74b33d51aSBarry Smith */ 84b33d51aSBarry Smith #include <math.h> 94b33d51aSBarry Smith #include "tsimpl.h" /*I "ts.h" I*/ 104b33d51aSBarry Smith #include "pinclude/pviewer.h" 114b33d51aSBarry Smith 124b33d51aSBarry Smith 13*8b1af7b3SBarry Smith typedef struct { 14*8b1af7b3SBarry Smith Vec update; /* work vector where F(t[i],u[i]) is stored */ 15*8b1af7b3SBarry Smith } TS_Euler; 16*8b1af7b3SBarry Smith 17*8b1af7b3SBarry Smith static int TSSetUp_Euler(TS ts) 184b33d51aSBarry Smith { 19*8b1af7b3SBarry Smith TS_Euler *euler = (TS_Euler*) ts->data; 204b33d51aSBarry Smith int ierr; 214b33d51aSBarry Smith 22*8b1af7b3SBarry Smith ierr = VecDuplicate(ts->vec_sol,&euler->update); CHKERRQ(ierr); 234b33d51aSBarry Smith return 0; 244b33d51aSBarry Smith } 254b33d51aSBarry Smith 26*8b1af7b3SBarry Smith static int TSStep_Euler(TS ts,int *steps,Scalar *time) 274b33d51aSBarry Smith { 28*8b1af7b3SBarry Smith TS_Euler *euler = (TS_Euler*) ts->data; 29*8b1af7b3SBarry Smith Vec sol = ts->vec_sol,update = euler->update; 30*8b1af7b3SBarry Smith int ierr,i,max_steps = ts->max_steps; 314b33d51aSBarry Smith 32*8b1af7b3SBarry Smith *steps = -ts->steps; 33*8b1af7b3SBarry Smith ierr = (*ts->monitor)(ts,ts->steps,ts->ptime,sol,ts->monP); CHKERRQ(ierr); 344b33d51aSBarry Smith 35*8b1af7b3SBarry Smith for ( i=0; i<max_steps; i++ ) { 36*8b1af7b3SBarry Smith ierr = TSComputeRHSFunction(ts,ts->ptime,sol,update); CHKERRQ(ierr); 37*8b1af7b3SBarry Smith ierr = VecAXPY(&ts->time_step,update,sol); CHKERRQ(ierr); 38*8b1af7b3SBarry Smith ts->ptime += ts->time_step; 39*8b1af7b3SBarry Smith ts->steps++; 40*8b1af7b3SBarry Smith ierr = (*ts->monitor)(ts,ts->steps,ts->ptime,sol,ts->monP); CHKERRQ(ierr); 41*8b1af7b3SBarry Smith if (ts->ptime > ts->max_time) break; 42*8b1af7b3SBarry Smith } 434b33d51aSBarry Smith 44*8b1af7b3SBarry Smith *steps += ts->steps; 45*8b1af7b3SBarry Smith *time = ts->ptime; 464b33d51aSBarry Smith return 0; 474b33d51aSBarry Smith } 484b33d51aSBarry Smith /*------------------------------------------------------------*/ 49*8b1af7b3SBarry Smith static int TSDestroy_Euler(PetscObject obj ) 504b33d51aSBarry Smith { 514b33d51aSBarry Smith TS ts = (TS) obj; 52*8b1af7b3SBarry Smith TS_Euler *euler = (TS_Euler*) ts->data; 53*8b1af7b3SBarry Smith 54*8b1af7b3SBarry Smith VecDestroy(euler->update); 55*8b1af7b3SBarry Smith PetscFree(euler); 56*8b1af7b3SBarry Smith return 0; 57*8b1af7b3SBarry Smith } 58*8b1af7b3SBarry Smith /*------------------------------------------------------------*/ 59*8b1af7b3SBarry Smith 60*8b1af7b3SBarry Smith static int TSSetFromOptions_Euler(TS ts) 61*8b1af7b3SBarry Smith { 624b33d51aSBarry Smith 634b33d51aSBarry Smith return 0; 644b33d51aSBarry Smith } 654b33d51aSBarry Smith 66*8b1af7b3SBarry Smith static int TSPrintHelp_Euler(TS ts) 674b33d51aSBarry Smith { 68*8b1af7b3SBarry Smith 69*8b1af7b3SBarry Smith return 0; 70*8b1af7b3SBarry Smith } 71*8b1af7b3SBarry Smith 72*8b1af7b3SBarry Smith static int TSView_Euler(PetscObject obj,Viewer viewer) 73*8b1af7b3SBarry Smith { 74*8b1af7b3SBarry Smith return 0; 75*8b1af7b3SBarry Smith } 76*8b1af7b3SBarry Smith 77*8b1af7b3SBarry Smith /* ------------------------------------------------------------ */ 78*8b1af7b3SBarry Smith int TSCreate_Euler(TS ts ) 79*8b1af7b3SBarry Smith { 80*8b1af7b3SBarry Smith TS_Euler *euler; 81*8b1af7b3SBarry Smith 824b33d51aSBarry Smith ts->type = 0; 83*8b1af7b3SBarry Smith ts->setup = TSSetUp_Euler; 84*8b1af7b3SBarry Smith ts->step = TSStep_Euler; 85*8b1af7b3SBarry Smith ts->destroy = TSDestroy_Euler; 86*8b1af7b3SBarry Smith ts->printhelp = TSPrintHelp_Euler; 87*8b1af7b3SBarry Smith ts->setfromoptions = TSSetFromOptions_Euler; 88*8b1af7b3SBarry Smith ts->view = TSView_Euler; 89*8b1af7b3SBarry Smith 90*8b1af7b3SBarry Smith euler = PetscNew(TS_Euler); CHKPTRQ(euler); 91*8b1af7b3SBarry Smith ts->data = (void *) euler; 924b33d51aSBarry Smith 934b33d51aSBarry Smith return 0; 944b33d51aSBarry Smith } 95