Merge branch 'feature/navier-stokes' into 'master'

Add Navier-Stokes example

See merge request dominic/dune-perftool!197

cmake/modules/GeneratedSystemtests.cmake

@@ -67,10 +67,15 @@ function(dune_add_formcompiler_system_test)
 
       _add_test(NAME ${tname}
                 COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env ${SYSTEMTEST_SCRIPT}
-                      --exec ${tname}
-                      --ini "${CMAKE_CURRENT_BINARY_DIR}/${inifile}"
-                      --source ${CMAKE_CURRENT_SOURCE_DIR}
-                )
+                --exec ${tname}
+                --ini "${CMAKE_CURRENT_BINARY_DIR}/${inifile}"
+                --source ${CMAKE_CURRENT_SOURCE_DIR}
+                --mpi-exec "${MPIEXEC}"
+                --mpi-numprocflag=${MPIEXEC_NUMPROC_FLAG}
+                --mpi-preflags "${MPIEXEC_PREFLAGS}"
+                --mpi-postflags "${MPIEXEC_POSTFLAGS}"
+                --max-processors=${DUNE_MAX_TEST_CORES}
+               )
 
       set_tests_properties(${tname} PROPERTIES SKIP_RETURN_CODE 77)
       set_tests_properties(${tname} PROPERTIES TIMEOUT 60)

python/dune/perftool/compile.py

@@ -16,7 +16,9 @@ from ufl.algorithms.formfiles import interpret_ufl_namespace
 from dune.perftool.generation import (delete_cache_items,
                                       global_context,
                                       )
-from dune.perftool.options import get_option, initialize_options
+from dune.perftool.options import (get_option,
+                                   initialize_options,
+                                   )
 from dune.perftool.pdelab.driver import generate_driver
 from dune.perftool.pdelab.localoperator import (generate_localoperator_basefile,
                                                 generate_localoperator_file,

python/dune/perftool/options.py

@@ -42,6 +42,8 @@ class PerftoolOptionsArray(ImmutableRecord):
     explicit_time_stepping = PerftoolOption(default=False, helpstr="use explicit time stepping")
     exact_solution_expression = PerftoolOption(helpstr="name of the exact solution expression in the ufl file")
     compare_l2errorsquared = PerftoolOption(helpstr="maximal allowed l2 error squared of difference between numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)")
+    l2error_tree_path = PerftoolOption(default=None, helpstr="Tree pathes that should be considered for l2 error calculation. Default None means we take all of them into account.")
+    interactive = PerftoolOption(default=False, helpstr="whether the optimization process should be guided interactively (also useful for debugging)")
     print_transformations = PerftoolOption(default=False, helpstr="print out dot files after ufl tree transformations")
     print_transformations_dir = PerftoolOption(default=".", helpstr="place where to put dot files (can be omitted)")
     quadrature_order = PerftoolOption(_type=int, helpstr="Quadrature order used for all integrals.")
@@ -73,7 +75,8 @@ class PerftoolOptionsArray(ImmutableRecord):
     precompute_quadrature_info = PerftoolOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
     blockstructured = PerftoolOption(default=False, helpstr="Use block structure")
     number_of_blocks = PerftoolOption(default=1, helpstr="Number of sub blocks in one direction")
-
+    overlapping = PerftoolOption(default=False, helpstr="Use an overlapping solver and constraints. You still need to make sure to construct a grid with overlap! The parallel option will be set automatically.")
+    parallel = PerftoolOption(default=False, helpstr="Mark that this program should be run in parallel. If set to true the c++ code will check that there are more than 1 MPI-ranks involved and the error computation will use communication.")
 
 # Until more sophisticated logic is needed, we keep the actual option data in this module
 _options = PerftoolOptionsArray()
@@ -124,6 +127,9 @@ def process_options(opt):
     if opt.sumfact:
         opt = opt.copy(unroll_dimension_loops=True)
 
+    if opt.overlapping:
+        opt = opt.copy(parallel=True)
+
     return opt
 
 

python/dune/perftool/pdelab/driver/__init__.py

@@ -251,10 +251,24 @@ def name_mpihelper():
     return name
 
 
+@preamble
+def check_parallel_execution():
+    from dune.perftool.pdelab.driver.gridfunctionspace import name_leafview
+    gv = name_leafview()
+    return ["if ({}.comm().size()==1){{".format(gv),
+            '  std::cout << "This program should be run in parallel!"  << std::endl;',
+            "  return 1;",
+            "}"]
+
+
 def generate_driver(formdatas, data):
     # The driver module uses a global dictionary for storing necessary data
     set_driver_data(formdatas, data)
 
+    # Add check to c++ file if this program should only be used in parallel mode
+    if get_option("parallel"):
+        check_parallel_execution()
+
     # Entrypoint for driver generation
     if get_option("opcounter") or get_option("time_opcounter"):
         if get_option("time_opcounter"):

python/dune/perftool/pdelab/driver/constraints.py

@@ -24,14 +24,26 @@ def name_assembled_constraints():
     return name
 
 
+def has_dirichlet_constraints(is_dirichlet):
+    if isinstance(is_dirichlet, (list, tuple)):
+        return any(bool(d) for d in is_dirichlet)
+    else:
+        return bool(is_dirichlet)
+
+
 @preamble
 def assemble_constraints(name):
     element = get_trial_element()
     gfs = name_trial_gfs()
     is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
-    bctype_function = name_bctype_function(element, is_dirichlet)
-    return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
-                                                           gfs,
+    if has_dirichlet_constraints(is_dirichlet):
+        bctype_function = name_bctype_function(element, is_dirichlet)
+        return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
+                                                               gfs,
+                                                               name,
+                                                               )
+    else:
+        return "Dune::PDELab::constraints({}, {});".format(gfs,
                                                            name,
                                                            )
 

python/dune/perftool/pdelab/driver/error.py

@@ -9,7 +9,8 @@ from dune.perftool.pdelab.driver import (get_formdata,
                                          get_trial_element,
                                          preprocess_leaf_data,
                                          )
-from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
+                                                           name_trial_gfs,
                                                            name_trial_subgfs,
                                                            type_range,
                                                            )
@@ -104,14 +105,21 @@ def _accumulate_L2_squared(treepath):
 
     strtp = ", ".join(str(t) for t in treepath)
 
+    gv = name_leafview()
+    sum_error_over_ranks = ""
+    if get_option("parallel"):
+        sum_error_over_ranks = "  err = {}.comm().sum(err);".format(gv)
     return ["{",
             "  // L2 error squared of difference between numerical",
             "  // solution and the interpolation of exact solution",
             "  // for treepath ({})".format(strtp),
             "  typename decltype({})::Traits::RangeType err(0.0);".format(dsa),
             "  Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
+            sum_error_over_ranks,
             "  {} += err;".format(accum_error),
-            "  std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
+            "  if ({}.comm().rank() == 0){{".format(gv),
+            "    std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
+            "  }"
             "}",
             ]
 
@@ -139,8 +147,13 @@ def treepath_to_index(element, treepath, offset=0):
 def accumulate_L2_squared():
     element = get_trial_element()
     if isinstance(element, MixedElement):
-        for i in range(element.value_size()):
-            _accumulate_L2_squared(get_treepath(element, i))
+        tree_pathes = (True,) * element.value_size()
+        if get_option("l2error_tree_path") is not None:
+            tree_pathes = list(map(int, get_option("l2error_tree_path").split(',')))
+            assert len(tree_pathes) == element.value_size()
+        for i, path in enumerate(tree_pathes):
+            if path:
+                _accumulate_L2_squared(get_treepath(element, i))
     else:
         _accumulate_L2_squared(())
 
@@ -160,12 +173,15 @@ def name_accumulated_L2_error():
 @preamble
 def compare_L2_squared():
     accumulate_L2_squared()
+    gv = name_leafview()
 
     accum_error = name_accumulated_L2_error()
     fail = name_test_fail_variable()
     return ["using std::abs;",
             "using std::isnan;",
-            "std::cout << \"\\nl2errorsquared: \" << {} << std::endl << std::endl;".format(accum_error),
+            "if ({}.comm().rank() == 0){{".format(gv),
+            "  std::cout << \"\\nl2errorsquared: \" << {} << std::endl << std::endl;".format(accum_error),
+            "}",
             "if (isnan({0}) or abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
             "  {} = true;".format(fail)]
 

python/dune/perftool/pdelab/driver/gridfunctionspace.py

@@ -367,6 +367,18 @@ def type_orderingtag(leaf):
         return "Dune::PDELab::EntityBlockedOrderingTag"
 
 
+@preamble
+def typedef_overlapping_dirichlet_constraintsassembler(name):
+    include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
+    return "using {} = Dune::PDELab::ConformingDirichletConstraints;".format(name)
+
+
+@preamble
+def typedef_p0parallel_constraintsassembler(name):
+    include_file("dune/pdelab/constraints/p0.hh", filetag="driver")
+    return "using {} = Dune::PDELab::P0ParallelConstraints;".format(name)
+
+
 @preamble
 def typedef_dirichlet_constraintsassembler(name):
     include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
@@ -380,14 +392,21 @@ def typedef_no_constraintsassembler(name):
 
 def type_constraintsassembler(is_dirichlet):
     assert isinstance(is_dirichlet, bool)
-    if is_dirichlet:
+    overlapping = get_option("overlapping")
+    if is_dirichlet and not overlapping:
         name = "DirichletConstraintsAssember"
         typedef_dirichlet_constraintsassembler(name)
-        return name
+    elif is_dirichlet and overlapping:
+        name = "OverlappingConformingDirichletConstraints"
+        typedef_overlapping_dirichlet_constraintsassembler(name)
+    elif not is_dirichlet and overlapping:
+        name = "P0ParallelConstraints"
+        typedef_p0parallel_constraintsassembler(name)
     else:
+        assert not is_dirichlet and not overlapping
         name = "NoConstraintsAssembler"
         typedef_no_constraintsassembler(name)
-        return name
+    return name
 
 
 def name_trial_subgfs(treepath):

