dune-community / dune-gdt

Passive mirror of https://zivgitlab.uni-muenster.de/ag-ohlberger/dune-community/dune-gdt

# This file is part of the dune-gdt project:
#   https://zivgitlab.uni-muenster.de/ag-ohlberger/dune-community/dune-gdt
# Copyright 2010-2021 dune-gdt developers and contributors. All rights reserved.
# License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
#      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
#          with "runtime exception" (http://www.dune-project.org/license.html)
# Authors:
#   Felix Schindler (2010, 2013 - 2014, 2016 - 2017)
#   Rene Milk       (2017 - 2018)
#   Tobias Leibner  (2021)

dune-gdt is a DUNE module which provides a generic discretization toolbox for grid-based numerical methods. It contains building blocks - like local operators, local evaluations, local assemblers - for discretization methods and suitable discrete function spaces.

New users may best try out this module by using the git supermodule dune-gdt-super. Experienced DUNE users may go ahead. As usual, you will have to use dunecontrol, working examples are located in 'dune/gdt/test/'...

If you want to start hacking go ahead and fork us on github or gitlab!

dune-gdt's Issues

remove disable/reenable warnings include guards

see wwu-numerik/dune-stuff#20

[test.operators_projections_l2] add test for FV

[projections.lagrange] merge with oswald interpolation operator

The way the Lagrange projection is currently implemented may be problematic. The same holds for the Lagrange prolongation, but that will not be an issue after #54.
The problem arises when applying the Lagrange projection to a localizable function that is discontinuous across the faces of the grid: we carry out the projection locally by assigning the values of the local DoF vector on each entity. Since there are Lagrange points that are shared by more than entity, the value of the projected function is not unique in that point and only depends on the ordering of the elements in the grid.

I propose to merge the oswald interpolation operator into the Lagrange projection. That operator basically takes the average value at each Lagrange point, which is probably what we want (and at least deterministic). In addition we might want to statically drop to the old behaviour if the source is contained in a continuous space for performance reasons, since the current approach is correct in those cases.

[test.operators_products_swipdgpenalty] add dimensions 2 and 3

[test.operators_projections_l2] add test for RT

tmp storage provider

This regards the question of providing common tmp storage for operators/functionals/localevaluations for use on each entity. We currently do so by means of the TmpMatrix and TmpVector provider (or something like that) and by having (/polluting) the evaluate and apply signatures of everyone with tmp_matrices. Since we never benchmarked this my personal take on this was to drop any of this at some point since it became quite annoying. I talked to Steffen, however, and they are convinced that it is really bad to have dynamic memory allocation on each entity. One aspect I never considered was the shared memory parallel case where this could lead to a potential bottleneck if allocation is global and each thread hast to wait on each entity. There are four possible ways to go:

• do not have any common tmp storage at all (this might be problematic, see above), drop the tmp stuff from the signatures.
• keep common storage as a mutable member in each object. This might not be such a bad idea, if we ensure that each object is duplicated for each thread.
• have a cleaver common tmp storage provider (further detailed below).
• do nothing and leave everything as it is. I do not like the current situation as it is cumbersome to work with and does by far not cover all cases.

For the third option we would need some kind of TmpStorageProvider that can be easily created in the assembler before the grid walk and then given to anyone needing it on each entity. Therefore, local evaluations, local operators and local functionals would have to statically define what kind of storage they need. This would probably be of the form std::tuple< FieldVector< d >, FieldMatrix< d, r >, ... > for each of the three. The TmpStorageProvider would have to take all these tuples and flatten and unique them (if such a thing is possible). It would then need be able to ask each object how many objects of each kind it requires, e.g. some templatized method that defaults to 0 if we do not need anything (this is kind of done in numTmpObjectsRequired) and returns a number for each element in the tuple. In the end we could create one such TmpStorageProvider object in the assembler which would (in the best case) be given to any functor on each entity. Thus we would only have one additional argument polluting the evalaute/apply signatures but gain common tmp storage for everything we need in the evaluations.

This might be quite a lot to implement so I am also seriously considering options one or two.

[test] move headers to dune/gdt/test

Since they might be of use in dependent modules.

[operators.elliptic-swipdg] update `make_elliptic_swipdg`

[local operators and functionals] update interfaces

This is a reminder to myself to refactor the interfaces for local operators and functionals (or just drop them), after the pending refactor of the (global) operators and functionals. The current interfaces only allow for assembling into a matrix and not for matrix-free operators. In addition, they are not used anywhere atm...

add typetraits.hh for Dune::AlwaysFalse where needed

