[petsc-users] TS question 1: how to stop explicit methods because they do not use SNES(VI)?
Ed Bueler
elbueler at alaska.edu
Tue Feb 14 23:44:22 CST 2017
> Jed suggested having any number of "RHS" functions. I don't think we
need more than 2, 1 for left hand side
> and 1 for right. If that ends up being not enough we can change to have
any number of them. Just to be clear.
> I suggest three functions
> IFunction which defaults to u_t plus TSSetMassMatrix() which changes
the default IFunction to M u_t
> LHS function which defaults to 0, if provided defaults to being treated
implicitly by IMEX
> RHS function which defaults to 0, if provided defaults to being treated
explicitly by IMEX
> Then a TSSetStiffMatrix(ts,Mat L) (horrible name) that provides u_t -Lu =
g() + Lu
At my current level of understanding I get this trifecta (quad-fecta?).
In my ice sheet case u_t = f(t,u) + g(t,u), I would not set IFunction at
all, or I would use a constant. The divergence of flux f(t,u) = div(q)
would go into the LHSFunction and the elevation-dependent mass balance
g(t,u) would go into the RHSFunction. I would not use TSSetStiffMatrix()
at all.
Barry's suggestion makes sense for this and the few other usage cases I
have already experienced ... kind of a lame endorsement, but, on such a
weak basis one casts a vote ...
Ed
PS I would not use TSSetStiffMatrix() because none of the
linearizations/simplfications known at the completion of a TS step could
"capture" the stiff part without being derived a la the Jacobian anyway.
To see what I mean, note that VINEWTONRSLS has a different *dimension* at
each Newton iteration. I don't know the (in-active-constraint) dimension
of the solution I'll get at each time step, much less how to scale the
various stiff-equation-parts for those in-active dimensions. In these
problems there is a lot going on in the implicit solve, besides just
handling stiffness.
On Tue, Feb 14, 2017 at 6:20 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> > On Feb 14, 2017, at 8:56 PM, Emil Constantinescu <emconsta at mcs.anl.gov>
> wrote:
> >
> > On 2/14/17 4:10 PM, Barry Smith wrote:
> >> Ok, you don't recompile but forcing that into user code is still
> disgusting. With my api the user code is
> >>
> >>>>> TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
> >>>>> TSSetLHSFunction(ts,NULL,LHSFunction,&ptype[0]);
> >>>>> TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
> >>>>> TSSetLHSJacobian(ts,Jac,Jac,LHSJacobian,&ptype[0]);
> >> and -ts_type xxx works correctly for ALL methods, implicit, explicit
> and imex without requiring any special command line options for different
> methods.
> >
> > Is this a viable solution? Growing the API to fix this situation will
> just put a burden with each new TS method after we refactor it in the
> current landscape.
>
> No just the opposite, the TS implementations will talk to functions who
> will put things together for it. So All implicit methods will call
> something like TSBuildImplicitFunction(), all explicit methods will call
> something like TSBuildExplicitFunction() and then IMEX methods will call
> both of these. In fact likely we can refactor to make things a little
> better than today. Depending on options and flagsTSBuildExplicitFunction()
> would build out of all the user provided functions what it needs etc.
>
> One problem with the current code is the TS methods call things with the
> same names as the user API. So implicit methods call TSComputeIFunction()
> while explicit methods call TSComputeRHSFunction(). This is wrong because
> implicit methods actually also use the rhs function provided by the user.
>
> The function below absolutely should not be called TSComputeIFunction()!
> It does not just compute IFunction()
> PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec
> Y,PetscBool imex)
> {
> PetscErrorCode ierr;
> TSIFunction ifunction;
> TSRHSFunction rhsfunction;
> void *ctx;
> DM dm;
>
> PetscFunctionBegin;
> PetscValidHeaderSpecific(ts,TS_CLASSID,1);
> PetscValidHeaderSpecific(U,VEC_CLASSID,3);
> PetscValidHeaderSpecific(Udot,VEC_CLASSID,4);
> PetscValidHeaderSpecific(Y,VEC_CLASSID,5);
>
> ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
> ierr = DMTSGetIFunction(dm,&ifunction,&ctx);CHKERRQ(ierr);
> ierr = DMTSGetRHSFunction(dm,&rhsfunction,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,Udot,Y);CHKERRQ(ierr);
> if (ifunction) {
> PetscStackPush("TS user implicit function");
> ierr = (*ifunction)(ts,t,U,Udot,Y,ctx);CHKERRQ(ierr);
> PetscStackPop;
> }
> if (imex) {
> if (!ifunction) {
> ierr = VecCopy(Udot,Y);CHKERRQ(ierr);
> }
> } else if (rhsfunction) {
> if (ifunction) {
> Vec Frhs;
> ierr = TSGetRHSVec_Private(ts,&Frhs);CHKERRQ(ierr);
> ierr = TSComputeRHSFunction(ts,t,U,Frhs);CHKERRQ(ierr);
> ierr = VecAXPY(Y,-1,Frhs);CHKERRQ(ierr);
> } else {
> ierr = TSComputeRHSFunction(ts,t,U,Y);CHKERRQ(ierr);
> ierr = VecAYPX(Y,-1,Udot);CHKERRQ(ierr);
> }
> }
> ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr);
> PetscFunctionReturn(0);
> }
>
> The current code entangles too much of the user API to the methods, this
> can be fixed.
>
> > If the user experiments with different ways of splitting the solution
> they would have to define RHS and IF or RHS and LHS in different ways
> (according to the splittings they experiment with). It may look disgusting,
> but I don't see another way around it unless you allow for a list of
> operators to be defined and then the user to assign them to LHS or RHS.
>
> Jed suggested having any number of "RHS" functions. I don't think we
> need more than 2, 1 for left hand side and 1 for right. If that ends up
> being not enough we can change to have any number of them. Just to be
> clear. I suggest three functions
>
> IFunction which defaults to u_t plus TSSetMassMatrix() which changes
> the default IFunction to M u_t
> LHS function which defaults to 0, if provided defaults to being treated
> implicitly by IMEX
> RHS function which defaults to 0, if provided defaults to being treated
> explicitly by IMEX
>
> Then a TSSetStiffMatrix(ts,Mat L) (horrible name) that provides u_t -Lu
> = g() + Lu
>
> None of the TS implementations will every directly know about what the
> user has provided. They will call the wrapper functions I mention above.
>
> I think Jed and Emil may be too enamored with the reductionist model of
> only IFunction() and RHSFunction() to see that though it encompasses
> everything it may not be the best user API.
>
>
>
> >
> >>
> >>> Adding all that logic to keep track of left sides and right sides for
> academic examples is likely not the best development.
> >> I don't think it is "just academic examples", it is all examples
> without a mass matrix.
> >>
> >> Once the user has decided with ts_type to use for production if it is
> fully implicit or explicit then they can depending on the type selected,
> write just a left hand side, just a right hand side for higher efficiency
> (less update of ghost points, fewer iterations over loops etc).
> >>
> >> With a constant mass matrix we can have TSSetMassMatrix() and then
> TSSetIFunction() is reserved for when it is absolutely needed.
> >
> > As much as I would disagree with growing the API at the level of
> defining the problem, I think TSSetMassMatrix() would let us do more things
> in the solvers. Also it would be useful to know if the mass matrix is
> singular or not for efficiency reasons.
> >
> > Emil
> >
> >> Barry
> >>
>
>
--
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)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170214/164ae87d/attachment-0001.html>
More information about the petsc-users
mailing list