 |
Peano
|
|
Peano
|
MGHyPE's DG realisations all support cell-dependent operators, i.e.
PDEs where the weak form over each cell differs. Probably the simplest equation of that type is
\( - \nabla \left( \epsilon(x) \nabla \right) u = f \)
with a space-dependent material parameter \( \epsilon \). The environments' vanilla configuration is prepared for such setups, even though, by default, it uses an invariant cell operator, i.e. the same operator per cell. The other way round, the code does not exploit any knowledge about invariant operators either, which means that we ignore some very low-hanging fruits regarding runtime optimisations.
(image from the SPE10 benchmark with jumping coefficients; original URL: https://www.spe.org/web/csp/datasets/set02.htm)
Each run of a DG solver generates two classes: An abstract solver class and the particular implementation. The abstract base class is regenerated every time you rerun the Python script, as it mirrors all the entries from the script such as polynomial orders or the matrices from the Python class' constructor. The particular subclass is meant to be the place where you tailor the PDE towards your needs for all features that are not available directly through the Python script.
By default, the base class befills all interesting assembly routines with the matrices from the Python script, and the only thing the user has to write is the function
The two routines that are relevant for the discussion in this chapter are however those accessing the assembly and approximate system matrices.
By default, they work with the matrices passed in through the Python script. If the argument there is None, the generated superclass routine will be abstract, and you will have to implement these two routines yourself. If you hand in some matrix, you still can overwrite and redefine these routines yourself later on.
Yet to be written.
By default, MGHyPE takes the cell operator and inverts this one per cell with a LU decomposition. LU might be very slow. If the operator in the cell is known beforehand and is the same in all cells, the whole LU business generates a lot of MFLOPS without any sense: It is way better to invert the cell operator once, to memorise this outcome and then to reuse it.
To do so, we open the header of MySolver and declare that we will overwrite the cell operator:
Before we continue, we quickly check that the autogenerated application routine in abstract base class uses this getter too:
Hence, we only need to provide a faster implementation of getInvertedApproxSystemMatrix() and should get significantly faster performance. Obviously, we can also manually tune applyInvertedApproxSystemMatrix(), but that's for the future.
The simplest way to optimise getInvertedApproxSystemMatrix() is to cache the LU decomposition's outcome. In the example below, I assume that the we know the system matrix, but that it lacks a proper scaling with a material parameter:
In the example above, we use the base class' inverted operator on the unit cube. This is for a scaled Poisson operator, so we can use this one directly as long as we work in two dimensions. In three dimensions, we have to rescale.
This is obvious one of the simplest variants one can think of, but it illustrates the fundamental idea of a very simple opimisation: We intert the operator per cell only once and then apply this one from hereon.
We can now do exactly the same "trick" for the cell operator, too. We tell MGHyPE that we will overwrite its built-in assembly
and then we implement them in our user class:
The abstract version (study the generated code in there) has an alternative implementation of applyLocalAssemblyMatrix(), where we use the helper routine mghype::applyMatrixFromHWeightedLinearCombination() which accumulates the outcome in-situ. If we can precompute the operator, we can directly call the mat-vec and rely on the compiler to make this thing fast, hoping that the compiler recognises that the matrix is static. Alternatively, we can also cache the outcome in applyLocalAssemblyMatrix() via an additional static variable. In this case, we obtain
If you have a more complicated operator, where you can pre-assemble parts of the matrices, you effectively end up with a hybrid of this method and the vanilla one. All of our tweaks here rely on the fact that the abstract base class has the operator matrices available, so we can build our code using them. If you pass in None into the Python script, the base class routines will be abstract, i.e. you cannot use them in your bespoke code and you have to write all operators yourself.
1.9.8