See 774af6c#commitcomment-7369692

[operators.projections] derive from base

And thus get a localizable and assemblable implementation

[prolongations] simplify implementation

Currently we duplicate quite a bit of code, regarding prolongation and projection operators. The key difference between a projection and a prolongation operator is what we expect of the source:

• A projection expects source to be localizable with respect to the given grid_view. In fact, it is more or less expected that source is localizable with respect to any grid_view.
• A prolongation knows that source is a discrete function and thus only localizable with respect to a specific grid_view (namely the one of the space the source belongs to). We expect the sources grid_view to be either a coarser version of ranges grid_view or the ranges grid_view to span a smaller (physical) domain than the sources grid_view.

Prolongations thus have to invest some effort on making the source localizable with respect to the ranges grid_view. Once this is done, projections and prolongations to the same work to project the source onto the range (which is where code is duplicated).

It would be far more elegant to let the prolongations produce a wrapped source that was localizable with respect to the ranges grid_view (using the techniques we already have implemented) and then forward that wrapped source and the range to the standard projections.

[test.operators_prolongations] add

Take the old tests and put everything in the new framework.

add locking to assembler

I just talked to Steffen and he argued that Joe implemented a per-entity locking mechanism in the assembler and they could not observe performance loss. If we ever run into problems we should investigate this (but probably directly on the level of the grid walker in dune-stuff).

clean up evaluate() and order() in local evaluations

All local evaluations have to implement one order and one evaluate method. All other should be made private and perhaps renamed, e.g. redirect_to_order or the like.

[test.stuff-la-solver] refactor

This test

• does not use the new testing framework yet
• takes too long
• could perhaps be (partly) moved to dune-stuff

[test.elliptic] refactor

The elliptic testcases, the discretizations and the expected resulsts are all properly implemented in dune-hdd, using higher level constructs from dune-hdd. We should move as much as possible to dune-gdt in order to minimize code duplication and especially double recordings of expected results.

[spaces] add Lagrange interface

And derive current CG and DG interface...

[spaces] implement copy

There are several reasons why copying spaces is not trivial, i.e 6b3ff6d. Either the backend_ spaces are not copy constructable (as in dune-fem) or one needs additional objects which have to be copied along (as with the fe_map_ for dune-pdelab). In addition we do not want to have several spaces to hold the same shared_ptr to the backend_ (this makes shared memory paralellization a pain or even impossible). Thus we should implement all copy ctors in the following way:

• copy the grid_layer of the imcoming space
• recreate the backend_ and all needed and dependent objects just as in the main ctor

This will lead to spaces that are nor really a copy of each other. On the other hand, since the grid parts and views are const in dune-gdt, if a recreated copy of one of those backend_ spaces on the same grid part or view does not behave exactly as the original one there is doom and destruction pending for the corresponding discretization module.

make use of pdelabs EntityIndexCache in our mapper

I just talked to Steffen and he argued that our hacks to obtain the correct mapping from the pdelab clasess could be streamlined by using the EntityIndexCache from dune/pdelab/gridfunctionspace/entityindexcache.hh. It is also more general and faster.

rename OS2014 related stuff

Since the publication is now 2015 either rename it to OS2015 or to another more meaningfull name. And reference the publication...

investigate bad performance of ESV07 test case

See http://users.dune-project.org/issues/273:

After letting the EOC run start on a finer grid level, EOCs are down by about 50%. SIPDG and SWIPDG do not show that problem.

cache evaluations of basefunctions

I just talked to Steffen and he argued that the evaluation of shapefunctions for higher order DG is quite expensive. I think this is also what we see in our tests for higher order DG. He suggested to implemented a caching of the evaluations using a hashed map (FieldVectors are hashable by now, this would happen in our basefunctionset). We should create a simple test case to measure possible performance impact for lower order spaces and possible performance gains for higher order spaces.

[spaces] get rid of grid_view() and grid_part()

Since we dropped shared_ptr and the spaces now hold the grid view or grid part in question per value it is expected that this

const auto grid_vw = space.grid_view()

is valid and creates a copy of the grid view held by the space. Moreover it is possible to do something like this:

OtherSpaceType other(space.grid_view());

The latter of course may fail spectactularly if space is a dune-fem space and OtherType is any grid view based type. In this case the grid_view() returned by space is just a fem wrapper of the grid_part() returned by space and not guaranteed to be valid, if space (and the grid part it holds) is destroid.

Since we will not drop dune-fem support in the near or far future (not open for discussion in this issue) I propose to replace grid_view() in SpaceInterface by grid_layer(), the type of which is either a grid part or a grid view (as is already automatically done with the whole grid provider use).

