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)
