diff --git a/python/dune/perftool/pdelab/driver/error.py b/python/dune/perftool/pdelab/driver/error.py
index dbbb0b9c4f13327e26a48f8ef8f87f1e1cf75693..6e7fb9128e5437650baf140191311ee4a02eb5f7 100644
--- a/python/dune/perftool/pdelab/driver/error.py
+++ b/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)]
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
index b439ded1b1465dd87fa839329a7238505409c7e5..2ab2234736db0d761efa1b1dfb29298d2b815983 100644
--- a/python/dune/perftool/pdelab/driver/gridfunctionspace.py
+++ b/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))
diff --git a/python/dune/perftool/pdelab/driver/instationary.py b/python/dune/perftool/pdelab/driver/instationary.py
index 7cefddb50d6a5cb2907c228dd328acbdc1db6170..4d8064d6286edf6b5ed59f6bea46832ea093ea1e 100644
--- a/python/dune/perftool/pdelab/driver/instationary.py
+++ b/python/dune/perftool/pdelab/driver/instationary.py
@@ -1,11 +1,35 @@
-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)
diff --git a/python/dune/perftool/pdelab/driver/interpolate.py b/python/dune/perftool/pdelab/driver/interpolate.py
index 626896c11e68495ab0337ea5bcde5a8a993f4a57..3e92f4b9e839f9986829af4967e862b81c612dbc 100644
--- a/python/dune/perftool/pdelab/driver/interpolate.py
+++ b/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)
diff --git a/python/dune/perftool/pdelab/driver/vtk.py b/python/dune/perftool/pdelab/driver/vtk.py
index 15c39c3d22f2a7542c0c08e59cc8ac1c9838ab31..a4adbff493351b7530833bea708114a3257beb5a 100644
--- a/python/dune/perftool/pdelab/driver/vtk.py
+++ b/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