• With C++-11 most classes do not really care if they get a view or a part, those needing static information can use dune/stuff/grid/entity.hh and dune/stuff/grid/intersection.hh to generically obtain those types. This is for instance true for the assembler and grid walker.
• The first example would leave the user with a usable grid layer (which is valid even if the space is destroid):

const auto grid_layer = space.grid_layer()

• The second example

OtherSpaceType other(space.grid_layer());

would work if space and OtherSpaceType either use both grid view or both use grid parts but would not compile if i.e. OtherSpaceType requires a grid view and space yields a grid part, which would be a huge benefit. The user can then decide to call space.grid_layer().gridView() manually if he is sure that everything is right.

[projections.l2] implement using local operator and functional

We currently still do the integration manually, which is error prone and a duplication of effort. We should just use the local product operator and functional instead.

move to dune-community

This will change project and repository URLs.

[spaces] update interface and introduce intermediate interfaces

The current situation with e.g. ContinuousLagrangeBase is undesirable. It has virtual methods that are always instatiated and only work for polynomial order 1. Thus (and for other reasons, s.b.) we need another CRTP interface for CG, DG, FV. RT spaces.

• For everything except CG this will be an empty interface for the moment.
• The CG intarface should provide the current implementation of lagrange_points and dirichlet_DoFs which can be pulled in by any base.

Overall this will (hopefully) allow to match against a CG/DG interface (regardless of the backend) to get rid of all specializations, e.g. in operators/projections.hh and the like.

In addition I propose that the SpaceInterface (and the new intermediate interfaces) should get a signature like

template< class Traits, int dimDomain, int dimRange, int dimRangeCols = 1 >
class SpaceInterface;

instead of the current

template< class Traits >
class SpaceInterface;

This has the drawback that current code like

template< class S >
void do_something_with_a_space(const SpaceInterface< S >& space);

become a bit longer, e.g.

template< class S, int d, int r, int rC >
void do_something_with_a_space(const SpaceInterface< S, d, r, rC >& space);

On the other hand there are many dimension specific things in dune-gdt which currently have to check for the space dimensions manually inside a method and use a static_assert there (or a runtime switch). The added template parameters would allow for easy template matching, like

class Projection
{
public:
  template< class SourceType, class T, int d >
  void apply(const SourceType& source, const SpaceInterface< T, d, 1 >& scalar_space)
  {
    // do something
  }
}; // class Projection

do not duplicate pdelabs GridFunctionSpace

I just talked to Steffen and he argued that copying the GridFunctionSpace for every thread might get very costly. Instead we should consider to just duplicate the local function space (since the entity binding is done here) instead. I am not sure if we really suffer from this, so this is mainly a reminder to investigate.

[playground]

The following stuff works:

