Refactor Operators #198

ftalbrecht · 2020-12-08T10:19:52Z

Remaining todos:

update functionals accordingly
make all c++ tests pass
update Python bindings

This completely rewrites the (global) operators to address several issues we've had over the time. This is the main change I want to have in dune-gdt before we write the paper. Once this is in, we can have a first release of a wheel.

Notes:

the Python bindings are not adapted yet and will not compile, I am first aiming for passing headercheck and C++ tests
the current documentation is not up to date, but I do not really don't know where to put it
I was playing around with tutorials in Jupyter Notebook using xeus-cling for interactive C++ code and surrounding Markdown cells, which would be my favorite solutions (however, I get lots of segfaults as our code seems to be a bit too much still); doxygen feels clumsy and not-read-by-many to me. I will just start to explain some of the concepts here if it would help the rest of you and we'll see where this ends up, though.
I started to adapt @tobiasleibner advection-with-reconstruction operators, but did not try to compile them, all other operators should be working by now

The main changes are:

there is a cascade of three interfaces:
- BilinearFormInterface, allowing for
  - apply2(v, u) and
  - norm(u), where u and v are grid functions;
- ForwardOperatorInterface : public BilinearFormInterface, additionally allowing for
  - apply(u, range_vector), reading a grid function u and writing into a range_vector (plus convenience methods for DiscreteFunction and the like); and
- OperatorInterface : public ForwardOperatorInterface, additionally allowing for
  - apply(source_vector, range_vector), reading a source_vector and writing into a range_vector,
    which then allows for
  - jacobian(...)
    and
  - apply_inverse(...)
    all based on vectors, plus the usual convenience methods.

(The main difference between ForwardOperatorInterface and OperatorInterface is, that the former acts on grid functions only and does not require a discrete source space, while the latter requires discrete source space but may additionally implement apply for grid functions.)

For each of the three interfaces there is exactly one default class

BilinearForm, allowing to += local element/intersection bilinear forms
ForwardOperator, allowing to += local element and intersection operators
Operator, allowing to += local element and intersection operators

which then implements the rest of the respective interface (including jacobian and apply_inverse using dampened newton).
Each of the three are appendable to the GridWalker, given all additional required information.

Each more general class allows to be converted to a more specific one given additional information, e.g.

the BilinearForm can be converted to a MatrixOperator given a matrix or a discrete source space and
the ForwardOperator can be converted to an Operator given a discrete source space.

Thus, there is now only one well-defined place where things happen. I.e., to assemble a bilinear form into a matrix, one cannot append the local bilinear forms to the MatrixOperator, but creates the BilinearForm first and then converts this to a MatrixOperator. In particular, one can still use the defined bilinear form as a product for arbitrary functions.

Some examples:

products can be easily built from bilinear forms as usual
all bilinear forms and forward operators can now be applied to everything a GridFunction can handle, see use of norms
an IPDG bilinear form assembled into a matrix (not the most elegant way of bilinear form to matrix conversion for efficient test code though)
arbitrary nonlinear operators can be built by appending local operators as usual

- Operator -> ForwardOperator - DiscreteOperator -> Operator

tobiasleibner · 2020-12-08T12:58:37Z

Looks sensible on a first (very short) look. Unfortunately, I won't have time to have a closer look at this in the next few weeks, but feel free to merge this already (with the hyperbolic tests disabled if necessary, I will then fix these in a follow-up PR).

renefritze

Very happy with the overall new design.

General question: Can we make not be necessary to pass two seperate objects around controlling the logger state?

renefritze · 2020-12-08T13:20:58Z

dune/gdt/algorithms/newton.hh

+      try {
+        tried_linear_solvers.push_back(linear_solver_type);
+        jacobian_solver.apply(
+            residual, update, {{"type", linear_solver_type}, {"precision", XT::Common::to_string(0.1 * precision)}});


0.1 is a magic number

Yes, do you not like magic numbers? :) As in: this should be configurable via the options?

I mean I like magic, just not in numbers maybe? 😉
If you didn't find it necessary to change this yet, maybe just make it a meaningfully named constant?

renefritze · 2020-12-08T13:22:00Z

dune/gdt/algorithms/newton.hh

+      } catch (const XT::LA::Exceptions::linear_solver_failed&) {
+      }
+    }
+    DUNE_THROW_IF(!linear_solve_succeeded,


I would not expect that a newton algorithms internally tries a number of not-configurable linear solvers.

Mhm, but this is here for a reason, Felix from a year ago could have told you! :)

For arbitrary nonlinear operators, it happens quite often that systems are barely solvable and that the choice of the linear solver might really make a difference. How should we proceed then?

