> On Feb 13, 2017, at 3:08 PM, Ed Bueler <[email protected]> wrote:
>
> Barry --
>
> > Sounds like a bug to me. The methods should be checking if an
> > IFunction is being provided and error out in that case.
>
> I don't think it is that simple. I.e. having an IFunction and an explicit
> scheme is not, by itself, a bug. I think.
>
> If the user has provided IFunction
> F(t,u,u_t) = u_t - f(t,u)
> which is the convenient form for e.g. TSARKIMEX,
> and RHSFunction
> G(t,u)
> then I guess I assumed that explicit methods like TSRK should,for
> unconstrained cases, get their RHS by callback like this:
> g(t,u) = - F(t,u,0) + G(t,u)
> so that the ODE is in form
> u_t = g(t,u) = f(t,u) + g(t,u)
I guess in theory TSRK could do this but looking at the code
static PetscErrorCode TSStep_RK(TS ts)
{
TS_RK *rk = (TS_RK*)ts->data;
RKTableau tab = rk->tableau;
const PetscInt s = tab->s;
const PetscReal *A = tab->A,*c = tab->c;
PetscScalar *w = rk->work;
Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS;
TSAdapt adapt;
PetscInt i,j;
PetscInt rejections = 0;
PetscBool stageok,accept = PETSC_TRUE;
PetscReal next_time_step = ts->time_step;
PetscErrorCode ierr;
PetscFunctionBegin;
rk->status = TS_STEP_INCOMPLETE;
while (!ts->reason && rk->status != TS_STEP_COMPLETE) {
PetscReal t = ts->ptime;
PetscReal h = ts->time_step;
for (i=0; i<s; i++) {
rk->stage_time = t + h*c[i];
ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr);
ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr);
for (j=0; j<i; j++) w[j] = h*A[i*s+j];
ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr);
ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr);
ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);
ierr =
TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr);
if (!stageok) goto reject_step;
ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr);
etscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
{
PetscErrorCode ierr;
TSRHSFunction rhsfunction;
TSIFunction ifunction;
void *ctx;
DM dm;
PetscFunctionBegin;
PetscValidHeaderSpecific(ts,TS_CLASSID,1);
PetscValidHeaderSpecific(U,VEC_CLASSID,3);
PetscValidHeaderSpecific(y,VEC_CLASSID,4);
ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
ierr = DMTSGetRHSFunction(dm,&rhsfunction,&ctx);CHKERRQ(ierr);
ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr);
if (!rhsfunction && !ifunction)
SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call
TSSetRHSFunction() and / or TSSetIFunction()");
ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);CHKERRQ(ierr);
if (rhsfunction) {
PetscStackPush("TS user right-hand-side function");
ierr = (*rhsfunction)(ts,t,U,y,ctx);CHKERRQ(ierr);
PetscStackPop;
} else {
ierr = VecZeroEntries(y);CHKERRQ(ierr);
}
ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);CHKERRQ(ie
it doesn't do that. It simply ignores the ifunction the user provided. We
should not do that since it will produce the wrong answer.
Regarding the "poor man's way" of handling the constraint with projection
methods. Our intention is that you would use TSSetPostStep() to provide a
function that takes the computed solution by TSRK and "projects" it as you
wish. You can try this now but note you will need to pack your ifunction()
(with the u_t) and your rhsfunction() into a single right hand side function*.
Barry
*Yes, eventually we should support the user not having to provide the function
in a different form but we don't have anybody working on "useful and needed
improvements to TS" we only have people working on "publishable and fundable
cool new stuff" with TS.
>
> I think that is the behavior I am seeing (i.e. on my problem, using PETSc
> master).
>
> The "nonsense" I am referring to is only from the non-enforcement of the
> constraint; it would be o.k. for a pure ODE.
>
> I would love to have some projection mechanism to try, for comparing explicit
> methods with projection to the SNESVI way (i.e. the right way), but that is
> asking for a lot of PETSc refactoring, I think. For now I just want to
> error-out if the method is not going to call the SNES at each time step.
>
> Ed
>
>
>
> On Mon, Feb 13, 2017 at 11:26 AM, Barry Smith <[email protected]> wrote:
>
> > On Feb 13, 2017, at 1:53 PM, Ed Bueler <[email protected]> wrote:
> >
> > Dear Petsc --
> >
> > This is the first of two short TS usage questions.
> >
> > My problem is both stiff (diffusive PDE) and constrained, so I require
> >
> > -snes_type vinewton{rs|ss}ls
> >
> > *and* I split my ODE system into IFunction and RHSFunction
> >
> > F(t,u,u_t) = G(t,u)
> >
> > where F(t,u,u_t) = u_t + f(t,u) in my case (i.e. no mass matrix needed),
> > and the stiff part goes in f(t,u).
> >
> > With this arrangement TS types beuler, theta, bdf, arkimex all work quite
> > well. However, the program runs and produces nonsense with type rk and
> > ssp, that is, explicit methods.
>
> Sounds like a bug to me. The methods should be checking if an IFunction is
> being provided and error out in that case.
>
> Barry
>
> >
> > So my question is, how do I ask the TS (at run time) whether the chosen TS
> > type will or will not call its SNES at each step? If SNES is not going to
> > be used then I want to SETERRQ and stop. That is, I want to error-out if
> > the *method* is fully explicit.
> >
> > Note the constraints are enforced by the SNESVI, as they should be, not ad
> > hoc projection. Also, as a technical matter, I cannot require my iterates
> > to be feasible inside my IFunction evaluation because that would break
> > VINEWTONSSLS.
> >
> > Neither TSProblemType (mine is NONLINEAR) nor TSEquationType (mine is
> > IMPLICIT I guess) seem to address this? My problem is indeed nonlinear and
> > has stiff parts, but it is not a DAE and it *is* constrained.
> >
> > Thanks!
> >
> > Ed
> >
> > PS I'd prefer not to enumerate the existing TS types and error on the bad
> > ones. It is not nicely-maintainable.
> >
> >
> >
> > --
> > Ed Bueler
> > Dept of Math and Stat and Geophysical Institute
> > University of Alaska Fairbanks
> > Fairbanks, AK 99775-6660
> > 301C Chapman and 410D Elvey
> > 907 474-7693 and 907 474-7199 (fax 907 474-5394)
>
>
>
>
> --
> Ed Bueler
> Dept of Math and Stat and Geophysical Institute
> University of Alaska Fairbanks
> Fairbanks, AK 99775-6660
> 301C Chapman and 410D Elvey
> 907 474-7693 and 907 474-7199 (fax 907 474-5394)
