[petsc-users] help about my first petsc program

[Xuan YU]

I got the correct results! Thanks!

Could you please tell me whether I can use TS to solve ODE without  
providing Jacobian matrix?

On Jun 22, 2010, at 7:19 PM, Jed Brown wrote:

> 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);
>

Xuan YU (俞烜)
xxy113 at psu.edu




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