dune-pdelab issueshttps://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues2023-05-15T14:56:33Zhttps://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/196Dune Assembler integration2023-05-15T14:56:33ZSantiago Ospina De Los Ríossospinar@gmail.comDune Assembler integration## Description
During last meeting with @christi and @peter, we agreed to bring latest advances from `dune-assembler` into `dune-pdelab`.
The main idea is to unify the `dune-functions` interface of the function spaces with the requireme...## Description
During last meeting with @christi and @peter, we agreed to bring latest advances from `dune-assembler` into `dune-pdelab`.
The main idea is to unify the `dune-functions` interface of the function spaces with the requirements of PDELab and to re-implement on top of that the assemblers in a more composable form. The main reason to bring new function spaces (rather than `dune-functions` directly) is that the function spaces in `dune-assembler` are essentially a re-implementation of `PDELab::GridFunctionSpace` with extensive testing and extra features, meaning that current PDELab functionality could still be achieved when needed (e.g. local index sets for parallel grid based communications, multi-domain spaces, entity blocking, etc).
Thus, the following is a high level road map of how to achieve that:
## Road Map
### From `Dune::Assembler` into `Dune::PDELab::inline Experimental`.
The idea is to bring the objects, concepts, and helpers functions from `Dune::Assembler` into the `Dune::PDELab::inline Experimental` namespace. During the process, we can refine, document, and improve the current implementation in `dune-assembler`.
- [ ] Basic Infrastructure !603
- [ ] Google testing
- [ ] Perfetto tracing
- [ ] Threading libraries
- [ ] Common utilities
- [ ] Property Tree !606
- [ ] Entity sets definition
- [ ] Basis functions !604
- [ ] Concept definition
- [ ] New implementation
- [ ] Wrapper for `PDELab::GridFunctionSpace`
- [ ] Wrapper for `Dune::Functions`
- [ ] Assemblers & Assembler infrastructure
- [ ] Hierarchical Pattern creation (based on !561) #181
- [ ] Abstract base classes for `Forward` and `Inverse` operators
- [ ] One step operator & Mass/Stiffness assembler
- [ ] Newton operator
- [ ] PDELab local operator wrapper for new assemblers
### Integrate current PDELab practices with new infrastructure
Once, the experimental objects are in place, we can start to integrate them into PDELab current infrastructure:
- [ ] Unify boundary conditions treatment
- [ ] Use new basis concept in current assemblers
- [ ] Wrap current solvers & grid operators into `Forward` and `Inverse` operatorsSantiago Ospina De Los Ríossospinar@gmail.comSantiago Ospina De Los Ríossospinar@gmail.comhttps://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/195Support for multi-type containers2023-03-06T12:12:35ZSantiago Ospina De Los Ríossospinar@gmail.comSupport for multi-type containersWe need to discuss if is possible to implement multi-typed containers in PDELab. This comes from the discussion on https://gitlab.dune-project.org/core/dune-istl/-/issues/105#note_125741.We need to discuss if is possible to implement multi-typed containers in PDELab. This comes from the discussion on https://gitlab.dune-project.org/core/dune-istl/-/issues/105#note_125741.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/194Please consider making a release2023-06-01T13:03:49ZYuri VicPlease consider making a releaseArch linux thinks that the version is 2.8 [here](https://aur.archlinux.org/cgit/aur.git/tree/PKGBUILD?h=dune-pdelab).
Is this correct?
The last tag is from 4 years ago.Arch linux thinks that the version is 2.8 [here](https://aur.archlinux.org/cgit/aur.git/tree/PKGBUILD?h=dune-pdelab).
Is this correct?
The last tag is from 4 years ago.PDELab 2.9Santiago Ospina De Los Ríossospinar@gmail.comSantiago Ospina De Los Ríossospinar@gmail.comhttps://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/191Commit messages say that the version is 2.8 but the last tag is v2.6.0-rc12022-05-22T03:56:23ZYuri VicCommit messages say that the version is 2.8 but the last tag is v2.6.0-rc1https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/187Unnecessary setReuse warnings in NewtonMethod2021-12-09T18:32:05ZValentin WüstUnnecessary setReuse warnings in NewtonMethodNewtonMethod prints a warning every time the `setReuse` flag of the linear solver is changed (this is only relevant for AMG, as this can keep the built hierarchy through successive calls of `solve`). They way Newton's method is currenty ...NewtonMethod prints a warning every time the `setReuse` flag of the linear solver is changed (this is only relevant for AMG, as this can keep the built hierarchy through successive calls of `solve`). They way Newton's method is currenty implemented, it will calculate a new Jacobian for the first step. If `ReassembleThreshold` is set to a value larger than zero, it will then continue using this Jacobian for the next steps until the achieved reduction falls below `ReassembleThreshold`; if we have `setReuse=true` until now it doesn't need to modify that. However, once it calculates a new Jacobian it has to set `setReuse` to false to build a new hierarchy, and back again to true if the Jacobian is afterwards again kept for the next steps. So even within a single call so `apply`, `setReuse` might have to be changed multiple times, so I don't really see the point of notifying the user about it every time.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/186Parallelization of MaxNorm in NewtonMethod2021-12-09T13:16:48ZValentin WüstParallelization of MaxNorm in NewtonMethodNewtonMethod calculates the current defect, i.e. the norm of the current residual, as
```C++
if (_useMaxNorm){
auto rankMax = Backend::native(_residual).infinity_norm();
_result.defect = _gridOperator.testGridFunctionS...NewtonMethod calculates the current defect, i.e. the norm of the current residual, as
```C++
if (_useMaxNorm){
auto rankMax = Backend::native(_residual).infinity_norm();
_result.defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
}
else
_result.defect = _linearSolver.norm(_residual);
```
This causes an issue if the residual in the overlap is not consistent across processors (which does not seem to be a requirement for NewtonMethod with ISTLBackend_BCGS_AMG_ILU0). Using the regular norm works, since `_linearSolver.norm` uses `disjointDot`, which masks the overlap from the scalar product. However, `infinity_norm` includes the overlap, and therefore the max norm does not produce the correct result.
Is this the intended behaviour? Or should the residual that is produced by the grid operator anyway be consistent across processors, and this is an issue of my implementation? This could explain why e.g. ISTLBackend_OVLP_GMRES_ILU0 for me does not converge in parallel. On a related note, is the correct behaviour of residuals, Jacobians etc. in terms of consistency across processors, and the respective requirements of the different linear solvers, documented anywhere? I did not manage to find anything beyond the occasional comment in the code.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/185Ordering interface change in !555 not fulfilled by all ordering implementations2023-05-12T15:50:00ZChristian EngwerOrdering interface change in !555 not fulfilled by all ordering implementationsIn !555 a new methods `.size(prefix)` and `.size(prefix,entityindex)` were introduced on the orderings.
(I'm not sure how they exactly relate, it seems that `.size(prefix)` is implemented in `gridviewordering.hh` and that all other imple...In !555 a new methods `.size(prefix)` and `.size(prefix,entityindex)` were introduced on the orderings.
(I'm not sure how they exactly relate, it seems that `.size(prefix)` is implemented in `gridviewordering.hh` and that all other implementations require `.size(prefix,entityindex)`)
In any case ... not all ordering implementations were updated accordinly. At least the `ChunkedBlockOrdering` is not covered.
A quick `grep` suggests that
```
chunkedblockordering.hh
interleavedordering.hh
leaforderingbase.hh
singlecodimleafordering.hh
subordering.hh
```
should be checked.PDELab 2.9Santiago Ospina De Los Ríossospinar@gmail.comSantiago Ospina De Los Ríossospinar@gmail.comhttps://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/182PDELab tutorials CI are broken2023-05-12T15:49:56ZSantiago Ospina De Los Ríossospinar@gmail.comPDELab tutorials CI are brokenIt seems to be some time since the PDELab tutorials are broken in this CI: https://gitlab.dune-project.org/pdelab/dune-pdelab/-/pipelines?page=1&scope=all&ref=master
* The CI does not report an error on configuration failure of the modu...It seems to be some time since the PDELab tutorials are broken in this CI: https://gitlab.dune-project.org/pdelab/dune-pdelab/-/pipelines?page=1&scope=all&ref=master
* The CI does not report an error on configuration failure of the module.
* The error looks like this:
```
CMake Error at overview/CMakeLists.txt:1 (add_latexmk_document):
Unknown CMake command "add_latexmk_document".
```PDELab 2.9https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/181Improve ISTL backend to have sparse-sparse patterns2023-05-12T15:43:53ZSantiago Ospina De Los Ríossospinar@gmail.comImprove ISTL backend to have sparse-sparse patterns## Summary
Currently, the ISTL Backend is able to create nested `BCRSMatrix` blocks by using the backend tag `bcrs`. However, only the deepest level of the blocks is treated as sparse by the `BCRSPattern`. This is only always true. For ...## Summary
Currently, the ISTL Backend is able to create nested `BCRSMatrix` blocks by using the backend tag `bcrs`. However, only the deepest level of the blocks is treated as sparse by the `BCRSPattern`. This is only always true. For example,
```c++
using DefaultVBE = PDELab::ISTL::VectorBackend<>;
using LeafGFS = PDELab::GridFucntionSpace<...,DefaultVBE>;
// Sparse vector backend for higher nodes
using SparseVBE = PDELab::ISTL::VectorBackend<PDELab::ISTL::Blocking::bcrs>;
using PGFS1 = PDELab::PowerGridFunctionSpace<k1, LeafGFS, SparseVBE, PDELab::EntityBlockedOrderingTag>;
using PGFS2 = PDELab::PowerGridFunctionSpace<k2, PGFS, SparseVBE>;
using Domain = Dune::PDELab::Backend::Vector<GFSU,double>; // -> BlockVector<BlockVector<FieldVector<double,1>>>
using Range = Dune::PDELab::Backend::Vector<GFSV,double>; // -> BlockVector<BlockVector<FieldVector<double,1>>>
using MB = PDELab::ISTL::BCRSMatrixBackend<>;
using Jacobian = Dune::PDELab::Backend::Matrix<MB,Domain,Range,JF>; // -> BCRSMatrix<BCRSMatrix<FieldMatrix<double,1,1>>>
```
The problem with the resulting `Jacobian` is that the outer block is sparse because it represents the sparsity links between entities, while the inner one represents the sparsity of the `PGFS1` nodes. Both may be sparse!
### Proposal
Modify `BCRSPattern` to hold inner sparsity patterns, and when a new link is created, a sub-pattern should also be created. I will show this in a coming MR.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/179SingleCodimMapper is broken2022-02-11T18:22:20ZSantiago Ospina De Los Ríossospinar@gmail.comSingleCodimMapper is brokenThe `OrderingTag` for single codimension `SingleCodimMapper` (i.e. Fast Finite Volume) is entirely broken.
With a very quick check on `test-transport-ccfv.cc` by adding `SingleCodimMapper` to the grid function space definition it is eas...The `OrderingTag` for single codimension `SingleCodimMapper` (i.e. Fast Finite Volume) is entirely broken.
With a very quick check on `test-transport-ccfv.cc` by adding `SingleCodimMapper` to the grid function space definition it is easy to realize that:
* The `LFSIndexCacheBase` specialization lacks definitions for `DOFIndex` and `ContainerIndex`.
* The `SimpleContainerIndex` has no hash and cannot be used in the unordered maps for constraint containers.
* The `SimpleDOFIndexAccessor` has no methods for `geometryIndex` and `entityIndex`.
* The `SingleCodimLeafOrdering` has no `_gfs_data`.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/178Bug in ISTL container index depth calculation2021-07-27T11:55:44ZSantiago Ospina De Los Ríossospinar@gmail.comBug in ISTL container index depth calculation### Summary
The `MultiIndex<std::size_t,depth>` which serves to find vector values in a nested vector container (a.k.a. `ContainerIndex`), [is calculated](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/blob/4e0e7bbeb456e181855bf5b...### Summary
The `MultiIndex<std::size_t,depth>` which serves to find vector values in a nested vector container (a.k.a. `ContainerIndex`), [is calculated](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/blob/4e0e7bbeb456e181855bf5bd964c8db3d714d717/dune/pdelab/ordering/transformations.hh#L58) using the typetree algorithm for accumulating values. That is, the `accummulate_value<...>`, which is a meta-algorithm that statically traverses *and accumulates a value* on the tree with a depth-first traversal with post-order accesses to apply functors (follow blue dots `#2A7FFF`). See https://en.wikipedia.org/wiki/Tree_traversal#Post-order
![image](/uploads/4f33fb87e1879caff318c967d833505d/image.png)
The calculation looks like this:
```c++
static const std::size_t ci_depth =
TypeTree::AccumulateValue<RootGFS,
extract_max_container_depth,
TypeTree::max<std::size_t>,
0,
TypeTree::plus<std::size_t>
>::result + 1;
```
where the value is initialized at 0, a `max` functor is used for _sibling-to-sibling_* reduction and a `plus` functor is used for child-to-parent. Finally, the `extract_max_container_depth` returns `1` or `0` depending on whether the node uses blocking or not. Here the idea is to find the maximum depth of the blocked nodes so that the `ContainerIndex` is able to represent all possible multi-indices in the nested vector structure.
The problem is that these functors and algorithm do not accumulate the desired depth of the blocking.
To see this, let's make an example: assume that we are working with the tree shown above. Nodes `B`, `G`, and `I` return `1` meaning that these nodes are blocked, otherwise `0`. Looking at the tree, it is clear that the `CointainerIndex` depth should be `2`. However, when the tree traversal makes the _sibling-to-sibling_ reduction between nodes `B` and `H`, the algorithm accumulates `1` instead of `0`, then and the overall depth accumulates `3` instead of the expected `2`.
*Note: I said "sibling-to-sibling" because this is what the documentation says, but this is indeed "sibling-to-next-leaf".
### Proposed ways to solve this
Two possible ways:
* Add a parent-to-child reduction in `dune-typetree` (https://gitlab.dune-project.org/copasi/dune-typetree/-/commit/d5f7e62a6bed5781b12a8ba6a3f332e1aac9f342) and add a add a first value functor for that reduction in this call in pdelab (https://gitlab.dune-project.org/copasi/dune-pdelab/-/commit/72f0a14e809328f27d3a909e9908c925ebb8da19)
* Change the accumulation algorithm in `dune-typetree` that is able to perform a complete depth-first traversal at compile time (https://gitlab.dune-project.org/staging/dune-typetree/-/merge_requests/93) and change to use it here.
In either case, to my understanding, this requires a change in `dune-typetree`...
### Related issues
https://gitlab.dune-project.org/staging/dune-typetree/-/merge_requests/93https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/177Remove AMG code in non-AMG ISTL solvers2021-05-04T15:02:56ZSantiago Ospina De Los Ríossospinar@gmail.comRemove AMG code in non-AMG ISTL solversThe base class for non overlapping ISTL solvers `ISTLBackend_NOVLP_BASE_PREC` uses the method [`createIndexSetAndProjectForAMG`](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/blob/875d3f3e8b5b61eb2b9c03c76da3fd3dc20fda10/dune/pdel...The base class for non overlapping ISTL solvers `ISTLBackend_NOVLP_BASE_PREC` uses the method [`createIndexSetAndProjectForAMG`](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/blob/875d3f3e8b5b61eb2b9c03c76da3fd3dc20fda10/dune/pdelab/backend/istl/novlpistlsolverbackend.hh#L766) to create a parallel index set for AMG. This line of code [dates back](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/commit/d13f7641dc663289ca44eb3a5dbdf90b0c783695#1f8d9e16146e46e11e85c226340c358b98cc2976_0_635) to the time where the class `ISTLBackend_AMG_NOVLP` did not exist. Nowadays, (I think) support for AMG should not be added to the non-AMG solvers.
This is even necessary when the underlying code does not hold [the assumptions](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/blob/875d3f3e8b5b61eb2b9c03c76da3fd3dc20fda10/dune/pdelab/backend/istl/parallelhelper.hh#L264) that the `createIndexSetAndProjectForAMG()` function has:
```c++
* \warning This function silenty assumes that the matrix only has a single level
* of blocking and will not work correctly otherwise. Also note that AMG
* will only work correctly for P1 discretisations.
```https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/176Bug in nonlinearconvectiondiffusionfem.hh2021-02-18T12:45:30ZRené HeßBug in nonlinearconvectiondiffusionfem.hh@chaiyod.kamthorncharoen found a bug in nonlinearconvectiondiffusionfem.hh:
The documentation talks about solving
q(x,u) - D(x) v(u) \nabla w(u) = f(u)
but the code seems to be doing
q(x, w(u)) - D(x) v(w(u)) \nabla w(u) = f(w(u))
T...@chaiyod.kamthorncharoen found a bug in nonlinearconvectiondiffusionfem.hh:
The documentation talks about solving
q(x,u) - D(x) v(u) \nabla w(u) = f(u)
but the code seems to be doing
q(x, w(u)) - D(x) v(w(u)) \nabla w(u) = f(w(u))
The implementation is quite confusing as the the names for u and w are mixed up but I think this is what happens.
Now the question is: Is this a documentation bug or a bug in the localoperator? What is the equation that this local operator is supposed to solve?
PS: There are some obvious documentation bugs in the description of the Neumann boundary condition that should definitely be fixed. But there it is clear what to do.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/171Pass entity and intersection on pattern methods2023-05-12T15:52:37ZSantiago Ospina De Los Ríossospinar@gmail.comPass entity and intersection on pattern methods### Summary
Pattern methods have an interface that does not allow to check where in the grid the method is being called, that is, the get local function space objects but they don't get the entity or intersection. This constraint local ...### Summary
Pattern methods have an interface that does not allow to check where in the grid the method is being called, that is, the get local function space objects but they don't get the entity or intersection. This constraint local operator implementors to make decisions on the pattern based on its local context. Moreover, it does not allow combine objects that join several local operators to systematically branch computations based on the context of the local calls (e.g. !531).
* I encountered this for multidomain settings and had to implement [awful workarounds](https://gitlab.dune-project.org/copasi/dune-copasi/-/blob/84da61987aa8b1dc8ac3e1b7051f81d41209d590/dune/copasi/local_operator/diffusion_reaction/multidomain.hh#L268)
* @dominic has also found the same issue while working on skeleton terms in part of the domain in an C0-interior penalty method.
* @linus.seelinger applications may also get some benefit from this ([conditional pattern assembly](https://gitlab.dune-project.org/pdelab/dune-pdelab/-/blob/add0d49c446c25057493ce46ba81f70274992064/dune/pdelab/backend/istl/geneo/localoperator_ovlp_region.hh#L61))
Passing entity and intersection is not a big deal as this information is in the engines anyways and no extra performance cost is suffered if these are not used.
### Proposal
Change the pattern interface to one that passes the entity and intersection. e.g., move from this
```c++
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_volume( const LFSU& lfsu, const LFSV& lfsv, LocalPattern& pattern) const;
```
to this
```c++
template<typename EG, typename LFSU, typename LFSV, typename LocalPattern>
void pattern_volume(const EG& eg, const LFSU& lfsu, const LFSV& lfsv, LocalPattern& pattern) const;
```
### Compatibility
The idea would be to deprecate the methods with no entity/intersection information. Although every local operator will be affected by this, most people use the full pattern methods that PDELab provides, so most people wouldn't even notice the change. Those who have custom patterns would have one or two releases to move to the new interface.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/169Move Matrix from NewtonMethod to StationaryLinearProblemSolver2023-05-12T21:06:51ZRené HeßMove Matrix from NewtonMethod to StationaryLinearProblemSolverIn !467 we discussed some changes to the `StationaryLinearProblemSolver` the `NewtonMethod` and the solver backends. @christi had the idea to move the matrix from `NewtonMethod` to the `StationaryLinearProblemSolver` class, add a new met...In !467 we discussed some changes to the `StationaryLinearProblemSolver` the `NewtonMethod` and the solver backends. @christi had the idea to move the matrix from `NewtonMethod` to the `StationaryLinearProblemSolver` class, add a new method that forces reassemblation of the matrix and then use this to solve the linear system in the `NewtonMethod`.
See also the discussion in !467, option 3.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/159Recipe TODOs2020-06-05T16:00:28ZLinus SeelingerRecipe TODOsIn order to keep track of recipe development, I'll post the list we came up during a group meeting a while ago. A lot of them have been finished already, some as part of dune-grid recipes.
For some topics we still need volunteers; Anyon...In order to keep track of recipe development, I'll post the list we came up during a group meeting a while ago. A lot of them have been finished already, some as part of dune-grid recipes.
For some topics we still need volunteers; Anyone is highly welcome to contribute! (Writing a recipe comes down to a simple doxygen page and an example code it draws code snippets from, see !405)
Volunteers needed:
* [ ] Build systems of function spaces, orderings
* [ ] TypeTree for systems (short intro, retrieving children etc, link to systems)
* [ ] EntitySet
Planned:
* [ ] Calculate nonzeros per row estimate (mention blocking) - Dominic
* [ ] Calculating errors (H1, L2, link to "Integrating grid functions") - Peter
* [ ] Setting up constraints in PDELab (which classes exist (as a list), also for parts of systems) - Chaiyod
* [ ] ParameterTree (parsing from file/console, get/set) - Ole
Done:
* [x] Operator splitting - Michal
* [x] Setting up a solver
* [x] - Core - Linus
* [x] - PDELab - Linus
* [x] Grid functions - Linus
* [x] CollectiveCommunication - Michal
* [x] More detail recipes vs tutorials - Linus
* [x] Blocking - Anne
* [x] Mesh transformation - Anne
* [x] Loop over grid (refer to it in "Integrating grid functions") - Peter
* [x] Use quadrature rule (iterate etc) - Peter
Later:
Setting up grids (UG, Yasp, gmsh, factory, ... show all options)
Codegenhttps://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/157Get rid of gridexamples.hh in dune/pdelab/test2021-03-22T15:02:17ZDominic Kempfdominic.kempf@iwr.uni-heidelberg.deGet rid of gridexamples.hh in dune/pdelab/testThe header provides grid implementations fro PDELab tests, however most of them are redundant implementations of the `StructuredGridFactory`.The header provides grid implementations fro PDELab tests, however most of them are redundant implementations of the `StructuredGridFactory`.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/156Unguarded use of SuperLU backend2020-10-14T15:17:48ZDominic Kempfdominic.kempf@iwr.uni-heidelberg.deUnguarded use of SuperLU backendWhile running tests for the PDELab release on my home desktop machine, I again stumbled over some tests (here `testinstationary`) using the SuperLU solver backend without checking that it exists. This results in compilation failure inste...While running tests for the PDELab release on my home desktop machine, I again stumbled over some tests (here `testinstationary`) using the SuperLU solver backend without checking that it exists. This results in compilation failure instead of test skipping.https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/154testtimedependentboundary_ovlpqk.cc: should test be run in parallel?2020-05-27T10:30:34ZAnsgar Burchardtansgar.burchardt@tu-dresden.detesttimedependentboundary_ovlpqk.cc: should test be run in parallel?I noticed a new testtimedependentboundary_ovlpqk.cc test. It's only run sequentially, but the comment at the beginning mentions a bug in an overlapping scheme.
Shouldn't the test be run in parallel as well?I noticed a new testtimedependentboundary_ovlpqk.cc test. It's only run sequentially, but the comment at the beginning mentions a bug in an overlapping scheme.
Shouldn't the test be run in parallel as well?https://gitlab.dune-project.org/pdelab/dune-pdelab/-/issues/153Refactor LocalMatrix2023-05-12T21:23:03ZSantiago Ospina De Los Ríossospinar@gmail.comRefactor LocalMatrix### Description
Local matrices are something which users are exposed to work with when implementing the local operators. However, it seems to me that it is not clear how to work with them beyond the `accumulate` method. For instance, It...### Description
Local matrices are something which users are exposed to work with when implementing the local operators. However, it seems to me that it is not clear how to work with them beyond the `accumulate` method. For instance, It took me a long time to realize one receives a `WeightedMatrixAccumulationView` instead of a `LocalMatrix` and it took me even longer to realize that in some cases one might receive `AliasedMatrixView` instead. The local matrix interface that a user should follow to fit the different types is not documented nor the way they are used within PDELab.
Now, I want to add yet another one: `SparseLocalMatrix`. If I just add it without organizing the problems I described above, the problem will become more apparent and it will just make switching between the different types even more difficult.
### Proposal
* Rename `LocalMatrix` to `DenseLocalMatrix`. Of course, using the corresponding deprecation warnings.
* Define an interface suitable for the different implementations that a local matrix should follow and document it. This can be enforced using concepts which in case of a non-conforming interface should give a descriptive message to the user.
* Separate the weighted view from the `LocalMatrix` implementation. This is an unnecessary redundancy on the current implementation. Furthermore, doing so may allow using the weighted view regardless of the actual local matrix implementation.
* Extend the `AliasedMatrixView` to also work with non-entity-blocked global matrices. This needs discussion with the FastDG guys (@dominic and @rhess) to see whether this is possible at all so that performance is not affected :red_car: :dash:.
* Leave the user to choose the type of local matrix to use. This ideally would happen at `GridOperator` level but we already have quite a bulky signature so I hear offers for other options :ear:.
### How to test the implementation?
* All current tests should keep passing after the refactoring.