We could add an option to the newton cfg about the linear solver, where something like "auto" is the current behavior. It's not that easy however, since we would also like to be able to pass the linear solver options as a sub cfg in the newton cfg. If the user selects a single linear solver this is easy.

However, a realistic scenario would be to specify my favorite direct solver and iterative solvers, say "superlu, amg", but how do I pass the respective opts as sub cfgs to the newton?

renefritze · 2020-12-08T15:47:33Z

dune/gdt/operators/bilinear-form.hh

    local_source_in_->bind(element);
    local_range_in_->bind(element);
-    for (auto& data : element_data_) {
+    for (auto& data : bilinear_form_.element_data()) {
      const auto& local_bilinear_form = *std::get<0>(data);


If you don't mind shifting the deref to the passing site, there's an alternate way to write this:

const auto&& [local_bilinear_form, filter] = data;

??? Do we now have tuple unpacking in C++? That would be amazing!

Well, yes and no. Look up structured bindings, got introduced with 17.

renefritze · 2020-12-08T15:48:56Z

dune/gdt/operators/bilinear-form.hh

 }


 } // namespace GDT
 } // namespace Dune

+#include "matrix.hh"


If this is not here accidentally a comment explaining the reason would be good for the future.

Agree! Here on purpose. There is a comment at the top where the forward is, but another comment here will do.

renefritze · 2020-12-08T16:00:14Z

dune/gdt/operators/advection-fv-entropybased.hh

-    return operator_.source_space();
+    this->assert_matching_range(range_vector);
+    VectorType u_update = range_vector.copy();
+    u_update.set_all(0.);


Could this be changed to init with 0 and correct size? Also in a few other places I think.

Yes, agreed. We also have something like zeros_like I think.

renefritze · 2020-12-08T16:24:17Z

dune/gdt/operators/interfaces.hh

+    DUNE_THROW_IF(!opts.has_key("type"), Exceptions::operator_error, "opts = " << opts);
+    const auto type = opts.get<std::string>("type");
+    if (type == "zero") {
+      LOG_(info) << "not adding zero jacobian ..." << std::endl;


Allow adding zero operators here, but throwing an error otherwise is intentional? If there's a structural reason for this, maybe add a note?

This is a default implementation in the interface for convenience. Ops with zero jacobian only need to implement available..., not jacobian(). So yes, some docs should be there.

renefritze · 2020-12-08T16:24:27Z

dune/gdt/operators/interfaces.hh

-            1.);
-    return ret;
-  }
+    DUNE_THROW(Exceptions::operator_error, "This DiscreteOperator does not support jacobian(source_vector)!");


Unreachable code?

renefritze · 2020-12-08T16:29:06Z

dune/gdt/operators/lincomb.hh

+    this->assert_matching_range(range_vector);
+    LOG_(info) << "applying " << this->num_ops() << " operators ..." << std::endl;
+    range_vector.set_all(0);
+    auto tmp_range_vector = range_vector;


renefritze · 2020-12-08T16:31:11Z

