[petsc-users] help about my first petsc program

[XUAN YU]

I got my program to run. But the results doesn't make any sence.

At t =          0  y =1.000000e+00  0.000000e+00  0.000000e+00, 
At t =       0.01  y =9.996000e-01  4.000000e-04  0.000000e+00, 
At t =       0.02  y =9.992002e-01  -4.720016e-02  4.800000e-02, 
At t =       0.03  y =7.722397e-01  -6.681768e+02  6.684045e+02, 
At t =       0.04  y =-4.466124e+07  -1.338934e+11  1.339381e+11, 
At t =       0.05  y =-1.793342e+24  -5.376439e+27  5.378233e+27, 
At t =       0.06  y =-2.891574e+57  -8.668938e+60  8.671830e+60, 
At t =       0.07  y =-7.517556e+123  -2.253763e+127  2.254515e+127, 
At t =       0.08  y =-5.081142e+256  -1.523326e+260  1.523834e+260, 
At t =       0.09  y =-inf  nan  inf, 
At t =        0.1  y =nan  nan  nan, 
At t =       0.11  y =nan  nan  nan, 
At t =       0.12  y =nan  nan  nan, 
At t =       0.13  y =nan  nan  nan, 
At t =       0.14  y =nan  nan  nan, 

Would you please help me check what's wrong with my code?

#include "petscts.h"

/*The problem is from
 * chemical kinetics, and consists of the following three rate
 * equations:         
 *    dy1/dt = -.04*y1 + 1.e4*y2*y3
 *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
 *    dy3/dt = 3.e7*(y2)^2
 * on the interval from t = 0.0 to t = 4.e10, with initial
 * conditions: y1 = 1.0, y2 = y3 = 0.*/


typedef struct {
    PetscInt m;
}AppCtx;

extern PetscErrorCode 
FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*),Monitor(TS
ts,PetscInt step,PetscReal time,Vec u,void*ctx),
     FormFunction(TS,PetscReal,Vec,Vec,void*);

int main(int argc,char **argv)
{
    TS ts;
    Vec u;
    PetscScalar    *u_localptr;
    Mat J;
    AppCtx user;
    PetscInt its,m;
    PetscReal dt,ftime;
    PetscErrorCode ierr;
    PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);
    m=3;
    dt=0.01;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);
  user.m=m;
  ierr=VecCreateSeq(PETSC_COMM_SELF,m,&u);CHKERRQ(ierr);
  ierr=VecGetArray(u,&u_localptr);CHKERRQ(ierr);
  u_localptr[0]=1.0;
  u_localptr[1]=0.0;
  u_localptr[2]=0.0;
  ierr = VecRestoreArray(u,&u_localptr);CHKERRQ(ierr);    
  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
  ierr = TSMonitorSet(ts,Monitor,&user,PETSC_NULL);
  ierr = TSSetSolution(ts,u);CHKERRQ(ierr);
  ierr = TSSetRHSFunction(ts,FormFunction,&user);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_SELF,&J);CHKERRQ(ierr);
  ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(ierr);
  ierr = MatSetFromOptions(J);
  ierr = TSSetRHSJacobian(ts,J,J,FormJacobian,&user);CHKERRQ(ierr);
  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
  ierr = TSSetDuration(ts,1000,8.e-1);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);

  ierr = TSSetUp(ts);CHKERRQ(ierr);
  ierr = TSStep(ts,&its,&ftime);CHKERRQ(ierr);
  
  printf("Number of timesteps = %d final time %4.2e\n",(int)its,ftime);
  ierr = VecDestroy(u);CHKERRQ(ierr);
  ierr = MatDestroy(J);CHKERRQ(ierr);
  ierr = TSDestroy(ts);CHKERRQ(ierr);
  ierr = PetscFinalize();CHKERRQ(ierr);

  return 0;
}



PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
{ 
  PetscErrorCode ierr;
  PetscScalar *x,*f;
  PetscInt n;
  n=3;
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
        f[0]=-0.04*x[0]+1.0e4*x[1]*x[2];
    f[1]=0.04*x[0]-1.0e4*x[1]*x[2]-3*1.0e7*x[1]*x[1];
    f[2]=3*1.0e7*x[1]*x[1];
    ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  return 0;
  
}
PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure
*flag,void *ptr)
{
        Mat jac=*J;
        PetscScalar v[3],*x;
        PetscInt row,col;
        PetscErrorCode ierr;
       ierr = VecGetArray(X,&x);CHKERRQ(ierr);
       v[0]=-0.04;
       v[1]=0.04;
       v[2]=0.0;
       row=0;
       col=0;
       ierr = MatSetValues(jac,3,&row,1,&col,v,INSERT_VALUES);CHKERRQ(ierr);
           v[0]=1.0e4*x[2];
       v[1]=-1.0e4*x[2]-6.0e7*x[1];
       v[2]=6.0e7*x[1];
       col=1;
       ierr = MatSetValues(jac,3,&row,1,&col,v,INSERT_VALUES);CHKERRQ(ierr);
       v[0]=1.0e4*x[1];
       v[1]=-1.0e4*x[1];
       v[2]=0.0;
       col=2;
       ierr = MatSetValues(jac,3,&row,1,&col,v,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  return 0;

}
PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u, void *ctx)
{
  PetscScalar *y;
  
  VecGetArray(u,&y);
  PetscPrintf(PETSC_COMM_WORLD,"At t =%11g  y =%e  %e  %e,
\n",time,y[0],y[1],y[2]);
  VecRestoreArray(u,&y);
  return 0;
}


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100622/1c621853/attachment.htm>