On Sat, Feb 7, 2015 at 11:00 AM, Barry Smith <bsm...@mcs.anl.gov> wrote:
> > 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? > Clone g...@bitbucket.org:knepley/magma-dynamics-problems.git Build it, make magma.dbg Run some tests, make run_stokes Let me know if you have problems. I think I can get it to fail on a small PETSc example, but I will not be able to confirm before Tuesday. Margarete is in London and I have Maria until then. Thanks, Matt > 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 <knep...@gmail.com> wrote: > > > > On Fri, Feb 6, 2015 at 9:35 PM, Jed Brown <j...@jedbrown.org> wrote: > > Barry Smith <bsm...@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 > > -- 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