dune/gdt/operators/lincomb.hh

  {
-    DUNE_THROW_IF(type != this->jacobian_options().at(0), Exceptions::operator_error, "type = " << type);
+    std::vector<XT::Common::Configuration> ret(1);


Could use plain Config here and return vector(cfg) instead.

renefritze · 2020-12-08T16:35:53Z

dune/gdt/operators/operator.hh

+
+  /// \}
+
+  const std::list<std::pair<std::unique_ptr<LocalElementOperatorType>, std::unique_ptr<ElementFilterType>>>&


auto return?

Do we have a view on this (when to use auto)? The old 'auto is easier' vs. 'I can see what the type is' argument. Currently, I use auto when it gets too complicated, so this could apply here. Do we want to use it more often in general?

Maybe use something like "more than one nesting level --> use auto" as a rule of thumb? If one can just use a type alias I still think it unnecessary/silly do use auto, yeah.

ftalbrecht · 2020-12-11T09:52:41Z

General question: Can we make not be necessary to pass two seperate objects around controlling the logger state?

Yes please! Lets discuss this here.

ftalbrecht · 2020-12-11T12:09:08Z

still awesome

renefritze · 2020-12-11T12:54:42Z

still awesome

I think the mail-in thingy misplaced you comment

ftalbrecht · 2020-12-11T21:16:06Z

still awesome

I think the mail-in thingy misplaced you comment

True. Was referring to the tuple-unpack...

ftalbrecht · 2021-02-07T19:05:27Z

I've just merged the required changes in dune-xt and updated the dune-gdt-super. So we can restart CI here once the new images are there and then start to tackle actual errors from this pr...

tobiasleibner · 2021-02-08T09:14:23Z

I restarted the pipeline, the first error is that dune/gdt/operators/localizable-operator.hh is still included in several places though it seems to be gone/renamed.

tobiasleibner · 2021-11-24T10:37:46Z

Moved to gitlab.

ftschindler added 30 commits November 9, 2020 10:02

[interpolations.rt] pass param to functions and bases

f13335d

[interpolations.rt] add some variants to match the others

2132902

[interpolations.default] allow to specify param

582be1c

[interpolations] fix delegation

aca0f00

[operators] introduce new interfaces

49f0331

[operators.constant] update to new interface

6e02b4a

[operators.matrix] update to new interface

8b86f18

[operators.bilinear-form] update to new interface

39122bc

[operators.operator] add new Operator

54d19ee

[operators.discrete-operator] add new DiscreteOperator

5ff22c5

update logging

99c9438

[operators] drop non-discrete jacobian() variant

de0486a

[operators] simplify logging

7708abf

[operators] rename

fb091c5

- Operator -> ForwardOperator - DiscreteOperator -> Operator

[operators] Operator -> ForwardOperator

dad8382

[operators] DiscreteOperator -> Operator

678756d

[operators.constant] rename

6331d4e

[operators.matrix] update names

5e8e982

[operators.bilinear-form] minor fixes

585c206

[operators.bilinear-form] minor updates

6f32410

[operators.constant] be a bit more silent

5edb448

[operators.identity] update to new interface

f97e8fc

[operators.matrix] be a bit more silent

10aad90

[operators.lincomb] update to new interface, add numeric operators

db2fb12

[operators.operator] add static_assert

8cf7e92

[operators.discrete-operator] whitespace

0fe7408

[operators] matrix-based.hh -> matrix.hh

5bdbb9f

[operators] operator.hh -> forward-operator.hh

7755686

[operator] discrete-operator.hh -> operator.hh

5c591ea

[interpolations.rt] fix to actually compile

bce4fb1

ftschindler added 7 commits December 7, 2020 16:46

[test...hyperbolic...] update advection op usage

997f9a7

[test.instationary...] updates

7332915

[operators.identity] fix factoy

63ffd24

[operators.lincomb] fix segfault in concatenation of multiple tmp ops

6dd6b68

[test...] minor

dcd21a0

fix docs

56b3a8e

[operators.advection] update to new API

16d1172

ftalbrecht requested review from renefritze and tobiasleibner December 8, 2020 10:19

ftalbrecht self-assigned this Dec 8, 2020

renefritze reviewed Dec 8, 2020

View reviewed changes

tobiasleibner added 7 commits August 23, 2021 12:10

[python] adapt to changes in dune-xt

67b2584

Add compile_commands.json to .gitignore

2fcb89c

Fix python errors and some warnings

49ec4c3

[ci] Update

9666d0f

Update CI again

689ef60

[github] fix autosquash action (uses main instead of master now)

a777040

Another try to fix CI

421b631

ftalbrecht mentioned this pull request Aug 26, 2021

Adapt to changes in dune-xt #201

Merged

Merge branch 'adapt_to_dune_xt_changes' into refactor-operators

992c17e

tobiasleibner mentioned this pull request Sep 28, 2021

Prepare license update #205

Merged

tobiasleibner closed this Nov 24, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Refactor Operators #198

Refactor Operators #198

ftalbrecht commented Dec 8, 2020

tobiasleibner commented Dec 8, 2020

renefritze left a comment

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020

renefritze Dec 11, 2020

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020

renefritze Dec 11, 2020

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020 •

edited

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020

renefritze Dec 8, 2020

renefritze Dec 8, 2020

renefritze Dec 8, 2020

renefritze Dec 8, 2020

ftalbrecht Dec 11, 2020

renefritze Dec 11, 2020

ftalbrecht commented Dec 11, 2020

ftalbrecht commented Dec 11, 2020 via email

renefritze commented Dec 11, 2020

ftalbrecht commented Dec 11, 2020

ftalbrecht commented Feb 7, 2021

tobiasleibner commented Feb 8, 2021

tobiasleibner commented Nov 24, 2021


		/// \}

		const std::list<std::pair<std::unique_ptr<LocalElementOperatorType>, std::unique_ptr<ElementFilterType>>>&

Refactor Operators #198

Refactor Operators #198

Conversation

ftalbrecht commented Dec 8, 2020

tobiasleibner commented Dec 8, 2020

renefritze left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ftalbrecht Dec 11, 2020 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ftalbrecht commented Dec 11, 2020

ftalbrecht commented Dec 11, 2020 via email

renefritze commented Dec 11, 2020

ftalbrecht commented Dec 11, 2020

ftalbrecht commented Feb 7, 2021

tobiasleibner commented Feb 8, 2021

tobiasleibner commented Nov 24, 2021

ftalbrecht Dec 11, 2020 •

edited