>
> Your Fortran is hurting my eyes :)
>
> Is solver%rVec2 == 0?
Should be OK. Here is a fuller listing appended.
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:
ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr);
if (snes->vec_func == snes->vec_sol)
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function
vector");
if (snes->vec_rhs == snes->vec_sol)
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand
side vector");
if (!snes->vec_sol_update /* && snes->vec_sol */) {
ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr);
ierr =
PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr);
}
It seems like first, the (snes->vec_rhs == snes->vec_sol) test should be moved
to SNESSolve where these are apparently set.
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.
Mark
! Extract global and local vectors from DM; then duplicate for remaining
! vectors that are the same types
call MatGetVecs(solver%KKTmat,solver%xVec2,solver%bVec2,ierr)
call VecDuplicate(solver%bVec2,solver%rVec2,ierr)
call SNESSetDM(solver%snes,solver%da,ierr)
! call SNESSetApplicationContext(solver%snes,solver,ierr)
! call SNESSetApplicationContext(solver%snes,grid,ierr)
! Set function evaluation routine and vector
call SNESSetFunction(solver%snes,solver%rVec2,FormFunction,solver,ierr)
call
SNESSetJacobian(solver%snes,solver%KKTmat,solver%KKTmat,FormJacobian,solver,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Customize nonlinear solver; set runtime options
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type
<type>)
call SNESSetFromOptions(solver%snes,ierr)
call SNESSetUp(solver%snes,ierr) ! pre setup, same as (old) KSP
