> 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]> wrote: > > > On Feb 14, 2017, at 8:56 PM, Emil Constantinescu <[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)
