> 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
>>
