[petsc-users] help about my first petsc program

[Jed Brown]

On Tue, 22 Jun 2010 19:09:32 -0400, "XUAN YU" <xxy113 at psu.edu> wrote:
> 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, 

There is a reason people don't use forward Euler for stiff systems.  Run
with -ts_type beuler or -ts_type theta, after applying the following
patch (which fixes your assembly bug).

Jed

diff --git a/ex1.c b/ex1.c
index ed78cd7..e4bc7f2 100644
--- a/ex1.c
+++ b/ex1.c
@@ -88,25 +88,24 @@ PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure
 {
   Mat jac=*J;
   PetscScalar v[3],*x;
-  PetscInt row,col;
+  PetscInt row[3] = {0,1,2},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);
+  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);
+  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 = 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);