Leaf gfs have to be unique isntances in composite grid function spaces!

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

@@ -9,13 +9,18 @@ from dune.perftool.pdelab.driver import (get_formdata,
                                          get_trial_element,
                                          preprocess_leaf_data,
                                          )
-from dune.perftool.pdelab.driver.gridfunctionspace import name_gfs
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+                                                           name_trial_subgfs,
+                                                           type_range,
+                                                           )
 from dune.perftool.pdelab.driver.interpolate import (interpolate_vector,
                                                      name_boundary_function,
                                                      )
 from dune.perftool.pdelab.driver.solve import (define_vector,
+                                               dune_solve,
                                                name_vector,
                                                )
+from ufl import MixedElement, TensorElement, VectorElement
 
 
 @preamble
@@ -29,43 +34,18 @@ def name_test_fail_variable():
     return name
 
 
-def name_exact_solution():
-    name = "solution"
-    define_vector(name, get_formdata())
-    interpolate_exact_solution(name)
-    return name
-
-
-def interpolate_exact_solution(name):
+def name_exact_solution_gridfunction(treepath):
     element = get_trial_element()
-    is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
-    gfs = name_gfs(element, is_dirichlet)
-    sol = preprocess_leaf_data(element, "exact_solution")
-    func = name_boundary_function(element, sol)
-    interpolate_vector(func, gfs, name)
-
-
-@preamble
-def compare_dofs():
-    v = name_vector(get_formdata())
-    solution = name_exact_solution()
-    fail = name_test_fail_variable()
-
-    return ["",
-            "// Maximum norm of difference between dof vectors of the",
-            "// numerical solution and the interpolation of the exact solution",
-            "using Dune::PDELab::Backend::native;",
-            "double maxerror(0.0);",
-            "for (std::size_t i=0; i<native({}).size(); ++i)".format(v),
-            "  if (std::abs(native({})[i]-native({})[i]) > maxerror)".format(v, solution),
-            "    maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
-            "std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
-            "if (maxerror>{})".format(get_option("compare_dofs")),
-            "  {} = true;".format(fail)]
+    func = preprocess_leaf_data(element, "exact_solution")
+    if isinstance(element, MixedElement):
+        index = treepath_to_index(element, treepath)
+        func = (func[index],)
+        element = element.extract_component(index)[1]
+    return name_boundary_function(element, func)
 
 
 def type_discrete_grid_function(gfs):
-    return "DGF_{}".format(gfs.upper())
+    return "{}_DGF".format(gfs.upper())
 
 
 @preamble
@@ -82,53 +62,93 @@ def name_discrete_grid_function(gfs, vector_name):
 
 
 @preamble
-def typedef_difference_squared_adapter(name):
-    sol = name_exact_solution()
+def typedef_difference_squared_adapter(name, treepath):
+    sol = name_exact_solution_gridfunction(treepath)
     vector = name_vector(get_formdata())
+    gfs = name_trial_subgfs(treepath)
     dgf = name_discrete_grid_function(gfs, vector)
 
     return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, sol, dgf)
 
 
-def type_difference_squared_adapter():
-    name = 'DifferenceSquaredAdapter'
-    typedef_difference_squared_adapter(name)
+def type_difference_squared_adapter(treepath):
+    name = 'DifferenceSquaredAdapter_{}'.format("_".join(str(t) for t in treepath))
+    typedef_difference_squared_adapter(name, treepath)
     return name
 
 
 @preamble
-def define_difference_squared_adapter(name):
-    t = type_difference_squared_adapter()
-    sol = name_exact_solution()
+def define_difference_squared_adapter(name, treepath):
+    t = type_difference_squared_adapter(treepath)
+    sol = name_exact_solution_gridfunction(treepath)
     vector = name_vector(get_formdata())
+    gfs = name_trial_subgfs(treepath)
     dgf = name_discrete_grid_function(gfs, vector)
 
     return '{} {}({}, {});'.format(t, name, sol, dgf)
 
 
-def name_difference_squared_adapter():
-    name = 'differencesquared_adapter'
-    define_difference_squared_adapter(name)
+def name_difference_squared_adapter(treepath):
+    name = 'dsa_{}'.format("_".join(str(t) for t in treepath))
+    define_difference_squared_adapter(name, treepath)
     return name
 
 
 @preamble
-def accumulate_L2_squared():
-    dsa = name_difference_squared_adapter()
+def _accumulate_L2_squared(treepath):
+    dsa = name_difference_squared_adapter(treepath)
     accum_error = name_accumulated_L2_error()
 
     include_file("dune/pdelab/gridfunctionspace/gridfunctionadapter.hh", filetag="driver")
     include_file("dune/pdelab/common/functionutilities.hh", filetag="driver")
 
