<font size="2"><span style="font-family: Arial;">I want to learn petsc by
solving a small time-dependent problem:<br> The problem is from
chemical kinetics, and consists of the following three rate
equations: <br>
dy1/dt = -.04*y1 + 1.e4*y2*y3<br> dy2/dt = .04*y1 - 1.e4*y2*y3 -
3.e7*(y2)^2<br> dy3/dt = 3.e7*(y2)^2<br> on the interval from
t = 0.0 to t = 4.e10, with initial<br> conditions: y1 = 1.0, y2 = y3 =
0.<br><br>I passed the compiling, and generated an executable file. But I got a
lot of - Error Message - when I run the file. I appreciate your help and
suggestions! Many thanks!<br><br>My program is :<br><br>#include
"petscts.h"<br>/*The problem is from<br> * chemical kinetics, and consists
of the following three rate<br> *
equations:
<br> * dy1/dt = -.04*y1 +
1.e4*y2*y3<br> * dy2/dt = .04*y1 - 1.e4*y2*y3 -
3.e7*(y2)^2<br> * dy3/dt = 3.e7*(y2)^2<br> * on the
interval from t = 0.0 to t = 4.e10, with initial<br> * conditions: y1 =
1.0, y2 = y3 = 0.<br> */<br>typedef struct {<br>
PetscInt m;<br>}AppCtx;<br><br>extern PetscErrorCode
FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*),<br> FormFunction(TS,PetscReal,Vec,Vec,void*);<br><br>int main(int argc,char **argv)<br>{<br> TS ts;<br> Vec u,r;<br> PetscScalar *u_localptr;<br> Mat J;<br> AppCtx user;<br> PetscInt its,m;<br> PetscReal dt,ftime;<br> PetscErrorCode ierr;<br> PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);<br> m=3;<br> dt=0.1;<br> ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);<br> user.m=m;<br> u_localptr[0]=1.0;<br> u_localptr[1]=0.0;<br> u_localptr[2]=0.0;<br> ierr = VecRestoreArray(u,&u_localptr);CHKERRQ(ierr); <br> <br> ierr = VecCreateSeq(PETSC_COMM_SELF,m,&u);CHKERRQ(ierr);<br!
><br> ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr);<br> ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);<br> ierr = TSSetSolution(ts,u);CHKERRQ(ierr);<br> ierr = TSSetRHSFunction(ts,FormFunction,&user);CHKERRQ(ierr);<br> ierr = TSSetRHSJacobian(ts,J,J,FormJacobian,&user);CHKERRQ(ierr);<br><br> ierr = TSSetType(ts,TSPSEUDO);CHKERRQ(ierr);<br> ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);<br> ierr = TSPseudoSetTimeStep(ts,TSPseudoDefaultTimeStep,0);CHKERRQ(ierr);<br> ierr = TSSetFromOptions(ts);CHKERRQ(ierr);<br><br> ierr = TSSetUp(ts);CHKERRQ(ierr);<br> ierr = TSStep(ts,&its,&ftime);CHKERRQ(ierr);<br> <br> printf("Number of pseudo timesteps = %d final time %4.2e\n",(int)its,ftime);<br> ierr = VecDestroy(u);CHKERRQ(ierr);<br> ierr = VecDestroy(r);CHKERRQ(ierr);<br> ierr = MatDestroy(J);CHKERRQ(ierr);<br> ierr = TSDestroy(ts);CHKERRQ(ie!
rr);<br> ierr = PetscFinalize();CHKERRQ(ierr);<br><br>&n!
bsp; ret
urn 0;<br>}<br><br><br><br>PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)<br>{<br><br> PetscErrorCode ierr;<br> PetscScalar *x,*f;<br> PetscInt n;<br> n=3;<br> ierr = VecGetArray(X,&x);CHKERRQ(ierr);<br> ierr = VecGetArray(F,&f);CHKERRQ(ierr);<br> f[0]=-0.04*x[0]+1.0e4*x[1]*x[2];<br> f[1]=0.04-1.0e4*x[1]*x[2]-3*1.0e7*x[1]*x[1];<br> f[2]=3*1.0e7*x[1]*x[1];<br> ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);<br> ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);<br> return 0;<br> <br>}<br>PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)<br>{<br> <br> Mat jac=*B;<br> PetscScalar v[3],*x;<br> PetscInt r!
ow,col;<br> PetscErrorCode ierr;<br> ierr = VecGetArray(X,&x);CHKERRQ(ierr);<br> v[0]=-0.04;<br> v[1]=0.04;<br> v[2]=0;<br> row=1;<br> col=1;<br> ierr = MatSetValues(jac,1,&row,1,&col,v,INSERT_VALUES);CHKERRQ(ierr);<br> v[0]=1.0e4*x[2];<br> v[1]=-1.0e4*x[2]-6.0e7*x[1];<br> v[2]=6.0e7*x[1];<br> col=2;<br> ierr = MatSetValues(jac,1,&row,1,&col,v,INSERT_VALUES);CHKERRQ(ierr);<br> v[0]=1.0e4*x[1];<br> v[1]=-1.0e4*x[1];<br> v[2]=0;<br> col=3;<br> &n!
bsp; ierr = MatSetValues(jac,1,&row,1,&col,v,INSERT_VA!
LUES);CH
KERRQ(ierr);<br> ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br> ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);<br> ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br> return 0;<br><br>}<br><br></span></font><br>