> On Feb 14, 2017, at 11:59 PM, Emil Constantinescu <[email protected]> > wrote: > > We seem to be needing more than two components and ways to pack them. > Allowing them to be dynamically assigned at run time would be really cool. > > We seem to like an API that allows us to specify M and u_t; but would also > like to keep support for f(t,u,u_t)=0.
Absolutely! No one is saying that IFunction should be removed! > > How about a costly experiment: build an API as discussed below (support the > above) that would map everything to the current one 1xRHS and 1xIFunction to > see how it would work? In other words, is there a scenario in which we could > introduce the new but also keep the old and be reasonable about it? Of > concern is whether the new one can be made friendly enough. > > Emil > > On 2/14/17 11:44 PM, Ed Bueler wrote: >>> 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 <[email protected] >> <mailto:[email protected]>> wrote: >> >> >> > On Feb 14, 2017, at 8:56 PM, Emil Constantinescu <[email protected] >> <mailto:[email protected]>> 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)