-    return ["// L2 error squared of difference between numerical",
-            "// solution and the interpolation of exact solution",
-            "  Dune::PDELab::integrateGridFunction({}, {}, 10);".format(dsa, accum_error),
+    strtp = ", ".join(str(t) for t in treepath)
+
+    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),
+            "  {} += err;".format(accum_error),
+            "  std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
+            "}",
             ]
 
 
+def get_treepath(element, index):
+    if isinstance(element, (VectorElement, TensorElement)):
+        return (index,)
+    if isinstance(element, MixedElement):
+        pos, rest = element.extract_subelement_component(index)
+        offset = sum(element.sub_elements()[i].value_size() for i in range(pos))
+        return (pos,) + get_treepath(element.sub_elements()[pos], index - offset)
+    else:
+        return ()
+
+
+def treepath_to_index(element, treepath, offset=0):
+    if len(treepath) == 0:
+        return offset
+    index = treepath[0]
+    offset = offset + sum(element.sub_elements()[i].value_size() for i in range(index))
+    subel = element.sub_elements()[index]
+    return treepath_to_index(subel, treepath[1:], offset)
+
+
+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))
+    else:
+        _accumulate_L2_squared(())
+
+
 @preamble
 def define_accumulated_L2_error(name):
-    return "double {}(0.0);".format(name)
+    t = type_range()
+    return "{} {}(0.0);".format(t, name)
 
 
 def name_accumulated_L2_error():
@@ -143,8 +163,10 @@ def compare_L2_squared():
 
     accum_error = name_accumulated_L2_error()
     fail = name_test_fail_variable()
