<div dir="ltr">On Wed, Oct 9, 2013 at 6:08 PM, Mark F. Adams <span dir="ltr"><<a href="mailto:mfadams@lbl.gov" target="_blank">mfadams@lbl.gov</a>></span> wrote:<br><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
><br>
> Your Fortran is hurting my eyes :)<br>
><br>
> Is solver%rVec2 == 0?<br>
<br>
Should be OK.  Here is a fuller listing appended.<br>
<br>
SNESSetUp seems wrong in that I don't see how SNESSetUp can work unless a user drills though the interface (or I'm missing some SMESSetXXX method or a side effect).  This code is in SNESSetUp:<br>
<br>
  ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr);<br>
  if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");<br>
  if (snes->vec_rhs  == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");<br>
<br>
  if (!snes->vec_sol_update /* && snes->vec_sol */) {<br>
    ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr);<br>
    ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr);<br>
  }<br>
<br>
<br>
It seems like first, the (snes->vec_rhs  == snes->vec_sol) test should be moved to SNESSolve where these are apparently set.<br></blockquote><div><br></div><div>This is not good enough since the next lines will fail if snes->vec_sol is NULL. I think we do not to substitute vec_func for vec_sol, but</div>
<div>we also need to move creation of the update vector since it must be the size of the solution. Will do it now.</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

Second can we use snes->vec_func instead of snes->vec_sol here?  This would remove all dependance on vec_rhs and vec_sol in SNESSetup.<br>
<br>
Mark<br>
<br>
  !  Extract global and local vectors from DM; then duplicate for remaining<br>
  !     vectors that are the same types<br>
  call MatGetVecs(solver%KKTmat,solver%xVec2,solver%bVec2,ierr)<br>
  call VecDuplicate(solver%bVec2,solver%rVec2,ierr)<br>
<br>
  call SNESSetDM(solver%snes,solver%da,ierr)<br>
<br>
  ! call SNESSetApplicationContext(solver%snes,solver,ierr)<br>
  ! call SNESSetApplicationContext(solver%snes,grid,ierr)<br>
<br>
  !  Set function evaluation routine and vector<br>
  call SNESSetFunction(solver%snes,solver%rVec2,FormFunction,solver,ierr)<br>
<br>
  call SNESSetJacobian(solver%snes,solver%KKTmat,solver%KKTmat,FormJacobian,solver,ierr)<br>
<br>
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
  !  Customize nonlinear solver; set runtime options<br>
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
  !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)<br>
  call SNESSetFromOptions(solver%snes,ierr)<br>
<br>
  call SNESSetUp(solver%snes,ierr) ! pre setup, same as (old) KSP<br>
<br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>