• basefunctionset/finitevolume.hh: rename to basefunctionset/finitevolume/default.hh
• functionals/*
• localevaluation/*
• mapper/*
• operators/*
• products/*
• spaces/block.hh
• spaces/continuouslagrange/fem.hh
• spaces/finitevolume/*
• spaces/raviartthomas/pdelab.hh

The following stuff should be removed:

• spaces/raviartthomas/fem-localfunction.hh: did never do what it was supposed to do

The following stuff should remain:

• spaces/continuouslagrange/fem-localfunctions.hh: probably does not work atm, but since @renemilk updated dune-fem-localfunctions we should give it a try again soon
• spaces/discontinuouslagrange/fem-localfunctions.hh: s.a
• spaces/discontinuousgalerkin/pdelab.hh: this does not really seem to work but we would like to have it there to be fixed soon

[products.base] -> [operators.base]

The functionality of the base products should be usable for the operators as well.

[operators.prolongation] derive from base

And thus get a localizable and assemblable implementation

[spaces.cg] enable dune-fem based for dimRange > 1

I think I had this running at some point but a quick test gave negative results (size() or mapper problem).

[spaces] drop shared_ptr wherever possible

All members of a space should be const values, whereever possible. Copy and move of a space should indeed make a copy of the members.

• In particular the backend (only exception probably the dune-fem backend space) and all dune-gdt mapper wrappers.
• grid views should also be const members, if grid parts support copy and move

[spaces.dg.pdelab] fails for higher polynomial orders

[test.operators_oswaldinterpolation] needs finish

[test.spaces_dg_ord3] takes too long

Investigate!

update tests

Use expectation macros whereever possible.

[space.constraints]

• introduce interface (this will reduce confusion in the grid walker)

[spaces.dg.pdelab] do not seem to work

Running test/operators_projections_l2_dg_pdelab.cc yields lots of errors for me, indicating together with #10 that the wrapper for DG spaces is totally broken for dune-pdelab.

[linearelliptic.discretizers.cg] add detection and handling of pure Neumann cases

We already have the DirichletDetector functor in dune-stuff. We only need to apply it and pass the information to the discretization container (so it can fix a DoF and subtract the mean afterwards).

[discretefunction.local]

• derive ConstLocalDoFVector from Stuff::LA::VectorInterface?
• drop ConstLocalDoFVector, let ConstLocalDiscreteFunction hold a vector from Stuff::LA::VectorInterface and map/copy in the ctor?

[operators.elliptic-cg] refactor

• base operators (after #33)
• add make_elliptic_cg

improve threaded assembly

While refactoring I stumbled across several places where we seem to treat the threaded assembly case differently. I think in the long run we need a clear concept of what is good and what is bad.

If I remember correctly, we do not want to dereference PerThreadValue on each grid entity (since it then has to determine which element to pick according to the current thread, and that may be slow). That was one of the reasons why we duplicate the spaces in the SystemAssembler (before we had a PerThreadValue vector of space backends in each space and dereference that on each call associated with an entity). Nevertheless, we have several wrappers in dune/gdt/assembler/wrapper.hh (the purpose of which is to wrap, e.g., a local assembler into a usable functor) that are constructed with the PerThreadValue of spaces and locally dereference that on each entity. That, however, is just as bad as having all the magic inside the spaces (if I am not mistaken about the current situation).

If that is true, we would have to duplicate the wrappers and functors within the add methods of the SystemAssembler. Most wrappers deal with local operators (and thus with local evaluations) where everything is const and should be thread safe (TM). If we drop the current tmp storage, however, (as decided in #49), those objects hold mutable tmp storage and need to be duplicated as well...

All in all I think we have a nice intertwined situation and should come up with a general concept of the good, the bad and the ugly in the context of threaded assembly. In particular at what point we want to duplicate which kind of objects, etc ...

[assemblers.system] should hold a PerThreadValue list of spaces

The system assembler should hold a PerThreadValue of test and ansatz space (does a copy of the given arguments in ctor), the wrappers should then receive a const & to those. Do we want to do the same in dune-stuff, i.e. the Grid::Walker should have a PerThreadValue of grid views?

[assembler.apply_on]

According to some thoughts we had the ApplyOn::BoundaryIntersections should return true if

!intersection.bounday() && !intersection.neighbor()

in the parallel case. Investigate! :)

[products] update source/range names

SourceType and RangeType (or SourceSpaceType and RangeSpaceType) are misleading and should be renamed to First and Second or Left and Right...

dirichlet constraints not thread safe

It may happen for certain grids that we do not detect all dirichlet DoFs of an entity. This is true, for instance, for a triangular grid of the following form, where the bottom models the domain boundary

x----b----x
 \ 1 | 2 /
  \  |  /
   \ | /
  0 \|/ 3
-----a-----
 boundary

The dirichlet DoF associated with the vertex a is detected correctly on entities 0 and 3, but not on entities 1 and 2. This is due to the fact that we identify boundaries by looking at intersections and neither entity 1 nor 2 have intersection that are on the boundary.

The result is that there are certain combinations of global DoF indices, where there exist nontrivial entries in the system matrix that are not deleted during application of the dirichlet DoFs. This is for instance the case with the global index corresponding to vertices a and b (since b is neither present in entity 0 nor in 3).

The remedy within the current approach is to check the neighbors of an entity (or rather their intersections) too, when we are trying to determine if the vertex of an entity lies on the boundary. Since we have no way to determine the maximum number of entities associated with a vertex this is a nasty task and can only be satisfactory solved for some grids, on others we would need some kind of heuristics.

I am currently implementing such an approach that should be sufficient for the near future but we might have to look into other approaches of finding the DoFs associated with the boundary. We would not have this problem if we would create a global container of dirichlet DoFs, for instance, but we might have others.

And yes, @renemilk, this did show up in those tests that have been broken forever.

add CG tests for higher orders

Implement higher orders for dune-fem based spaces

[test.elliptic_estimators] investigate new expected results

The expected results recorded by now and the actual results are off by a slow margin. This could be because of the change of the testing framework (different comparison of the results) but we need to check and make sure that this is the cause and that those are only rounding errors and nothing worse...

[discretefunction.default]

• let local_discrete_function() return a unique_ptr