<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sat, Feb 7, 2015 at 11:00 AM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><br>
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?<br></blockquote><div><br></div><div>Clone git@bitbucket.org:knepley/magma-dynamics-problems.git</div><div><br></div><div>Build it, make magma.dbg</div><div><br></div><div>Run some tests, make run_stokes</div><div><br></div><div>Let me know if you have problems.</div><div><br></div><div>I think I can get it to fail on a small PETSc example, but I will not be</div><div>able to confirm before Tuesday. Margarete is in London and I have</div><div>Maria until then.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
PetscErrorCode SNESSetSolution(SNES snes, Vec u)<br>
{<br>
DM dm;<br>
PetscErrorCode ierr;<br>
<br>
PetscFunctionBegin;<br>
PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);<br>
PetscValidHeaderSpecific(u, VEC_CLASSID, 2);<br>
ierr = PetscObjectReference((PetscObject) u);CHKERRQ(ierr);<br>
ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr);<br>
<br>
snes->vec_sol = u;<br>
<br>
ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);<br>
ierr = DMShellSetGlobalVector(dm, u);CHKERRQ(ierr);<br>
PetscFunctionReturn(0);<br>
}<br>
<br>
PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y)<br>
{<br>
PetscErrorCode ierr;<br>
DM dm;<br>
DMSNES sdm;<br>
<br>
PetscFunctionBegin;<br>
PetscValidHeaderSpecific(snes,SNES_CLASSID,1);<br>
PetscValidHeaderSpecific(x,VEC_CLASSID,2);<br>
PetscValidHeaderSpecific(y,VEC_CLASSID,3);<br>
PetscCheckSameComm(snes,1,x,2);<br>
PetscCheckSameComm(snes,1,y,3);<br>
ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);<br>
<br>
ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);<br>
ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);<br>
if (sdm->ops->computefunction) {<br>
ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr);<br>
ierr = VecLockPush(x);CHKERRQ(ierr);<br>
PetscStackPush("SNES user function");<br>
ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr);<br>
PetscStackPop;<br>
ierr = VecLockPop(x);CHKERRQ(ierr);<br>
ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr);<br>
} else if (snes->vec_rhs) {<br>
ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr);<br>
} else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");<br>
if (snes->vec_rhs) {<br>
ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr);<br>
}<br>
snes->nfuncs++;<br>
ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr);<br>
PetscFunctionReturn(0);<br>
<div class=""><div class="h5">}<br>
<br>
> On Feb 7, 2015, at 12:53 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
><br>
> On Fri, Feb 6, 2015 at 9:35 PM, Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>> wrote:<br>
> Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> writes:<br>
> > I am not proposing removing the solution from TSSolve() argument<br>
> > list, far from it. I am proposing removing the TSSetSolution(),<br>
> > keeping the solution in TSSolve() as the standard approach, but<br>
> > also allowing a callback for when the user wants the DM passed to<br>
> > TS to be in charge of creating all the vectors.<br>
><br>
> Oh, that's fine if/when you have a use case for it.<br>
><br>
> I needed SNESSetSolution() so I put it in. Here is the use case.<br>
><br>
> I have a TS (magma dynamics), for which I would like to test the solver. I<br>
> have an explicit solution when I turn some terms off, and I test all these.<br>
> Sometimes I turn off time, so I TSGetSNES() and SNESSetSolution(), and<br>
> then call SNESComputeFunction(snes, exactSol, res). If I do not SetSolution()<br>
> first, things are not setup correctly and it bombs.<br>
><br>
> Matt<br>
><br>
> --<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<br>
<br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">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></div>