-    return ["std::cout << \"l2errorsquared: \" << {} << std::endl;".format(accum_error),
-            "if (std::isnan({0}) or std::abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
+    return ["using std::abs;",
+            "using std::isnan;",
+            "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

@@ -175,27 +175,30 @@ def name_test_gfs():
     return name_gfs(element, is_dirichlet)
 
 
-def name_gfs(element, is_dirichlet):
+def name_gfs(element, is_dirichlet, treepath=()):
     if isinstance(element, (VectorElement, TensorElement)):
-        subel = element.sub_element()[0]
-        subgfs = name_gfs(subel, is_dirichlet[:subel.value_size()])
-        name = "{}_pow{}gfs".format(subgfs, element.num_sub_elements())
+        subel = element.sub_elements()[0]
+        subgfs = name_gfs(subel, is_dirichlet[:subel.value_size()], treepath=treepath + (0,))
+        name = "{}_pow{}gfs_{}".format(subgfs,
+                                       element.num_sub_elements(),
+                                       "_".join(str(t) for t in treepath))
         define_power_gfs(element, is_dirichlet, name, subgfs)
         return name
     if isinstance(element, MixedElement):
         k = 0
         subgfs = []
-        for subel in element.sub_elements():
-            subgfs.append(name_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
+        for i, subel in enumerate(element.sub_elements()):
+            subgfs.append(name_gfs(subel, is_dirichlet[k:k + subel.value_size()], treepath=treepath + (i,)))
             k = k + subel.value_size()
         name = "_".join(subgfs)
-        define_composite_gfs(element, is_dirichlet, name, subgfs)
+        name = "{}_{}".format(name, "_".join(str(t) for t in treepath))
+        define_composite_gfs(element, is_dirichlet, name, tuple(subgfs))
         return name
     else:
         assert isinstance(element, FiniteElement)
-        name = "{}{}_gfs".format(FEM_name_mangling(element).lower(),
-                                 "_dirichlet" if is_dirichlet else "",
-                                 )
+        name = "{}{}_gfs_{}".format(FEM_name_mangling(element).lower(),
+                                    "_dirichlet" if is_dirichlet[0] else "",
+                                    "_".join(str(t) for t in treepath))
         define_gfs(element, is_dirichlet, name)
         return name
 
@@ -214,7 +217,7 @@ def type_trial_gfs():
 
 def type_gfs(element, is_dirichlet):
     if isinstance(element, (VectorElement, TensorElement)):
-        subel = element.sub_element()[0]
+        subel = element.sub_elements()[0]
         subgfs = type_gfs(subel, is_dirichlet[:subel.value_size()])
         name = "{}_POW{}GFS".format(subgfs, element.num_sub_elements())
         typedef_power_gfs(element, is_dirichlet, name, subgfs)
@@ -226,12 +229,12 @@ def type_gfs(element, is_dirichlet):
             subgfs.append(type_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
             k = k + subel.value_size()
         name = "_".join(subgfs)
-        typedef_composite_gfs(element, name, subgfs)
+        typedef_composite_gfs(element, name, tuple(subgfs))
         return name
     else:
         assert isinstance(element, FiniteElement)
         name = "{}{}_GFS".format(FEM_name_mangling(element).upper(),
-                                 "_dirichlet" if is_dirichlet else "",
+                                 "_dirichlet" if is_dirichlet[0] else "",
                                  )
         typedef_gfs(element, is_dirichlet, name)
         return name
@@ -242,13 +245,15 @@ def define_gfs(element, is_dirichlet, name):
     gfstype = type_gfs(element, is_dirichlet)
     gv = name_leafview()
     fem = name_fem(element)
-    return "{} {}({}, {});".format(gfstype, name, gv, fem)
+    return ["{} {}({}, {});".format(gfstype, name, gv, fem),
+            "{}.name(\"{}\");".format(name, name)]
 
 
 @preamble
 def define_power_gfs(element, is_dirichlet, name, subgfs):
     gv = name_leafview()
     fem = name_fem(element.sub_elements()[0])
+    gfstype = type_gfs(element, is_dirichlet)
     return "{} {}({}, {});".format(gfstype, name, gv, fem)
 
 
@@ -263,7 +268,7 @@ def typedef_gfs(element, is_dirichlet, name):
     vb = type_vectorbackend(element)
     gv = type_leafview()
     fem = type_fem(element)
-    cass = type_constraintsassembler(bool(is_dirichlet))
+    cass = type_constraintsassembler(bool(is_dirichlet[0]))
     return "using {} = Dune::PDELab::GridFunctionSpace<{}, {}, {}, {}>;".format(name, gv, fem, cass, vb)
 
 
@@ -343,3 +348,30 @@ def type_constraintsassembler(is_dirichlet):
         name = "NoConstraintsAssembler"
         typedef_no_constraintsassembler(name)
         return name
+
+
+def name_trial_subgfs(treepath):
+    if len(treepath) == 0:
+        return name_trial_gfs()
+    else:
+        return name_subgfs(treepath)
+
+
+def name_subgfs(treepath):
+    gfs = name_trial_gfs()
+    name = "{}_{}".format(gfs, "_".join(str(t) for t in treepath))
+    define_subgfs(name, treepath)
+    return name
+
+
+@preamble
+def define_subgfs(name, treepath):
+    t = type_subgfs(treepath)
+    gfs = name_trial_gfs()
+    return "{} {}({});".format(t, name, gfs)
+
+
+def type_subgfs(tree_path):
+    include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
+    gfs = type_trial_gfs()
+    return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(str(t) for t in tree_path))

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

-from dune.perftool.generation import preamble
+from dune.perftool.generation import (include_file,
+                                      preamble,
+                                      )
 from dune.perftool.pdelab.driver import (get_formdata,
+                                         get_mass_formdata,
+                                         get_trial_element,
                                          is_linear,
+                                         name_initree,
+                                         preprocess_leaf_data,
                                          )
-from dune.perftool.pdelab.driver.gridoperator import (type_gridoperator,)
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+                                                           type_range,
+                                                           )
+from dune.perftool.pdelab.driver.gridoperator import (name_gridoperator,
+                                                      name_parameters,
+                                                      type_gridoperator,)
+from dune.perftool.pdelab.driver.constraints import (name_bctype_function,
+                                                     name_constraintscontainer,
+                                                     )
+from dune.perftool.pdelab.driver.interpolate import (interpolate_dirichlet_data,
+                                                     name_boundary_function,
+                                                     )
 from dune.perftool.pdelab.driver.solve import (print_matrix,
                                                print_residual,
+                                               name_stationarynonlinearproblemsolver,
+                                               name_vector,
+                                               type_stationarynonlinearproblemssolver,
+                                               type_vector,
                                                )
+from dune.perftool.pdelab.driver.vtk import (name_vtk_sequence_writer,
+                                             visualize_initial_condition,
+                                             )
 from dune.perftool.options import get_option
 
 
@@ -28,12 +52,6 @@ def solve_instationary():
 
     print_residual()
     print_matrix()
-    from dune.perftool.pdelab.driver.error import compare_dofs, compare_L2_squared
-    if get_option("exact_solution_expression"):
-        if get_option("compare_dofs"):
-            compare_dofs()
-        if get_option("compare_l2errorsquared"):
-            compare_L2_squared()
 
 
 @preamble
@@ -42,12 +60,14 @@ def time_loop():
     formdata = get_formdata()
     params = name_parameters(formdata)
     time = name_time()
-    expr = get_trial_element()
-    bctype = name_bctype_function(expr)
-    gfs = name_gfs(expr)
-    cc = name_constraintscontainer(expr)
+    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)
 
     # Choose between explicit and implicit time stepping
     explicit = get_option('explicit_time_stepping')
@@ -55,7 +75,8 @@ def time_loop():
         osm = name_explicitonestepmethod()
         apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
     else:
-        boundary = name_boundary_function(expr)
+        dirichlet = preprocess_leaf_data(element, "dirichlet_expression")
+        boundary = name_boundary_function(element, dirichlet)
         osm = name_onestepmethod()
         apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
 
@@ -95,7 +116,6 @@ def name_time():
     return "time"
 
 
-
 @preamble
 def typedef_timesteppingmethod(name):
     r_type = type_range()
@@ -131,8 +151,8 @@ def name_timesteppingmethod():
 @preamble
 def typedef_instationarygridoperator(name):
     include_file("dune/pdelab/gridoperator/onestep.hh", filetag="driver")
-    go_type = type_gridoperator(_driver_data['formdata'])
-    mass_go_type = type_gridoperator(_driver_data['mass_formdata'])
+    go_type = type_gridoperator(get_formdata())
+    mass_go_type = type_gridoperator(get_mass_formdata())
     explicit = get_option('explicit_time_stepping')
     if explicit:
         return "using {} = Dune::PDELab::OneStepGridOperator<{},{},false>;".format(name, go_type, mass_go_type)
@@ -148,8 +168,8 @@ def type_instationarygridoperator():
 @preamble
 def define_instationarygridoperator(name):
     igo_type = type_instationarygridoperator()
-    go = name_gridoperator(_driver_data['formdata'])
-    mass_go = name_gridoperator(_driver_data['mass_formdata'])
+    go = name_gridoperator(get_formdata())
+    mass_go = name_gridoperator(get_mass_formdata())
     return "{} {}({}, {});".format(igo_type, name, go, mass_go)
 
 
@@ -163,7 +183,7 @@ def typedef_onestepmethod(name):
     r_type = type_range()
     igo_type = type_instationarygridoperator()
     snp_type = type_stationarynonlinearproblemssolver(igo_type)
-    vector_type = type_vector(_driver_data['formdata'])
+    vector_type = type_vector(get_formdata())
     return "using {} = Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}>;".format(name, r_type, igo_type, snp_type, vector_type, vector_type)
 
 
@@ -192,7 +212,7 @@ def typedef_explicitonestepmethod(name):
     r_type = type_range()
     igo_type = type_instationarygridoperator()
     ls_type = type_linearsolver()
-    vector_type = type_vector(_driver_data['formdata'])
+    vector_type = type_vector(get_formdata())
     return "using {} = Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}>;".format(name, r_type, igo_type, ls_type, vector_type)
 
 

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

@@ -10,7 +10,7 @@ from dune.perftool.pdelab.driver import (FEM_name_mangling,
                                          is_stationary,
                                          preprocess_leaf_data,
                                          )
-from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
                                                            name_leafview,
                                                            )
 from dune.perftool.pdelab.driver.gridoperator import (name_parameters,)
@@ -28,8 +28,8 @@ def _do_interpolate(dirichlet):
 def interpolate_dirichlet_data(name):
     element = get_trial_element()
     is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
-    if _do_interpolate(is_dirichlet):
-        gfs = name_gfs(element, is_dirichlet)
+    if _do_interpolate(is_dirichlet) or not is_stationary():
+        gfs = name_trial_gfs()
         dirichlet = preprocess_leaf_data(element, "dirichlet_expression")
         bf = name_boundary_function(element, dirichlet)
         interpolate_vector(bf, gfs, name)

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

@@ -122,7 +122,6 @@ def visualize_initial_condition():
     include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
     vtkwriter = name_vtk_sequence_writer()
     element = get_trial_element()
-    define_gfs_names(element)
     gfs = name_trial_gfs()
     vector = name_vector(get_formdata())
     predicate = name_predicate()
@@ -130,23 +129,3 @@ def visualize_initial_condition():
     time = name_time()
     return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vector, predicate),
             "{}.write({}, Dune::VTK::appendedraw);".format(vtkwriter, time)]
-
-
-@preamble
-def _define_gfs_name(element, is_dirichlet, offset=0):
-    from ufl import MixedElement
-    if isinstance(element, MixedElement):
-        k = 0
-        for subel in element.sub_elements():
-            define_gfs_name(subel, is_dirichlet[k:k + subel.value_size()], offset=offset + k)
-            k = k + subel.value_size()
-        return []
-    else:
-        # This is a scalar leaf
-        gfs = name_gfs(element, is_dirichlet)
-        return "{}.name(\"data_{}\");".format(gfs, offset)
-
-
-def define_gfs_names(element):
-    is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
-    return _define_gfs_name(element, is_dirichlet)
\ No newline at end of file