python/dune/perftool/pdelab/driver/gridoperator.py

@@ -170,7 +170,8 @@ def define_parameters(name, formdata):
     return "{} {};".format(partype, name)
 
 
-def name_parameters(formdata):
+def name_parameters(formdata, define=True):
     name = form_name_suffix("params", formdata)
-    define_parameters(name, formdata)
+    if define:
+        define_parameters(name, formdata)
     return name

python/dune/perftool/pdelab/driver/instationary.py

@@ -14,7 +14,8 @@ from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
 from dune.perftool.pdelab.driver.gridoperator import (name_gridoperator,
                                                       name_parameters,
                                                       type_gridoperator,)
-from dune.perftool.pdelab.driver.constraints import (name_bctype_function,
+from dune.perftool.pdelab.driver.constraints import (has_dirichlet_constraints,
+                                                     name_bctype_function,
                                                      name_constraintscontainer,
                                                      )
 from dune.perftool.pdelab.driver.interpolate import (interpolate_dirichlet_data,
@@ -53,24 +54,35 @@ def time_loop():
     params = name_parameters(formdata)
     time = name_time()
     element = get_trial_element()
-    is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
-    bctype = name_bctype_function(element, is_dirichlet)
-    gfs = name_trial_gfs()
-    cc = name_constraintscontainer()
     vector_type = type_vector(formdata)
     vector = name_vector(formdata)
     interpolate_dirichlet_data(vector)
 
+    is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+    assemble_new_constraints = ""
+    if has_dirichlet_constraints(is_dirichlet):
+        bctype = name_bctype_function(element, is_dirichlet)
+        gfs = name_trial_gfs()
+        cc = name_constraintscontainer()
+        assemble_new_constraints = ("  // Assemble constraints for new time step\n"
+                                    "  {}.setTime({}+dt);\n"
+                                    "  Dune::PDELab::constraints({}, {}, {});\n"
+                                    "\n".format(params, time, bctype, gfs, cc)
+                                    )
+
     # Choose between explicit and implicit time stepping
     explicit = get_option('explicit_time_stepping')
     if explicit:
         osm = name_explicitonestepmethod()
         apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
     else:
+        osm = name_onestepmethod()
+    if has_dirichlet_constraints(is_dirichlet):
         dirichlet = preprocess_leaf_data(element, "interpolate_expression")
         boundary = name_boundary_function(element, dirichlet)
-        osm = name_onestepmethod()
         apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
+    else:
+        apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
 
     # Setup visualization
     visualize_initial_condition()
@@ -82,10 +94,7 @@ def time_loop():
             "int step_number(0);"
             "int output_every_nth = {}.get<int>(\"instat.output_every_nth\", 1);".format(ini),
             "while (time<T-1e-8){",
-            "  // Assemble constraints for new time step",
-            "  {}.setTime({}+dt);".format(params, time),
-            "  Dune::PDELab::constraints({}, {}, {});".format(bctype, gfs, cc),
-            "",
+            "{}".format(assemble_new_constraints),
             "  // Do time step",
             "  {} {}new({});".format(vector_type, vector, vector),
             "  {}".format(apply_call),

python/dune/perftool/pdelab/driver/solve.py

@@ -7,7 +7,12 @@ from dune.perftool.pdelab.driver import (form_name_suffix,
                                          is_linear,
                                          name_initree,
                                          )
-from dune.perftool.pdelab.driver.gridfunctionspace import name_trial_gfs
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+                                                           type_trial_gfs,
+                                                           )
+from dune.perftool.pdelab.driver.constraints import (type_constraintscontainer,
+                                                     name_assembled_constraints,
+                                                     )
 from dune.perftool.pdelab.driver.gridoperator import (name_gridoperator,
                                                       type_gridoperator,
                                                       )
@@ -104,7 +109,12 @@ def define_vector(name, formdata):
 @preamble
 def typedef_linearsolver(name):
     include_file("dune/pdelab/backend/istl.hh", filetag="driver")
-    return "using {} = Dune::PDELab::ISTLBackend_SEQ_SuperLU;".format(name)
+    if get_option('overlapping'):
+        gfs = type_trial_gfs()
+        cc = type_constraintscontainer()
+        return "using {} = Dune::PDELab::ISTLBackend_OVLP_BCGS_ILU0<{},{}>;".format(name, gfs, cc)
+    else:
+        return "using {} = Dune::PDELab::ISTLBackend_SEQ_SuperLU;".format(name)
 
 
 def type_linearsolver():
@@ -116,7 +126,12 @@ def type_linearsolver():
 @preamble
 def define_linearsolver(name):
     lstype = type_linearsolver()
-    return "{} {}(false);".format(lstype, name)
+    if get_option('overlapping'):
+        gfs = name_trial_gfs()
+        cc = name_assembled_constraints()
+        return "{} {}({}, {});".format(lstype, name, gfs, cc)
+    else:
+        return "{} {}(false);".format(lstype, name)
 
 
 def name_linearsolver():

python/dune/perftool/pdelab/parameter.py

@@ -42,7 +42,7 @@ def define_parameterclass(name):
 def name_paramclass():
     formdata = get_global_context_value("formdata")
     from dune.perftool.pdelab.driver.gridoperator import name_parameters
-    name = name_parameters(formdata)
+    name = name_parameters(formdata, define=False)
     define_parameterclass(name)
     return name
 

test/CMakeLists.txt

 add_subdirectory(hyperbolic)
 add_subdirectory(heatequation)
 add_subdirectory(laplace)
+add_subdirectory(navier-stokes)
 add_subdirectory(nonlinear)
 add_subdirectory(poisson)
 add_subdirectory(stokes)
 
 add_subdirectory(sumfact)
 
-add_subdirectory(blockstructured)
\ No newline at end of file
+add_subdirectory(blockstructured)

test/navier-stokes/CMakeLists.txt

+add_subdirectory(reference_program)
+
+dune_add_formcompiler_system_test(UFLFILE navierstokes_2d_dg_quadrilateral.ufl
+                                  BASENAME navierstokes_2d_dg_quadrilateral
+                                  INIFILE navierstokes_2d_dg_quadrilateral.mini
+                                  SCRIPT dune_execute_parallel.py
+                                  )
+
+dune_add_formcompiler_system_test(UFLFILE navierstokes_3d_dg_quadrilateral.ufl
+                                  BASENAME navierstokes_3d_dg_quadrilateral
+                                  INIFILE navierstokes_3d_dg_quadrilateral.mini
+                                  )

test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini

+__name = navierstokes_2d_dg_quadrilateral_{__exec_suffix}
+__exec_suffix = symdiff, numdiff | expand num
+
+cells = 16 16
+lowerleft = -1. -1.
+extension = 2. 2.
+periodic = true true
+
+printmatrix = false
+
+[wrapper.execute_parallel]
+numprocessors = 4
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0, 1 | expand num
+compare_l2errorsquared = 5e-5
+# Only calculate error for the velocity part
+l2error_tree_path = 1, 1, 0
+explicit_time_stepping = 0
+grid_offset = 1
+overlapping = 1
+
+[instat]
+T = 1e-2
+dt = 1e-3
+nth = 1
+
+# Disable numdiff tests
+{__exec_suffix} == numdiff | exclude

test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl

+# Taylor-Green vortex
+
+cell = quadrilateral
+degree = 2
+dim = 2
+
+x = SpatialCoordinate(cell)
+time = get_time(cell)
+
+P2 = VectorElement("DG", cell, degree)
+P1 = FiniteElement("DG", cell, degree-1)
+TH = P2 * P1
+
+v, q = TestFunctions(TH)
+u, p = TrialFunctions(TH)
+
+n = FacetNormal(cell)('+')
+
+rho = 1.0
+mu = 1.0/100.0
+
+g_v = as_vector((-exp(-2*pi*mu/rho*time)*cos(pi*x[0])*sin(pi*x[1]),
+                 exp(-2*pi*mu/rho*time)*sin(pi*x[0])*cos(pi*x[1])))
+g_p = -0.25*rho*exp(-4*pi*pi*mu/rho*time)*(cos(2*pi*x[0])+cos(2*pi*x[1]))
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+alpha = 1.0
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
+mass = rho*inner(u,v)*dx
+
+r = mu * inner(grad(u), grad(v))*dx \
+  - p*div(v)*dx \
+  - q*div(u)*dx \
+  + rho * inner(grad(u)*u,v)*dx \
+  + mu * inner(avg(grad(u))*n, jump(v))*dS \
+  + mu * gamma_int * inner(jump(u), jump(v))*dS \
+  - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
+  - avg(p)*inner(jump(v), n)*dS \
+  - avg(q)*inner(jump(u), n)*dS \
+
+forms = [mass,r]
+interpolate_expression = g_v, g_p
+exact_solution = g_v, g_p

test/navier-stokes/navierstokes_3d_dg_quadrilateral.mini

+__name = navierstokes_3d_dg_quadrilateral_{__exec_suffix}
+__exec_suffix = symdiff, numdiff | expand num
+
+cells = 4 4 4
+lowerleft = -1. -1. -1.
+extension = 2. 2. 2.
+
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0, 1 | expand num
+explicit_time_stepping = 0
+grid_offset = 1
+compare_l2errorsquared = 5e-4
+# Only calculate error for the velocity part
+l2error_tree_path = 1, 1, 1, 0
+
+[instat]
+T = 1e-1
+dt = 5e-2
+
+# Disable numdiff tests
+{__exec_suffix} == numdiff | exclude

test/navier-stokes/navierstokes_3d_dg_quadrilateral.ufl

+# Beltrami flow
+
+cell = hexahedron
+degree = 2
+dim = 2
+
+x = SpatialCoordinate(cell)
+time = get_time(cell)
+
+P2 = VectorElement("DG", cell, degree)
+P1 = FiniteElement("DG", cell, degree-1)
+TH = P2 * P1
+
+v, q = TestFunctions(TH)
+u, p = TrialFunctions(TH)
+
+# g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0, 0.0))
+# bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+# ds = ds(subdomain_id=1, subdomain_data=bctype)
+
+n = FacetNormal(cell)('+')
+
+rho = 1.0
+mu = 1.0
+
+a = pi/4
+d = pi/2
+g_v = as_vector((-a*exp(-d*d*time)*(exp(a*x[0])*sin(d*x[2]+a*x[1])+cos(d*x[1]+a*x[0])*exp(a*x[2])),
+                 -a*exp(-d*d*time)*(exp(a*x[0])*cos(d*x[2]+a*x[1])+exp(a*x[1])*sin(a*x[2]+d*x[0])),
+                 -a*exp(-d*d*time)*(exp(a*x[1])*cos(a*x[2]+d*x[0])+sin(d*x[1]+a*x[0])*exp(a*x[2]))))
+g_p = -0.5*a*a*rho*exp(-d*d*time)*  (  2*cos(d*x[1]+a*x[0])*exp(a*(x[2]+x[0]) )*sin(d*x[2]+a*x[1])  +  2*exp(a*(x[1]+x[0]))*sin(a*x[2]+d*x[0])*cos(d*x[2]+a*x[1])  +  2*sin(d*x[1]+a*x[0])*exp(a*(x[2]+x[1]))*cos(a*x[2]+d*x[0])  +  exp(2*a*x[2])  +  exp(2*a*x[1])  +  exp(2*a*x[0])  )
+
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
+mass = rho*inner(u,v)*dx
+
+r = mu * inner(grad(u), grad(v))*dx \
+  - p*div(v)*dx \
+  - q*div(u)*dx \
+  + rho * inner(grad(u)*u,v)*dx \
+  + mu * inner(avg(grad(u))*n, jump(v))*dS \
+  + mu * gamma_int * inner(jump(u), jump(v))*dS \
+  - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
+  - avg(p)*inner(jump(v), n)*dS \
+  - avg(q)*inner(jump(u), n)*dS \
+  - mu * inner(grad(u)*n, v)*ds \
+  + mu * gamma_ext * inner(u-g_v, v)*ds \
+  + mu * theta * inner(grad(v)*n, u-g_v)*ds \
+  + p*inner(v, n)*ds \
+  + q*inner(u-g_v, n)*ds
+
+# r = mu * inner(grad(u), grad(v))*dx \
+#   - p*div(v)*dx \
+#   - q*div(u)*dx \
+#   + rho * inner(grad(u)*u,v)*dx \
+#   + mu * inner(avg(grad(u))*n, jump(v))*dS \
+#   + mu / h_e * sigma * inner(jump(u), jump(v))*dS \
+#   - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
+#   - avg(p)*inner(jump(v), n)*dS \
+#   - avg(q)*inner(jump(u), n)*dS \
+#   - mu * inner(grad(u)*n, v)*ds \
+#   + mu / h_e * sigma * inner(u-g_v, v)*ds \
+#   + mu * theta * inner(grad(v)*n, u-g_v)*ds \
+#   + p*inner(v, n)*ds \
+#   + q*inner(u-g_v, n)*ds
+
+forms = [mass,r]
+interpolate_expression = g_v, g_p
+exact_solution = g_v, g_p
\ No newline at end of file

test/navier-stokes/reference_program/CMakeLists.txt

+add_executable(taylor-green taylor-green.cc)
+dune_symlink_to_source_files(FILES taylor-green.ini)
+set_target_properties(taylor-green PROPERTIES EXCLUDE_FROM_ALL 1)

test/navier-stokes/reference_program/taylor-green.cc

+// -*- tab-width: 2; indent-tabs-mode: nil -*-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <iostream>
+#include <vector>
+#include <map>
+#include <string>
+#include <random>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/pdelab/common/function.hh>
+#include <dune/pdelab/common/functionutilities.hh>
+#include <dune/pdelab/finiteelementmap/qkdg.hh>
+#include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include <dune/pdelab/gridfunctionspace/subspace.hh>
+#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/gridfunctionspace/interpolate.hh>
+#include <dune/pdelab/localoperator/dgnavierstokes.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/finiteelementmap/monomfem.hh>
+#include <dune/pdelab/common/function.hh>
+#include <dune/pdelab/common/vtkexport.hh>
+#include <dune/pdelab/constraints/p0.hh>
+#include<dune/pdelab/gridoperator/onestep.hh>
+#include<dune/pdelab/newton/newton.hh>
+#include "dune/perftool/vtkpredicate.hh"
+#include "dune/grid/io/file/vtk/vtksequencewriter.hh"
+
+#include "taylor-green.hh"
+
+
+#define PERIODIC
+// #define NORMALIZE_PRESSURE
+
+//===============================================================
+// Problem setup and solution
+//===============================================================
+template<typename GV, typename RF, int vOrder, int pOrder>
+void taylor_green (const GV& gv, const Dune::ParameterTree& configuration, std::string filename)
+{
+  // Some types
+  using ES = Dune::PDELab::AllEntitySet<GV>;
+  ES es(gv);
+  using DF = typename ES::Grid::ctype;
+  static const unsigned int dim = ES::dimension;
+  Dune::Timer timer;
+
+  // Create finite element maps
+  const int velocity_degree = 2;
+  const int pressure_degree = 1;
+  using FEM_V = Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, velocity_degree, dim>;
+  using FEM_P = Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, pressure_degree, dim>;
+  FEM_V fem_v;
+  FEM_P fem_p;
+
+  // Do not block anything and order it lexicographic
+  using VectorBackend_V = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none>;
+  using VectorBackend_P = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none>;
+  using VectorBackend = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none>;
+  using OrderingTag_V = Dune::PDELab::LexicographicOrderingTag;
+
+  // For periodic boundary conditions in Yasp grid we need an
+  // overlap. Therefore we run our program in parallel and need these
+  // constraints
+#ifdef PERIODIC
+  using Con = Dune::PDELab::P0ParallelConstraints;
+#else
+  using Con = Dune::PDELab::NoConstraints;
+#endif
+
+  // Velocity GFS
+  using GFS_V = Dune::PDELab::VectorGridFunctionSpace<
+    ES,FEM_V,dim,
+    VectorBackend,
+    VectorBackend_V,
+    Con,
+    OrderingTag_V
+    >;
+  GFS_V gfs_v(es,fem_v);
+  gfs_v.name("v");
+
+  // Pressure GFS
+  using GFS_P = Dune::PDELab::GridFunctionSpace<
+    ES,
+    FEM_P,
+    Con,
+    VectorBackend_P>;
+  GFS_P gfs_p(es,fem_p);
+  gfs_p.name("p");
+
+
+  // GFS
+  using OrderingTag = Dune::PDELab::LexicographicOrderingTag;
+  using GFS = Dune::PDELab::CompositeGridFunctionSpace<VectorBackend,OrderingTag,GFS_V,GFS_P>;
+  GFS gfs(gfs_v, gfs_p);
+  using namespace Dune::Indices;
+  gfs_v.child(_0).name("velocity_0");
+  gfs_v.child(_1).name("velocity_1");
+  gfs_p.name("pressure");
+  gfs.name("test");
+  gfs.update();
+  using CC = typename GFS::template ConstraintsContainer<double>::Type;
+  CC cc;
+  cc.clear();
+#ifdef PERIODIC
+  Dune::PDELab::constraints(gfs,cc);
+#endif
+  std::cout << "gfs with " << gfs.size() << " dofs generated  "<< std::endl;
+  std::cout << "cc with " << cc.size() << " dofs generated  "<< std::endl;
+
+  // Parameter functions
+  using FType = ZeroVectorFunction<ES,RF,dim>;
+  FType f(es);
+  using BType = BCTypeParamGlobalDirichlet;
+  BType b;
+  using VType = TaylorGreenVelocity<ES,RF,dim>;
+  VType v(es, configuration.sub("parameters"));
+  using PType = TaylorGreenPressure<ES,RF>;
+  PType p(es, configuration.sub("parameters"));
+  using PenaltyTerm = Dune::PDELab::DefaultInteriorPenalty<RF>;
+
+  // Local operator
+  using LOP_Parameters =
+    Dune::PDELab::DGNavierStokesParameters<ES,RF,FType,BType,VType,PType,true,false,PenaltyTerm>;
+  LOP_Parameters lop_parameters(configuration.sub("parameters"),f,b,v,p);
+  using LOP = Dune::PDELab::DGNavierStokes<LOP_Parameters>;
+  const int superintegration_order = 0;
+  LOP lop(lop_parameters,superintegration_order);
+  using LOP_M = Dune::PDELab::NavierStokesMass<LOP_Parameters>;
+  LOP_M lop_m(lop_parameters,1);
+
+  // Grid operator
+  using MBE = Dune::PDELab::istl::BCRSMatrixBackend<>;
+  MBE mbe(75); // Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+  using GO_R = Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,CC,CC>;
+  GO_R go_r(gfs,cc,gfs,cc,lop,mbe);
+  using GO_M = Dune::PDELab::GridOperator<GFS,GFS,LOP_M,MBE,RF,RF,RF,CC,CC>;
+  GO_M go_m(gfs,cc,gfs,cc,lop_m,mbe);
+  using IGO = Dune::PDELab::OneStepGridOperator<GO_R,GO_M>;
+  IGO igo(go_r,go_m);
+
+  // Create initial solution
+  using InitialVelocity = TaylorGreenVelocity<GV,RF,2>;
+  InitialVelocity initial_velocity(gv, configuration.sub("parameters"));
+  using InitialPressure = TaylorGreenPressure<GV,RF>;
+  InitialPressure initial_pressure(gv, configuration.sub("parameters"));
+  using InitialSolution = Dune::PDELab::CompositeGridFunction<InitialVelocity,InitialPressure>;
+  InitialSolution initial_solution(initial_velocity, initial_pressure);
+
+  // Make coefficent vector and initialize it from a function
+  using V = typename IGO::Traits::Domain;
+  V xold(gfs);
+  xold = 0.0;
+  Dune::PDELab::interpolate(initial_solution,gfs,xold);
+
+  // Solver
+#ifdef PERIODIC
+  using LinearSolver = Dune::PDELab::ISTLBackend_OVLP_BCGS_ILU0<GFS,CC>;
+  LinearSolver ls(gfs,cc);
+#else
+  using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_BCGS_ILU0;
+  LinearSolver ls;
+  // using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_UMFPack;
+  // LinearSolver ls(false);
+#endif
+  using PDESolver = Dune::PDELab::Newton<IGO,LinearSolver,V>;
+  PDESolver newton(igo,xold,ls);
+  // newton.setReassembleThreshold(0.0);
+  // newton.setVerbosityLevel(2);
+  // newton.setMaxIterations(50);
+  // newton.setLineSearchMaxIterations(30);
+
+  // Time stepping
+  // using TSM = Dune::PDELab::OneStepThetaParameter<RF>;
+  // TSM tsm(1.0);
+  using TSM = Dune::PDELab::Alexander2Parameter<RF>;
+  TSM tsm;
+  Dune::PDELab::OneStepMethod<RF,IGO,PDESolver,V,V> osm(tsm,igo,newton);
+  // osm.setVerbosityLevel(2);
+
+  // Set time
+  RF time = 0.0;
+  RF time_end = configuration.get<RF>("driver.time_end");
+  RF dt = configuration.get<RF>("driver.dt");
+  RF dt_min = 1e-8;
+
+  // Visualize initial condition
+  using VTKSW = Dune::VTKSequenceWriter<GV>;
+  using VTKWriter = Dune::SubsamplingVTKWriter<GV>;
+  VTKWriter vtkwriter(gv, 2);
+  VTKSW vtkSequenceWriter(std::make_shared<VTKWriter>(vtkwriter), filename);
+  CuttingPredicate predicate;
+  Dune::PDELab::addSolutionToVTKWriter(vtkSequenceWriter, gfs, xold, Dune::PDELab::vtk::defaultNameScheme(), predicate);
+  vtkSequenceWriter.write(time, Dune::VTK::appendedraw);
+
+  V x(gfs,0.0);
+
+#ifdef NORMALIZE_PRESSURE
+  // Pressure normalization
+  using PressureSubGFS = typename Dune::PDELab::GridFunctionSubSpace <GFS,Dune::TypeTree::TreePath<1> >;
+  PressureSubGFS pressureSubGfs(gfs);
+  using PDGF = Dune::PDELab::DiscreteGridFunction<PressureSubGFS,V>;
+  PDGF pdgf(pressureSubGfs,x);
+  typename PDGF::Traits::RangeType pressure_integral(0);
+
+  int elements = int(sqrt(gv.size(0)));
+  int pressure_index = elements * elements * dim * pow((velocity_degree + 1), dim);
+  using Dune::PDELab::Backend::native;
+  std::cout << std::endl;
+  std::cout << "info elements: " << elements << std::endl;
+  std::cout << "info pressure_index: " << pressure_index << std::endl;
+  std::cout << "info gfs.size(): " << gfs.size() << std::endl;
+  std::cout << "info native(x).size(): " << native(x).size() << std::endl;
+  std::cout << std::endl;
+#endif
+
+  // Time loop
+  int step = 0;
+  const int nth = configuration.get<RF>("driver.nth");
+  while (time < time_end - dt_min*0.5){
+    osm.apply(time,dt,xold,x);
+
+#ifdef NORMALIZE_PRESSURE
+    // Correct pressure after each step. Without this pressure
+    // correction the velocity will be ok but the pressure will be
+    // shifted by a constant.
+    Dune::PDELab::integrateGridFunction(pdgf,pressure_integral,2);
+    pressure_integral = gv.comm().sum(pressure_integral);
+    std::cout << gv.comm().rank() << " pressure_integral before normalization: " << pressure_integral << std::endl;
+
+    // Scale integral
+    pressure_integral = pressure_integral/4;
+    for (int i=pressure_index; i<gfs.size(); ++i){
+      native(x)[i] -= pressure_integral;
+    }
+    Dune::PDELab::integrateGridFunction(pdgf,pressure_integral,2);
+    pressure_integral = gv.comm().sum(pressure_integral);
+    std::cout << "pressure_integral after normalization: " << pressure_integral << std::endl;
+#endif
+
+    xold = x;
+    time += dt;
+    step++;
+
+    if(step%nth==0){
+      vtkSequenceWriter.write(time, Dune::VTK::appendedraw);
+    }
+  }
+}
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+int main(int argc, char** argv)
+{
+  try{
+    // Maybe initialize Mpi
+    Dune::MPIHelper::instance(argc, argv);
+
+    // Read ini file
+    Dune::ParameterTree configuration;
+    const std::string config_filename("taylor-green.ini");
+    std::cout << "Reading ini-file \""<< config_filename
+              << "\"" << std::endl;
+    Dune::ParameterTreeParser::readINITree(config_filename, configuration);
+
+    // Create grid
+    const int dim = 2;
+    const int cells_per_dir = configuration.get<double>("driver.cells_per_dir");
+    Dune::FieldVector<double,dim> lowerleft(-1.0);
+    Dune::FieldVector<double,dim> upperright(1.0);
+    std::array<int, dim> cells(Dune::fill_array<int, dim>(cells_per_dir));
+    std::bitset<dim> periodic(false);
+    int overlap = 0;
+#ifdef PERIODIC
+    periodic[0] = true;
+    periodic[1] = true;
+    overlap = 1;
+#endif
+    using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+    Grid grid(lowerleft, upperright, cells, periodic, overlap);
+
+    // Solve problem
+    using GV = Grid::LeafGridView;
+    const GV gv=grid.leafGridView();
+    Dune::dinfo.push(false);
+    taylor_green<GV,double,2,1>(gv,configuration,"taylor-green");
+    return 0;
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+}

test/navier-stokes/reference_program/taylor-green.hh

+#ifndef TAYLOR_GREEN_HH
+#define TAYLOR_GREEN_HH
+
+//===============================================================
+// Define parameter functions f,g,j and \partial\Omega_D/N
+//===============================================================
+
+// constraints parameter class for selecting boundary condition type
+class BCTypeParamGlobalDirichlet
+{
+public :
+  typedef Dune::PDELab::StokesBoundaryCondition BC;
+
+  struct Traits
+  {
+    typedef BC::Type RangeType;
+  };
+
+  BCTypeParamGlobalDirichlet() {}
+
+  template<typename I>
+  inline void evaluate (const I & intersection,   /*@\label{bcp:name}@*/
+                        const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
+                        BC::Type& y) const
+  {
+    y = BC::VelocityDirichlet;
+  }
+
+  template<typename T>
+  void setTime(T t){
+  }
+};
+
+
+template<typename GV, typename RF, int dim>
+class TaylorGreenVelocity :
+  public Dune::PDELab::AnalyticGridFunctionBase<
+  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
+  TaylorGreenVelocity<GV,RF,dim> >
+{
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
+  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TaylorGreenVelocity<GV,RF,dim> > BaseT;
+
+  typedef typename Traits::DomainType DomainType;
+  typedef typename Traits::RangeType RangeType;
+
+  TaylorGreenVelocity(const GV & gv, const Dune::ParameterTree& params) : BaseT(gv)
+  {
+    mu = params.get<RF>("mu");
+    rho = params.get<RF>("rho");
+    time = 0.0;
+  }
+
+  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
+  {
+    // TODO Get mu and rho from somewhere else!
+    RF pi = 3.14159265358979323846;
+    RF nu = mu/rho;
+    y[0] = -exp(-2.0*pi*pi*nu*time)*cos(pi*x[0])*sin(pi*x[1]);
+    y[1] = exp(-2.0*pi*pi*nu*time)*sin(pi*x[0])*cos(pi*x[1]);
+  }
+
+  template <typename T>
+  void setTime(T t){
+    time = t;
+  }
+
+private:
+  RF rho;
+  RF mu;
+  RF time;
+};
+
+
+template<typename GV, typename RF>
+class TaylorGreenPressure
+  : public Dune::PDELab::AnalyticGridFunctionBase<
+  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
+  TaylorGreenPressure<GV,RF> >
+{
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TaylorGreenPressure<GV,RF> > BaseT;
+
+  typedef typename Traits::DomainType DomainType;
+  typedef typename Traits::RangeType RangeType;
+
+  TaylorGreenPressure (const GV& gv, const Dune::ParameterTree& params) : BaseT(gv)
+  {
+    mu = params.get<RF>("mu");
+    rho = params.get<RF>("rho");
+    time = 0.0;
+  }
+
+  inline void evaluateGlobal (const typename Traits::DomainType& x,
+                              typename Traits::RangeType& y) const
+  {
+    RF pi = 3.14159265358979323846;
+    RF nu = mu/rho;
+    y = -0.25*rho*exp(-4.0*pi*pi*nu*time)*(cos(2.0*pi*x[0])+cos(2.0*pi*x[1]));
+  }
+
+  template<typename T>
+  void setTime(T t){
+    time = t;
+  }
+
+private:
+  RF rho;
+  RF mu;
+  RF time;
+};
+
+
+
+template<typename GV, typename RF, std::size_t dim_range>
+class ZeroVectorFunction :
+  public Dune::PDELab::AnalyticGridFunctionBase<
+  Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range>,
+  ZeroVectorFunction<GV,RF,dim_range> >,
+  public Dune::PDELab::InstationaryFunctionDefaults
+{
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range> Traits;
+  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, ZeroVectorFunction> BaseT;
+
+  typedef typename Traits::DomainType DomainType;
+  typedef typename Traits::RangeType RangeType;
+
+  ZeroVectorFunction(const GV & gv) : BaseT(gv) {}
+
+  inline void evaluateGlobal(const DomainType & x, RangeType & y) const
+  {
+    y=0;
+  }
+};
+
+template<typename GV, typename RF>
+class ZeroScalarFunction
+  : public ZeroVectorFunction<GV,RF,1>
+{
+public:
+
+  ZeroScalarFunction(const GV & gv) : ZeroVectorFunction<GV,RF,1>(gv) {}
+
+};
+
+#endif