[petsc-dev] Proposed changes to TS API

[Zhang, Hong]

On May 11, 2018, at 1:01 PM, Lisandro Dalcin <dalcinl at gmail.com<mailto:dalcinl at gmail.com>> wrote:

On Fri, 11 May 2018 at 19:34, Jed Brown <jed at jedbrown.org<mailto:jed at jedbrown.org>> wrote:

"Smith, Barry F." <bsmith at mcs.anl.gov<mailto:bsmith at mcs.anl.gov>> writes:


I assemble the combined system directly.

 How, two sets of calls to MatSetValues(), One for the scaled mass
matrix and one for the Jacobian entries? For a constant mass matrix does
this mean you are recomputing the mass matrix entries each call? Or are you
storing the mass matrix entries somewhere?  Or is your mass matrix diagonal
only?

 Or do you build element by element the M*shift + J element stiffness
and then insert it with a single MatSetValues() call?

It isn't even built separately at the element scale, just summed
contributions at quadrature points.

That's exactly the way I do it in PetIGA, and the way I would do it in any
other general-purpose FEM-like code. In high-order FEM, Jacobian assembly
may very well account from 50% to 80% of runtime (at least that's my
experience with PetIGA). IMHO, forcing users to have to do TWO global
matrix assemblies per time step is simply unacceptable, both in terms of
memory and runtime.

We are not forcing users to do two matrix assemblies per time step. For most cases, there is even no need to update dF/dUdot at all. For extreme cases that the application requires frequent update on dF/dUdot and assembly is expensive, one can always assemble the single-matrix Jacobian and throw it to SNES directly.


TS used to be an unusable pile of crap until Jed proposed the marvelous
IJacobian interface. Suddenly COMPLEX fully-implicit DAE problems become
SIMPLE to express, and a single IFunction/IJacobian pair reusable for
different timestepper implementations a reality. And we got an added
bounus: this was efficient, it involved a SINGLE global matrix assembly. If
the issue is in supporting simpler problems, then we should go for an
alternative interface with broken Jacobians, just as Stefano propossed in a
previous email. I'm totally in favor of an additional broken Jacobians API,
and totally againt the removal of the single-matrix IJacobian interface.

The current IJacobian is essentially SNESJacobian. And the single-matrix SNESJacobian interface is always there. Power users could set up the SNESJacobian directly if we pass a read-only shift parameter to them. So we are by no means prohibiting power users from doing what they want.

IJacobian with shift mixes TS Jacobian and SNES Jacobian. This is the issue we need to fix.

Thanks,
Hong

I don't buy the argument that IJacobian with shift is ugly, and that such
API drives users away from TS. At best, this is a documentation problem.
Come on, this is just chain rule, should be kindergarden-level stuff for
PETSc users. If we simplify things for the sake of making things simple for
newcomers and beginners and make them annoyingly slow for power users
solving complex problems, we will do a very bad business. This is not
politically correct, but I'm much worried about loosing power users, you
know, those that can eventually make a meaningful contributions to science
and software projects. Beginners and newcomers eventually learn and benefit
for common-sense software design, and will eventually appreciate it. I
really hope populism to not win this battle :-)
--
Lisandro Dalcin
============
Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa<http://www.kaust.edu.sa/>

Office Phone: +966 12 808-0459

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20180511/6b7736f2/attachment-0001.html>