[petsc-dev] TSSolve web page JED and EMIL READ THIS!

[Barry Smith]

  Hmm, we should definitely debug this because SNESComputeFunction() should be usable as a stand-alone function. In a quick look at the code I cannot see why a SNESSetSolution() call should be needed, nor would I want it to be needed in order to use SNESComputeFunction().  How can I reproduce your problem?

PetscErrorCode SNESSetSolution(SNES snes, Vec u)
{
  DM             dm;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
  PetscValidHeaderSpecific(u, VEC_CLASSID, 2);
  ierr = PetscObjectReference((PetscObject) u);CHKERRQ(ierr);
  ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr);

  snes->vec_sol = u;

  ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
  ierr = DMShellSetGlobalVector(dm, u);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode  SNESComputeFunction(SNES snes,Vec x,Vec y)
{
  PetscErrorCode ierr;
  DM             dm;
  DMSNES         sdm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
  PetscValidHeaderSpecific(x,VEC_CLASSID,2);
  PetscValidHeaderSpecific(y,VEC_CLASSID,3);
  PetscCheckSameComm(snes,1,x,2);
  PetscCheckSameComm(snes,1,y,3);
  ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);

  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
  ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);
  if (sdm->ops->computefunction) {
    ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr);
    ierr = VecLockPush(x);CHKERRQ(ierr);
    PetscStackPush("SNES user function");
    ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr);
    PetscStackPop;
    ierr = VecLockPop(x);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr);
  } else if (snes->vec_rhs) {
    ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr);
  } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
  if (snes->vec_rhs) {
    ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr);
  }
  snes->nfuncs++;
  ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

> On Feb 7, 2015, at 12:53 AM, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Fri, Feb 6, 2015 at 9:35 PM, Jed Brown <jed at jedbrown.org> wrote:
> Barry Smith <bsmith at mcs.anl.gov> writes:
> >    I am not proposing removing the solution from TSSolve() argument
> >    list, far from it. I am proposing removing the TSSetSolution(),
> >    keeping the solution in TSSolve() as the standard approach, but
> >    also allowing a callback for when the user wants the DM passed to
> >    TS to be in charge of creating all the vectors.
> 
> Oh, that's fine if/when you have a use case for it.
> 
> I needed SNESSetSolution() so I put it in. Here is the use case.
> 
> I have a TS (magma dynamics), for which I would like to test the solver. I
> have an explicit solution when I turn some terms off, and I test all these.
> Sometimes I turn off time, so I TSGetSNES() and SNESSetSolution(), and
> then call SNESComputeFunction(snes, exactSol, res). If I do not SetSolution()
> first, things are not setup correctly and it bombs.
> 
>    Matt
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener