Fwd: [petsc-maint #38493] TSSetIJacobian()

[Barry Smith]

    Any ideas?


Begin forwarded message:

> From: Dennis Zeki Kaya <Dennis.Kaya.3319 at student.uu.se>
> Date: December 4, 2009 5:06:35 AM CST
> To: petsc-maint Maintenance <petsc-maint at mcs.anl.gov>
> Subject: Re: [petsc-maint #38493] TSSetIJacobian()
> Reply-To: petsc-maint at mcs.anl.gov, Dennis Zeki Kaya <Dennis.Kaya.3319 at student.uu.se 
> >
>
> Hello,
>
> I'm using the TSTHETA solver to solve my DAE problem since I have a
> nonsingular Jacobian. I'm following the procedure in ex8.c in
> petsc-dev/src/ts/examples/tutorials as close as I can. I have attached
> a outline of my Petsc_Solve() routine to this email. The only big
> difference in my implementation to that used in ex8.c is that I create
> a global Jacobian matrix with my routine Petsc_JacCreate().
>
>
>
> Quoting Barry Smith <bsmith at mcs.anl.gov>:
>
>>
>>   Did you set the type to Backward Euler? The default Euler doesn't
>> use Jacobians.
>>
>>   Barry
>>
>> On Dec 3, 2009, at 6:28 PM, Dennis Zeki Kaya wrote:
>>
>>> Hello,
>>>
>>> I have a problem with a routine that evaluates a Jacobian matrix.  
>>> The
>>> problem is that PETSc never enters the RHSJacobian() routine. I know
>>> that PETSc enters the RHSFunction() rotine threw a print statement.
>>> The output is something like:
>>>
>>> At t[0] =      0.00e+00 u=       0.00e+00 at the center
>>> RHSFunction!
>>> RHSFunction!
>>> RHSFunction!
>>> ...
>>> RHSFunction!
>>> RHSFunction!
>>> RHSFunction!
>>> At t[10] =      1.00e-02 u=       0.00e+00 at the center
>>> RHSFunction!
>>> RHSFunction!
>>> ...
>>>
>>> I have taken bits and pieces for monitoring the solution procedure
>>> from the examples provided. Since the Jacobian is never computed the
>>> solution stays zero at each time-step, indifferent of parameters  
>>> like
>>> the size of time-step etc.
>>> This is the calls to the setting of the functions:
>>>
>>> /* Set user-provided RHSFunction and RHSJacobian */
>>> ierr = TSSetIFunction(ts, RHSFunction, ctx);CHKERRQ(ierr);
>>> ierr = TSSetIJacobian(ts, J, J, RHSJacobian, ctx);CHKERRQ(ierr);
>>>
>>> I can attach the whole Petsc_Solve() routine with the RHSFunction()
>>> and the RHSJacobian() routines if needed. I just thought that there
>>> maybe is an easy answer to my question regarding the time-stepper  
>>> not
>>> entering the RHSJacobian() routine.
>>>
>>> Thank you in advance,
>>> Dennis Kaya
>>>
>>>
>>>
>>>
>>>
>
>
>
> PetscErrorCode Petsc_Solve(/* Parameters for setting the DAE  
> userdefined context  */) {
>
>  PetscErrorCode ierr;
>  TS             ts; /* nonlinear solver */
>  PetscInt       time_steps=100, steps, iout, NOUT=1;
>  PetscReal      ftime;
>  Mat            J; /* Jacobian matrix  */
>  Vec            u; /* solution */
>  PC             pc; /* preconditioner */
>  DataCtx        *ctx; /* user-defined DAE-data context  */
>  PetscReal      dt;
>  PetscTruth     flg;
>  const TSType   tstype;
>  PetscViewer    viewer;
>  PetscViewer	 viewfile;
>  char           pcinfo[120], tsinfo[120];
>
>  SNES           ts_snes;
>  KSP            ksp;
>
>  PetscFunctionBegin;
>
>  /* Create Jacobian matrix  */
>  ierr = Petsc_JacCreate(&J, ctx);CHKERRQ(ierr);
>
>  /* Set the solution vector to have the same parallel layout as the  
> Jacobian matrix. */
>  ierr = MatGetVecs(J, &u, PETSC_NULL);CHKERRQ(ierr);
>
>  /* Create time-step context */
>  ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
>  ierr = TSSetProblemType(ts, TS_NONLINEAR);CHKERRQ(ierr);
>  ierr = TSSetType(ts, TSTHETA);CHKERRQ(ierr);
>
>  /* Set user-provided RHSFunction and RHSJacobian */
>  ierr = TSSetIFunction(ts, RHSFunction, ctx);CHKERRQ(ierr);
>  ierr = TSSetIJacobian(ts, J, J, RHSJacobian, ctx);CHKERRQ(ierr);
>  ierr = TSSetDuration(ts, time_steps, 0.3);CHKERRQ(ierr);
>
>  ierr = TSMonitorSet(ts, Monitor, ctx, PETSC_NULL);CHKERRQ(ierr);
>
>  /* Set up the preconditioner type to be used */
>  ierr = TSGetSNES(ts, &ts_snes);CHKERRQ(ierr);
>  ierr = SNESGetKSP(ts_snes, &ksp);CHKERRQ(ierr);
>  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
>  /* Set preconditioner to Jacobi, i.e diagonal scaling  
> preconditioning. */
>  ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr);
>
>  dt   = 0.001;
>  /* Set the initial solution vector to zero. */
>  ierr = Initial(0.0, u, ctx);CHKERRQ(ierr);
>  ierr = TSSetInitialTimeStep(ts, 0.0, dt);CHKERRQ(ierr);
>  ierr = TSSetSolution(ts, u);CHKERRQ(ierr);
>
>  /* Test TSSetPostStep() */
>  ierr = PetscOptionsHasName(PETSC_NULL,"- 
> test_PostStep",&flg);CHKERRQ(ierr);
>  if (flg) {
>    ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr);
>  }
>
>  /* Set up the time-stepping context from the the options given */
>  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
>
>  ierr = PetscOptionsGetInt(PETSC_NULL,"- 
> NOUT",&NOUT,PETSC_NULL);CHKERRQ(ierr);
>  for (iout=1; iout<=NOUT; iout++){
>    ierr = TSSetDuration(ts,time_steps,iout*1.0/NOUT);CHKERRQ(ierr);
>    ierr = TSStep(ts,&steps,&ftime);CHKERRQ(ierr);
>    ierr = TSSetInitialTimeStep(ts,ftime,dt);CHKERRQ(ierr);
>  }
>
>  ierr = PetscPrintf(PETSC_COMM_WORLD,"steps %D, ftime %G 
> \n",steps,ftime);CHKERRQ(ierr);
>
>  ierr = TSGetType(ts, &tstype);CHKERRQ(ierr);
>  ierr = PetscViewerStringOpen(PETSC_COMM_WORLD, tsinfo, 120,  
> &viewer);CHKERRQ(ierr);
>  ierr = TSView(ts, viewer);CHKERRQ(ierr);
>  ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
>  ierr = PetscViewerStringOpen(PETSC_COMM_WORLD, pcinfo, 120,  
> &viewer);CHKERRQ(ierr);
>  ierr = PCView(pc, viewer);CHKERRQ(ierr);
>  ierr = PetscPrintf(PETSC_COMM_WORLD, "%d Procs,%s TSType, %s  
> Preconditioner\n",
> 		     size, tsinfo, pcinfo);CHKERRQ(ierr);
>  ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
>
>  ierr = PetscOptionsHasName(PETSC_NULL,"- 
> matlab_view",&flg);CHKERRQ(ierr);
>  if (flg) { /* print solution into a MATLAB file */
>    ierr = TSGetSolution(ts, &u);CHKERRQ(ierr);
>    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",  
> &viewfile);CHKERRQ(ierr);
>    ierr = PetscViewerSetFormat(viewfile,  
> PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
>    ierr = VecView(u, viewfile);CHKERRQ(ierr);
>    ierr = PetscViewerDestroy(viewfile);CHKERRQ(ierr);
>  }
>
>  /* Free memory used by this routine. */
>  ierr = TSDestroy(ts);CHKERRQ(ierr);
>  ierr = VecDestroy(u);CHKERRQ(ierr);
>  ierr = PetscFree(ctx);CHKERRQ(ierr);
>
>  PetscFunctionReturn(0);
> }/*.Petsc_Solve.*/

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20091204/e2568379/attachment.html>