From 0868a1fd0ce89e058bfb99c06e175d0af87ab588 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Wed, 2 Aug 2017 14:09:50 +0200
Subject: [PATCH] Rewrite driver generation
Do not annotate FiniteElements anymore, further split driver code
---
.../dune/perftool/pdelab/driver/__init__.py | 1331 +----------------
.../perftool/pdelab/driver/constraints.py | 144 ++
python/dune/perftool/pdelab/driver/error.py | 195 +--
.../pdelab/driver/gridfunctionspace.py | 345 +++++
.../perftool/pdelab/driver/gridoperator.py | 173 +++
.../perftool/pdelab/driver/instationary.py | 213 +++
.../perftool/pdelab/driver/interpolate.py | 115 ++
python/dune/perftool/pdelab/driver/solve.py | 245 +++
python/dune/perftool/pdelab/driver/vtk.py | 151 ++
python/dune/perftool/ufl/execution.py | 58 -
test/poisson/poisson.mini | 1 -
test/poisson/poisson.ufl | 5 +-
12 files changed, 1503 insertions(+), 1473 deletions(-)
create mode 100644 python/dune/perftool/pdelab/driver/constraints.py
create mode 100644 python/dune/perftool/pdelab/driver/gridfunctionspace.py
create mode 100644 python/dune/perftool/pdelab/driver/gridoperator.py
create mode 100644 python/dune/perftool/pdelab/driver/instationary.py
create mode 100644 python/dune/perftool/pdelab/driver/interpolate.py
create mode 100644 python/dune/perftool/pdelab/driver/solve.py
create mode 100644 python/dune/perftool/pdelab/driver/vtk.py
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index dcc3ec4a..c59330ce 100644
--- a/python/dune/perftool/pdelab/driver/__init__.py
+++ b/python/dune/perftool/pdelab/driver/__init__.py
@@ -47,6 +47,10 @@ def get_dimension():
return _driver_data['form'].ufl_cell().geometric_dimension()
+def get_cell():
+ return _driver_data['form'].ufl_cell().cellname()
+
+
def get_test_element():
return _driver_data['form'].arguments()[0].ufl_element()
@@ -70,29 +74,8 @@ def form_name_suffix(name, formdata):
return name + '_' + form_name
-def has_constraints(element):
- from ufl import MixedElement, VectorElement
- if isinstance(element, MixedElement) and not isinstance(element, VectorElement):
- return any(has_constraints(el) for el in element.sub_elements())
- return hasattr(element, "dirichlet_constraints")
-
-
-# Those symbol generators that depend on the meta data we hacked into
-# the ufl.FiniteElement should use fem_metadata_dependent_symbol instead of symbol.
-def get_constraints_metadata(element):
- """
- Normally, we suppress the constraints information on a finite element.
- However, we might want to use it explicitly in those case where it matters
- This function does the job!
- """
- dc = getattr(element, 'dirichlet_constraints', None)
- de = getattr(element, 'dirichlet_expression', None)
- return (element, (dc, de))
-
-
-fem_metadata_dependent_symbol = generator_factory(item_tags=("symbol",),
- cache_key_generator=get_constraints_metadata,
- )
+def get_object(name):
+ return _driver_data['data'].object_by_name.get(name, None)
def mass_form_index(formdatas, data):
@@ -104,8 +87,9 @@ def mass_form_index(formdatas, data):
continue
-def is_linear(form):
- '''Test if form is linear in test function'''
+def is_linear():
+ '''Test if form is linear in trial function'''
+ form = get_formdata().original_form
from ufl import derivative
from ufl.algorithms import expand_derivatives
jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
@@ -120,20 +104,26 @@ def isLagrange(fem):
return fem._short_name is 'CG'
-def isSimplical(fem):
- return any(fem._cell._cellname in x for x in ["triangle", "tetrahedron"])
+def isSimplical(cell):
+ # Cells can be identified through strings *or* ufl objects
+ if not isinstance(cell, str):
+ cell = cell.cellname()
+ return cell in ["vertex", "interval", "triangle", "tetrahedron"]
-def isQuadrilateral(fem):
- return any(fem._cell._cellname in x for x in ["vertex", "interval", "quadrilateral", "hexahedron"])
+def isQuadrilateral(cell):
+ # Cells can be identified through strings *or* ufl objects
+ if not isinstance(cell, str):
+ cell = cell.cellname()
+ return cell in ["vertex", "interval", "quadrilateral", "hexahedron"]
def isPk(fem):
- return isLagrange(fem) and isSimplical(fem)
+ return isLagrange(fem) and isSimplical(fem.cell())
def isQk(fem):
- return isLagrange(fem) and isQuadrilateral(fem)
+ return isLagrange(fem) and isQuadrilateral(fem.cell())
def isDG(fem):
@@ -164,6 +154,40 @@ def FEM_name_mangling(fem):
raise NotImplementedError("FEM NAME MANGLING")
+def _flatten_list(l):
+ if isinstance(l, (tuple, list)):
+ for i in l:
+ for ni in _flatten_list(i):
+ yield ni
+ else:
+ yield l
+
+
+def preprocess_leaf_data(element, data):
+ data = get_object(data)
+ from ufl import MixedElement
+ if isinstance(element, MixedElement):
+ # Dirichlet is None -> no dirichlet boundaries
+ if data is None:
+ return (0,) * element.value_size()
+ # Dirichlet for MixedElement is not iterable -> Same
+ # constraint on all the leafs.
+ elif not isinstance(data, (tuple, list)):
+ return (data,) * element.value_size()
+ # List sizes do not match -> flatten list
+ elif len(data) != element.value_size():
+ flattened = [i for i in _flatten_list(data)]
+ assert len(flattened) == element.value_size()
+ return flattened
+ else:
+ # Do not return lists for non-MixedElement
+ if not isinstance(data, (tuple, list)):
+ return (data,)
+ else:
+ assert len(data) == 1
+ return data
+
+
def name_inifile():
# TODO pass some other option here.
return "argv[1]"
@@ -183,945 +207,6 @@ def name_initree():
return "initree"
-@preamble
-def define_dimension(name):
- return "static const int {} = {};".format(name, get_dimension())
-
-
-def name_dimension():
- define_dimension("dim")
- return "dim"
-
-
-@preamble
-def typedef_grid(name):
- dim = name_dimension()
- if any(_driver_data['form'].ufl_cell().cellname() in x for x in ["vertex", "interval", "quadrilateral", "hexahedron"]):
- # For Yasp Grids the jacobi of the transformation is diagonal and constant on each cell
- from dune.perftool.options import set_option
- set_option('diagonal_transformation_matrix', True)
- set_option('constant_transformation_matrix', True)
-
- range_type = type_range()
- gridt = "Dune::YaspGrid<{}, Dune::EquidistantCoordinates<{},dim>>".format(dim, range_type)
- include_file("dune/grid/yaspgrid.hh", filetag="driver")
- else:
- if any(_driver_data['form'].ufl_cell().cellname() in x for x in ["triangle", "tetrahedron"]):
- # gridt = "Dune::UGGrid<{}>".format(dim)
- # include_file("dune/grid/uggrid.hh", filetag="driver")
- gridt = "Dune::ALUGrid<{}, {}, Dune::simplex, Dune::conforming>".format(dim, dim)
- include_file("dune/alugrid/grid.hh", filetag="driver")
- else:
- raise PerftoolCodegenError("Cant match your geometry with a DUNE grid. Please report bug.")
- return "using {} = {};".format(name, gridt)
-
-
-def type_grid():
- typedef_grid("Grid")
- return "Grid"
-
-
-@preamble
-def define_grid(name):
- include_file("dune/testtools/gridconstruction.hh", filetag="driver")
- ini = name_initree()
- _type = type_grid()
- return ["IniGridFactory<{}> factory({});".format(_type, ini),
- "std::shared_ptr<{}> grid = factory.getGrid();".format(_type)]
-
-
-def name_grid():
- define_grid("grid")
- return "grid"
-
-
-@preamble
-def typedef_leafview(name):
- grid = type_grid()
- return "using {} = {}::LeafGridView;".format(name, grid)
-
-
-def type_leafview():
- typedef_leafview("GV")
- return "GV"
-
-
-@preamble
-def define_leafview(name):
- _type = type_leafview()
- grid = name_grid()
- return "{} {} = {}->leafGridView();".format(_type, name, grid)
-
-
-def name_leafview():
- define_leafview("gv")
- return "gv"
-
-
-@preamble
-def typedef_vtkwriter(name):
- include_file("dune/grid/io/file/vtk/subsamplingvtkwriter.hh", filetag="driver")
- gv = type_leafview()
- return "using {} = Dune::SubsamplingVTKWriter<{}>;".format(name, gv)
-
-
-def type_vtkwriter():
- typedef_vtkwriter("VTKWriter")
- return "VTKWriter"
-
-
-@preamble
-def define_subsamplinglevel(name):
- ini = name_initree()
- return "int {} = {}.get<int>(\"vtk.subsamplinglevel\", 0);".format(name, ini)
-
-
-def name_subsamplinglevel():
- define_subsamplinglevel("sublevel")
- return "sublevel"
-
-
-@preamble
-def define_vtkwriter(name):
- _type = type_vtkwriter()
- gv = name_leafview()
- subsamp = name_subsamplinglevel()
- return "{} {}({}, {});".format(_type, name, gv, subsamp)
-
-
-def name_vtkwriter():
- define_vtkwriter("vtkwriter")
- return "vtkwriter"
-
-
-@preamble
-def typedef_domainfield(name):
- gridt = type_grid()
- return "using {} = {}::ctype;".format(name, gridt)
-
-
-def type_domainfield():
- typedef_domainfield("DF")
- return "DF"
-
-
-def basetype_range():
- if get_option('opcounter'):
- from dune.perftool.pdelab.driver.timings import setup_timer
- setup_timer()
- return "oc::OpCounter<double>"
- else:
- return "double"
-
-
-@preamble
-def typedef_range(name):
- return "using {} = {};".format(name, basetype_range())
-
-
-def type_range():
- typedef_range("RangeType")
- return "RangeType"
-
-
-@preamble
-def typedef_fem(expr, name):
- gv = type_leafview()
- df = type_domainfield()
- r = type_range()
- dim = name_dimension()
- if get_option("blockstructured"):
- include_file("dune/perftool/blockstructured/blockstructuredqkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, expr._degree)
- if isPk(expr):
- include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, expr._degree)
- if isQk(expr):
- include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, expr._degree)
- if isDG(expr):
- if isQuadrilateral(expr):
- include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
- # TODO allow switching the basis here!
- return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, df, r, expr._degree, dim)
- if isSimplical(expr):
- include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;".format(name, df, r, expr._degree, dim)
- raise NotImplementedError("Geometry type not known in code generation")
-
- raise NotImplementedError("FEM not implemented in dune-perftool")
-
-
-def type_fem(expr):
- name = "{}_FEM".format(FEM_name_mangling(expr).upper())
- typedef_fem(expr, name)
- return name
-
-
-@preamble
-def define_fem(expr, name):
- femtype = type_fem(expr)
- if isDG(expr):
- return "{} {};".format(femtype, name)
- else:
- gv = name_leafview()
- return "{} {}({});".format(femtype, name, gv)
-
-
-def name_fem(expr):
- name = "{}_fem".format(FEM_name_mangling(expr).lower())
- define_fem(expr, name)
- return name
-
-
-@preamble
-def define_blocksize(name, expr):
- assert isDG(expr)
- assert isQuadrilateral(expr)
- dimension = name_dimension()
- degree = expr._degree
- return "static const int {} = Dune::QkStuff::QkSize<{}, {}>::value;".format(name, degree, dimension)
-
-
-def name_blocksize(expr):
- name = "blocksize"
- define_blocksize(name, expr)
- return name
-
-
-@preamble
-def typedef_vectorbackend(name, expr):
- include_file("dune/pdelab/backend/istl.hh", filetag="driver")
- if get_option("fastdg"):
- blocksize = name_blocksize(expr)
- return "using {} = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed, {}>;".format(name, blocksize)
- else:
- return "using {} = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none, 1>;".format(name)
-
-
-def type_vectorbackend(expr):
- name = "VectorBackend"
- typedef_vectorbackend(name, expr)
- return name
-
-
-def type_orderingtag():
- return "Dune::PDELab::LexicographicOrderingTag"
-
-
-@preamble
-def typedef_dirichlet_constraintsassembler(name):
- include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
- return "using {} = Dune::PDELab::ConformingDirichletConstraints;".format(name)
-
-
-@preamble
-def typedef_no_constraintsassembler(name):
- return "using {} = Dune::PDELab::NoConstraints;".format(name)
-
-
-def type_constraintsassembler(dirichlet):
- if dirichlet[1]:
- typedef_dirichlet_constraintsassembler("DirichletConstraintsAssember")
- return "DirichletConstraintsAssember"
- else:
- typedef_no_constraintsassembler("NoConstraintsAssembler")
- return "NoConstraintsAssembler"
-
-
-@preamble
-def typedef_constraintscontainer(element, name):
- gfs = type_gfs(element)
- r = type_range()
- return "using {} = {}::ConstraintsContainer<{}>::Type;".format(name, gfs, r)
-
-
-def type_constraintscontainer(element):
- name = "{}_cc".format(FEM_name_mangling(element)).upper()
- typedef_constraintscontainer(element, name)
- return name
-
-
-@preamble
-def define_constraintscontainer(expr, name):
- cctype = type_constraintscontainer(expr)
- return ["{} {};".format(cctype, name), "{}.clear();".format(name)]
-
-
-# TODO this does not yet fully support having same elements with different constraints!
-@fem_metadata_dependent_symbol
-def name_constraintscontainer(expr):
- name = "{}{}_cc".format(FEM_name_mangling(expr), "_dirichlet" if has_constraints(expr) else "").lower()
- define_constraintscontainer(expr, name)
- return name
-
-
-@preamble
-def define_intersection_lambda(expression, name):
- from dune.perftool.ufl.execution import Expression
- from ufl.classes import Expr
- if expression is None:
- return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
- elif isinstance(expression, Expression):
- if expression.is_global:
- return "auto {} = [&](const auto& x){{ {} }};".format(name, expression.c_expr[0])
- else:
- return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr[0])
- elif isinstance(expression, Expr):
- # Set up a visitor
- with global_context(integral_type="exterior_facet", formdata=_driver_data["formdata"], driver=True):
- from dune.perftool.ufl.visitor import UFL2LoopyVisitor
- from dune.perftool.pdelab import PDELabInterface
- visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
- from pymbolic.mapper.c_code import CCodeMapper
- ccm = CCodeMapper()
- expr = visitor(expression)
- return "auto {} = [&](const auto& x){{ return {}; }};".format(name, ccm(expr))
-
- raise ValueError("Expression not understood")
-
-
-def name_bctype_lambda(name, dirichlet):
- name = name + "_lambda"
- define_intersection_lambda(dirichlet, name)
- return name
-
-
-@preamble
-def define_bctype_function(dirichlet, name):
- gv = name_leafview()
- bctype_lambda = name_bctype_lambda(name, dirichlet)
- include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
- return "auto {} = Dune::PDELab::makeBoundaryConditionFromCallable({}, {});".format(name,
- gv,
- bctype_lambda,
- )
-
-
-@preamble
-def define_powergfs_constraints(name, child, dim):
- include_file('dune/pdelab/constraints/common/constraintsparameters.hh', filetag='driver')
- return "Dune::PDELab::PowerConstraintsParameters<decltype({}), {}> {}({});".format(child, dim, name, child)
-
-
-@preamble
-def define_compositegfs_constraints(name, *children):
- include_file('dune/pdelab/constraints/common/constraintsparameters.hh', filetag='driver')
- return "Dune::PDELab::CompositeConstraintsParameters<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in children),
- name,
- ', '.join(c for c in children)
- )
-
-
-@cached
-def name_bctype_function(expr):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(expr, (VectorElement, TensorElement)):
- element, (dirichlet, _) = get_constraints_metadata(expr)
- # get the correct name from the corresponding UFL file
- from dune.perftool.generation import get_global_context_value
- name = get_global_context_value("data").object_names.get(id(dirichlet), "everywhere")
- define_bctype_function(dirichlet, name)
- pgfs_name = '{}_{}'.format(name, expr.num_sub_elements())
- define_powergfs_constraints(pgfs_name, name, expr.num_sub_elements())
- return pgfs_name
- if isinstance(expr, MixedElement):
- children = tuple(name_bctype_function(el) for el in expr.sub_elements())
- name = '_'.join(children)
- define_compositegfs_constraints(name, *children)
- return name
-
- # So, this seems to be a leaf!
- element, (dirichlet, _) = get_constraints_metadata(expr)
-
- # get the correct name from the corresponding UFL file
- from dune.perftool.generation import get_global_context_value
- name = get_global_context_value("data").object_names.get(id(dirichlet), "everywhere")
-
- define_bctype_function(dirichlet, name)
- return name
-
-
-@preamble
-def assemble_constraints(name, expr):
- bctype_function = name_bctype_function(expr)
- gfs = name_gfs(expr)
- return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
- gfs,
- name,
- )
-
-
-def name_assembled_constraints(expr):
- name = name_constraintscontainer(expr)
- if has_constraints(expr):
- assemble_constraints(name, expr)
- return name
-
-
-@preamble
-def typedef_gfs(element, dirichlet, name):
- vb = type_vectorbackend(element)
- from ufl import FiniteElement, MixedElement, VectorElement, EnrichedElement, RestrictedElement, TensorElement
- if isinstance(element, FiniteElement):
- gv = type_leafview()
- fem = type_fem(element)
- cass = type_constraintsassembler(dirichlet)
- return "using {} = Dune::PDELab::GridFunctionSpace<{}, {}, {}, {}>;".format(name, gv, fem, cass, vb)
- if isinstance(element, (VectorElement, TensorElement)):
- include_file("dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh", filetag="driver")
- gv = type_leafview()
- fem = type_fem(element.sub_elements()[0])
- # TODO implement symmetry here
- dim = element.num_sub_elements()
- cass = type_constraintsassembler(dirichlet)
- return "using {} = Dune::PDELab::VectorGridFunctionSpace<{}, {}, {}, {}, {}, {}>;".format(name, gv, fem, dim, vb, vb, cass)
- if isinstance(element, MixedElement):
- ot = type_orderingtag()
- args = ", ".join(type_gfs(e) for e in element.sub_elements())
- return "using {} = Dune::PDELab::CompositeGridFunctionSpace<{}, {}, {}>;".format(name, vb, ot, args)
- if isinstance(element, EnrichedElement):
- raise NotImplementedError("Dune does not support enriched elements!")
- if isinstance(element, RestrictedElement):
- raise NotImplementedError("Dune does not support restricted elements!")
-
-
-@fem_metadata_dependent_symbol
-def type_gfs(expr):
- element, dirichlet = get_constraints_metadata(expr)
- name = "{}{}_gfs".format(FEM_name_mangling(element), "_dirichlet" if dirichlet else "").upper()
- typedef_gfs(element, dirichlet, name)
- return name
-
-
-@preamble
-def define_gfs(element, dirichlet, name):
- gfstype = type_gfs(element)
- from ufl import FiniteElement, MixedElement, VectorElement, EnrichedElement, RestrictedElement, TensorElement
- if isinstance(element, FiniteElement):
- gv = name_leafview()
- fem = name_fem(element)
- return "{} {}({}, {});".format(gfstype, name, gv, fem)
- if isinstance(element, (VectorElement, TensorElement)):
- gv = name_leafview()
- fem = name_fem(element.sub_elements()[0])
- return "{} {}({}, {});".format(gfstype, name, gv, fem)
- if isinstance(element, MixedElement):
- args = ", ".join(name_gfs(childgfs) for childgfs in element.sub_elements())
- return "{} {}({});".format(gfstype, name, args)
- raise NotImplementedError("Dune does not support this type of finite element!")
-
-
-@fem_metadata_dependent_symbol
-def name_gfs(expr):
- element, dirichlet = get_constraints_metadata(expr)
- name = "{}{}_gfs".format(FEM_name_mangling(element), "_dirichlet" if dirichlet else "").lower()
- define_gfs(element, dirichlet, name)
- return name
-
-
-# TODO adjust for blockstructured
-@preamble
-def define_dofestimate(name):
- # Provide a worstcase estimate for the number of entries per row based on the given gridfunction space and cell geometry
- if isQuadrilateral(get_trial_element()):
- geo_factor = "4"
- else:
- geo_factor = "6"
- gfs = name_gfs(get_trial_element())
- ini = name_initree()
-
- # Assure that gfs in initialized
- formdata = _driver_data['formdata']
- x = name_vector(formdata)
- define_vector(x, formdata)
-
- return ["int generic_dof_estimate = {} * {}.maxLocalSize();".format(geo_factor, gfs),
- "int {} = {}.get<int>(\"istl.number_of_nnz\", generic_dof_estimate);".format(name, ini)]
-
-
-def name_dofestimate():
- define_dofestimate("dofestimate")
- return "dofestimate"
-
-
-@preamble
-def typedef_matrixbackend(name):
- include_file("dune/pdelab/backend/istl.hh", filetag="driver")
- return "using {} = Dune::PDELab::istl::BCRSMatrixBackend<>;".format(name)
-
-
-def type_matrixbackend():
- typedef_matrixbackend("MatrixBackend")
- return "MatrixBackend"
-
-
-@preamble
-def define_matrixbackend(name):
- mbtype = type_matrixbackend()
- dof = name_dofestimate()
- return "{} {}({});".format(mbtype, name, dof)
-
-
-def name_matrixbackend():
- define_matrixbackend("mb")
- return "mb"
-
-
-def type_parameters(formdata):
- from dune.perftool.generation import get_global_context_value
- data = get_global_context_value("data")
- from dune.perftool.pdelab.parameter import parameterclass_basename
- name = parameterclass_basename(formdata, data)
- return name
-
-
-@preamble
-def define_parameters(name, formdata):
- partype = type_parameters(formdata)
- return "{} {};".format(partype, name)
-
-
-def name_parameters(formdata):
- name = form_name_suffix("params", formdata)
- define_parameters(name, formdata)
- return name
-
-
-@preamble
-def typedef_localoperator(name, formdata):
- # No Parameter class here, yet
- # params = type_parameters()
- # return "using {} = LocalOperator<{}>;".format(name, params)
- ugfs = type_gfs(get_trial_element())
- vgfs = type_gfs(get_test_element())
- from dune.perftool.generation import get_global_context_value
- data = get_global_context_value("data")
- filename = get_option("operator_file")
- include_file(filename, filetag="driver")
- from dune.perftool.pdelab.localoperator import localoperator_basename
- lopname = localoperator_basename(formdata, data)
- range_type = type_range()
- return "using {} = {}<{}, {}, {}>;".format(name, lopname, ugfs, vgfs, range_type)
-
-
-def type_localoperator(formdata):
- name = form_name_suffix("LOP", formdata).upper()
- typedef_localoperator(name, formdata)
- return name
-
-
-@preamble
-def define_localoperator(name, formdata):
- loptype = type_localoperator(formdata)
- ini = name_initree()
- params = name_parameters(formdata)
- return "{} {}({}, {});".format(loptype, name, ini, params)
-
-
-def name_localoperator(formdata):
- name = form_name_suffix("lop", formdata)
- define_localoperator(name, formdata)
- return name
-
-
-@preamble
-def typedef_gridoperator(name, formdata):
- ugfs = type_gfs(get_trial_element())
- vgfs = type_gfs(get_test_element())
- lop = type_localoperator(formdata)
- ucc = type_constraintscontainer(get_trial_element())
- vcc = type_constraintscontainer(get_test_element())
- mb = type_matrixbackend()
- df = type_domainfield()
- r = type_range()
- if get_option("fastdg"):
- if not get_option("sumfact"):
- raise PerftoolCodegenError("FastDGGridOperator is only implemented for sumfactorization.")
- include_file("dune/pdelab/gridoperator/fastdg.hh", filetag="driver")
- return "using {} = Dune::PDELab::FastDGGridOperator<{}, {}, {}, {}, {}, {}, {}, {}, {}>;".format(name, ugfs, vgfs, lop, mb, df, r, r, ucc, vcc)
- else:
- include_file("dune/pdelab/gridoperator/gridoperator.hh", filetag="driver")
- return "using {} = Dune::PDELab::GridOperator<{}, {}, {}, {}, {}, {}, {}, {}, {}>;".format(name, ugfs, vgfs, lop, mb, df, r, r, ucc, vcc)
-
-
-def type_gridoperator(formdata):
- name = form_name_suffix("GO", formdata).upper()
- typedef_gridoperator(name, formdata)
- return name
-
-
-@preamble
-def define_gridoperator(name, formdata):
- gotype = type_gridoperator(formdata)
- # TODO use formdata insteat of _driver_data object
- ugfs = name_gfs(get_trial_element())
- ucc = name_assembled_constraints(get_trial_element())
- vgfs = name_gfs(get_test_element())
- vcc = name_assembled_constraints(get_test_element())
- lop = name_localoperator(formdata)
- mb = name_matrixbackend()
- return ["{} {}({}, {}, {}, {}, {}, {});".format(gotype, name, ugfs, ucc, vgfs, vcc, lop, mb),
- "std::cout << \"gfs with \" << {}.size() << \" dofs generated \"<< std::endl;".format(ugfs),
- "std::cout << \"cc with \" << {}.size() << \" dofs generated \"<< std::endl;".format(ucc)]
-
-
-def name_gridoperator(formdata):
- name = form_name_suffix("go", formdata)
- define_gridoperator(name, formdata)
- return name
-
-
-@preamble
-def typedef_vector(name, formdata):
- go_type = type_gridoperator(formdata)
- return "using {} = {}::Traits::Domain;".format(name, go_type)
-
-
-def type_vector(formdata):
- name = form_name_suffix("V", formdata).upper()
- typedef_vector(name, formdata)
- return name
-
-
-@preamble
-def define_vector(name, formdata):
- vtype = type_vector(formdata)
- gfs = name_gfs(get_trial_element())
- return ["{} {}({});".format(vtype, name, gfs), "{} = 0.0;".format(name)]
-
-
-@preamble
-def define_boundary_lambda(boundary, name):
- from dune.perftool.ufl.execution import Expression
- from ufl.classes import Expr
- if boundary is None:
- return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
- elif isinstance(boundary, Expression):
- if boundary.is_global:
- return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr[0])
- else:
- return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr[0])
- elif isinstance(boundary, Expr):
- # Set up a visitor
- with global_context(integral_type="exterior_facet", formdata=_driver_data["formdata"], driver=True):
- from dune.perftool.ufl.visitor import UFL2LoopyVisitor
- from dune.perftool.pdelab import PDELabInterface
- visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
- from dune.perftool.loopy.target import numpy_to_cpp_dtype
- return "auto {} = [&](const auto& x){{ return ({}){}; }};".format(name,
- numpy_to_cpp_dtype("float64"),
- visitor(boundary))
-
- raise ValueError("Expression not understood")
-
-
-def name_boundary_lambda(boundary, name):
- define_boundary_lambda(boundary, name + "lambda")
- return name + "lambda"
-
-
-@preamble
-def define_boundary_function(boundary, name):
- gv = name_leafview()
- lambdaname = name_boundary_lambda(boundary, name)
- include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
- if is_stationary():
- return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
- gv,
- lambdaname,
- )
- else:
- params = name_parameters(_driver_data['formdata'])
- return "auto {} = Dune::PDELab::makeInstationaryGridFunctionFromCallable({}, {}, {});".format(name,
- gv,
- lambdaname,
- params)
-
-
-@preamble
-def define_compositegfs_parameterfunction(name, *children):
- return "Dune::PDELab::CompositeGridFunction<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in children),
- name,
- ', '.join(children))
-
-
-def expression_splitter(expr, length):
- if expr is None:
- return (None,) * length
- else:
- from dune.perftool.ufl.execution import Expression, split_expression
- from ufl.classes import ListTensor
- if isinstance(expr, Expression):
- return split_expression(expr)
- elif isinstance(expr, ListTensor):
- return expr.ufl_operands
-
- raise TypeError("No idea how to split {} in components".format(type(expr)))
-
-
-@cached
-def _name_boundary_function(element, boundary, name=None):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(element, (VectorElement, TensorElement)):
- from dune.perftool.generation import get_global_context_value
- basename = get_global_context_value("data").object_names.get(id(boundary), "veczero")
- children = tuple(_name_boundary_function(el, exp, basename + '_' + str(i)) for el, exp, i in zip(element.sub_elements(), expression_splitter(boundary, element.num_sub_elements()), range(element.num_sub_elements())))
- define_compositegfs_parameterfunction(basename, *children)
- return basename
- if isinstance(element, MixedElement):
- children = tuple(name_boundary_function(el) for el in element.sub_elements())
- name = '_'.join(children)
- define_compositegfs_parameterfunction(name, *children)
- return name
-
- from dune.perftool.generation import get_global_context_value
- if name is None:
- name = get_global_context_value("data").object_names.get(id(boundary), "zero")
-
- define_boundary_function(boundary, name)
- return name
-
-
-def name_boundary_function(expr):
- # This is a scalar leaf of the ansatz function tree
- element, (_, boundary) = get_constraints_metadata(expr)
-
- return _name_boundary_function(element, boundary)
-
-
-@preamble
-def interpolate_vector(name, formdata):
- define_vector(name, formdata)
- element = get_trial_element()
- bf = name_boundary_function(element)
- gfs = name_gfs(element)
- return "Dune::PDELab::interpolate({}, {}, {});".format(bf,
- gfs,
- name,
- )
-
-
-def maybe_interpolate_vector(name, formdata):
- element = get_trial_element()
- if has_constraints(element):
- interpolate_vector(name, formdata)
- else:
- define_vector(name, formdata)
-
-
-def name_vector(formdata):
- name = form_name_suffix("x", formdata)
- maybe_interpolate_vector(name, formdata)
- return name
-
-
-@preamble
-def typedef_linearsolver(name):
- include_file("dune/pdelab/backend/istl.hh", filetag="driver")
- return "using {} = Dune::PDELab::ISTLBackend_SEQ_UMFPack;".format(name)
-
-
-def type_linearsolver():
- typedef_linearsolver("LinearSolver")
- return "LinearSolver"
-
-
-@preamble
-def define_linearsolver(name):
- lstype = type_linearsolver()
- return "{} {}(false);".format(lstype, name)
-
-
-def name_linearsolver():
- define_linearsolver("ls")
- return "ls"
-
-
-@preamble
-def define_reduction(name):
- ini = name_initree()
- return "double {} = {}.get<double>(\"reduction\", 1e-12);".format(name, ini)
-
-
-def name_reduction():
- define_reduction("reduction")
- return "reduction"
-
-
-@preamble
-def typedef_stationarylinearproblemsolver(name):
- include_file("dune/pdelab/stationary/linearproblem.hh", filetag="driver")
- gotype = type_gridoperator(_driver_data['formdata'])
- lstype = type_linearsolver()
- xtype = type_vector(_driver_data['formdata'])
- return "using {} = Dune::PDELab::StationaryLinearProblemSolver<{}, {}, {}>;".format(name, gotype, lstype, xtype)
-
-
-def type_stationarylinearproblemsolver():
- typedef_stationarylinearproblemsolver("SLP")
- return "SLP"
-
-
-@preamble
-def define_stationarylinearproblemsolver(name):
- slptype = type_stationarylinearproblemsolver()
- formdata = _driver_data['formdata']
- go = name_gridoperator(formdata)
- ls = name_linearsolver()
- x = name_vector(formdata)
- red = name_reduction()
- return "{} {}({}, {}, {}, {});".format(slptype, name, go, ls, x, red)
-
-
-def name_stationarylinearproblemsolver():
- define_stationarylinearproblemsolver("slp")
- return "slp"
-
-
-@preamble
-def typedef_stationarynonlinearproblemsolver(name, go_type):
- include_file("dune/pdelab/newton/newton.hh", filetag="driver")
- ls_type = type_linearsolver()
- x_type = type_vector(_driver_data['formdata'])
- return "using {} = Dune::PDELab::Newton<{}, {}, {}>;".format(name, go_type, ls_type, x_type)
-
-
-def type_stationarynonlinearproblemssolver(go_type):
- typedef_stationarynonlinearproblemsolver("SNP", go_type)
- return "SNP"
-
-
-@preamble
-def define_stationarynonlinearproblemsolver(name, go_type, go):
- snptype = type_stationarynonlinearproblemssolver(go_type)
- x = name_vector(_driver_data['formdata'])
- ls = name_linearsolver()
- return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
-
-
-def name_stationarynonlinearproblemsolver(go_type, go):
- define_stationarynonlinearproblemsolver("snp", go_type, go)
- return "snp"
-
-
-@preamble
-def typedef_timesteppingmethod(name):
- r_type = type_range()
- explicit = get_option('explicit_time_stepping')
- if explicit:
- return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
- else:
- return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
-
-
-def type_timesteppingmethod():
- typedef_timesteppingmethod("TSM")
- return "TSM"
-
-
-@preamble
-def define_timesteppingmethod(name):
- tsm_type = type_timesteppingmethod()
- explicit = get_option('explicit_time_stepping')
- if explicit:
- return "{} {};".format(tsm_type, name)
- else:
- return "{} {}(1.0);".format(tsm_type, name)
-
-
-def name_timesteppingmethod():
- define_timesteppingmethod("tsm")
- return "tsm"
-
-
-@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'])
- explicit = get_option('explicit_time_stepping')
- if explicit:
- return "using {} = Dune::PDELab::OneStepGridOperator<{},{},false>;".format(name, go_type, mass_go_type)
- else:
- return "using {} = Dune::PDELab::OneStepGridOperator<{},{}>;".format(name, go_type, mass_go_type)
-
-
-def type_instationarygridoperator():
- typedef_instationarygridoperator("IGO")
- return "IGO"
-
-
-@preamble
-def define_instationarygridoperator(name):
- igo_type = type_instationarygridoperator()
- go = name_gridoperator(_driver_data['formdata'])
- mass_go = name_gridoperator(_driver_data['mass_formdata'])
- return "{} {}({}, {});".format(igo_type, name, go, mass_go)
-
-
-def name_instationarygridoperator():
- define_instationarygridoperator("igo")
- return "igo"
-
-
-@preamble
-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'])
- return "using {} = Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}>;".format(name, r_type, igo_type, snp_type, vector_type, vector_type)
-
-
-def type_onestepmethod():
- typedef_onestepmethod("OSM")
- return "OSM"
-
-
-@preamble
-def define_onestepmethod(name):
- ilptype = type_onestepmethod()
- tsm = name_timesteppingmethod()
- igo_type = type_instationarygridoperator()
- igo = name_instationarygridoperator()
- snp = name_stationarynonlinearproblemsolver(igo_type, igo)
- return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)
-
-
-def name_onestepmethod():
- define_onestepmethod("osm")
- return "osm"
-
-
-@preamble
-def typedef_explicitonestepmethod(name):
- r_type = type_range()
- igo_type = type_instationarygridoperator()
- ls_type = type_linearsolver()
- vector_type = type_vector(_driver_data['formdata'])
- return "using {} = Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}>;".format(name, r_type, igo_type, ls_type, vector_type)
-
-
-def type_explicitonestepmethod():
- typedef_explicitonestepmethod("EOSM")
- return "EOSM"
-
-
-@preamble
-def define_explicitonestepmethod(name):
- eosm_type = type_explicitonestepmethod()
- tsm = name_timesteppingmethod()
- igo = name_instationarygridoperator()
- ls = name_linearsolver()
- return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)
-
-
-def name_explicitonestepmethod():
- define_explicitonestepmethod("eosm")
- return "eosm"
-
-
@preamble
def define_mpihelper(name):
include_file("dune/common/parallel/mpihelper.hh", filetag="driver")
@@ -1134,306 +219,6 @@ def name_mpihelper():
return name
-@preamble
-def dune_solve():
- # Test if form is linear in ansatzfunction
- linear = is_linear(_driver_data['formdata'].original_form)
-
- # Test wether we want to do matrix free operator evaluation
- matrix_free = get_option('matrix_free')
-
- # Get right solve command
- if linear and matrix_free:
- formdata = _driver_data['formdata']
- go = name_gridoperator(formdata)
- x = name_vector(formdata)
- include_file("dune/perftool/matrixfree.hh", filetag="driver")
- solve = "solveMatrixFree({},{});".format(go, x)
- elif linear and not matrix_free:
- slp = name_stationarylinearproblemsolver()
- solve = "{}.apply();".format(slp)
- elif not linear and matrix_free:
- # TODO copy of linear case and obviously broken, used to generate something ;)
- formdata = _driver_data['formdata']
- go = name_gridoperator(formdata)
- x = name_vector(formdata)
- include_file("dune/perftool/matrixfree.hh", filetag="driver")
- solve = "solveNonlinearMatrixFree({},{});".format(go, x)
- elif not linear and not matrix_free:
- go_type = type_gridoperator(_driver_data['formdata'])
- go = name_gridoperator(_driver_data['formdata'])
- snp = name_stationarynonlinearproblemsolver(go_type, go)
- solve = "{}.apply();".format(snp)
-
- if get_option('instrumentation_level') >= 2:
- from dune.perftool.pdelab.driver.timings import setup_timer, name_timing_stream
- setup_timer()
- from dune.perftool.generation import post_include
- post_include("HP_DECLARE_TIMER(solve);", filetag="driver")
-
- # Print times after solving
- from dune.perftool.generation import get_global_context_value
- formdatas = get_global_context_value("formdatas")
- print_times = []
- for formdata in formdatas:
- lop_name = name_localoperator(formdata)
- timestream = name_timing_stream()
- print_times.append("{}.dump_timers({}, {}, true);".format(lop_name, timestream, name_timing_identifier()))
-
- solve = ["HP_TIMER_START(solve);",
- "{}".format(solve),
- "HP_TIMER_STOP(solve);",
- "DUMP_TIMER(solve, {}, true);".format(timestream),
- ]
- if get_option('instrumentation_level') >= 3:
- solve.extend(print_times)
-
- return solve
-
-
-@preamble
-def define_vtkfile(name):
- ini = name_initree()
- include_file("string", filetag="driver")
- return "std::string {} = {}.get<std::string>(\"wrapper.vtkcompare.name\", \"output\");".format(name, ini)
-
-
-def name_vtkfile():
- define_vtkfile("vtkfile")
- return "vtkfile"
-
-
-@preamble
-def print_residual():
- ini = name_initree()
- formdata = _driver_data['formdata']
- n_go = name_gridoperator(formdata)
- v = name_vector(formdata)
- t_v = type_vector(formdata)
- include_file("random", system=True, filetag="driver")
-
- return ["if ({}.get<bool>(\"printresidual\", false)) {{".format(ini),
- " using Dune::PDELab::Backend::native;",
- " {} r({});".format(t_v, v),
- " // Setup random input",
- " std::size_t seed = 0;",
- " auto rng = std::mt19937_64(seed);",
- " auto dist = std::uniform_real_distribution<>(-1., 1.);",
- " for (auto& v : {})".format(v),
- " v = dist(rng);",
- " r=0.0;",
- " {}.residual({}, r);".format(n_go, v),
- " Dune::printvector(std::cout, native(r), \"residual vector\", \"row\");",
- "}"]
-
-
-@preamble
-def print_matrix():
- formdata = _driver_data['formdata']
- ini = name_initree()
- t_go = type_gridoperator(formdata)
- n_go = name_gridoperator(formdata)
- v = name_vector(formdata)
- t_v = type_vector(formdata)
-
- return ["if ({}.get<bool>(\"printmatrix\", false)) {{".format(ini),
- " using M = typename {}::Traits::Jacobian;".format(t_go),
- " M m({});".format(n_go),
- " {}.jacobian({},m);".format(n_go, v),
- " using Dune::PDELab::Backend::native;",
- " Dune::printmatrix(std::cout, native(m),\"global stiffness matrix\",\"row\",9,1);",
- "}"]
-
-
-@preamble
-def define_gfs_name(element, prefix=()):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(element, (VectorElement, TensorElement)):
- gfs = name_gfs(element)
- return "{}.name(\"data_{}\");".format(gfs, "_".join(str(p) for p in prefix))
- if isinstance(element, MixedElement):
- for (i, el) in enumerate(element.sub_elements()):
- define_gfs_name(el, prefix + (i,))
- return ""
- # This is a scalar leaf
- gfs = name_gfs(element)
- return "{}.name(\"data_{}\");".format(gfs, "_".join(str(p) for p in prefix))
-
-
-def type_predicate():
- include_file("dune/perftool/vtkpredicate.hh", filetag="driver")
- return "CuttingPredicate"
-
-
-@preamble
-def define_predicate(name):
- t = type_predicate()
- return "{} {};".format(t, name)
-
-
-def name_predicate():
- define_predicate("predicate")
- return "predicate"
-
-
-@preamble
-def vtkoutput():
- element = get_trial_element()
- define_gfs_name(element)
- include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
- vtkwriter = name_vtkwriter()
- gfs = name_gfs(element)
- vec = name_vector(_driver_data['formdata'])
- vtkfile = name_vtkfile()
- predicate = name_predicate()
- dune_solve()
- 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()
-
- return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate),
- "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]
-
-
-@preamble
-def define_time(name):
- return "double {} = 0.0;".format(name)
-
-
-def name_time():
- define_time("time")
- return "time"
-
-
-@preamble
-def typedef_vtk_sequence_writer(name):
- include_file("dune/grid/io/file/vtk/vtksequencewriter.hh", filetag="driver")
- gv_type = type_leafview()
- return "using {} = Dune::VTKSequenceWriter<{}>;".format(name, gv_type)
-
-
-def type_vtk_sequence_writer():
- typedef_vtk_sequence_writer("VTKSW")
- return "VTKSW"
-
-
-@preamble
-def define_vtk_sequence_writer(name):
- vtksw_type = type_vtk_sequence_writer()
- vtkw_type = type_vtkwriter()
- vtkw = name_vtkwriter()
- vtkfile = name_vtkfile()
- return "{} {}(std::make_shared<{}>({}), {});".format(vtksw_type, name, vtkw_type, vtkw, vtkfile)
-
-
-def name_vtk_sequence_writer():
- define_vtk_sequence_writer("vtkSequenceWriter")
- return "vtkSequenceWriter"
-
-
-@preamble
-def visualize_initial_condition():
- include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
- vtkwriter = name_vtk_sequence_writer()
- element = get_trial_element()
- define_gfs_name(element)
- gfs = name_gfs(element)
- vector = name_vector(_driver_data['formdata'])
- predicate = name_predicate()
- 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 time_loop():
- ini = name_initree()
- formdata = _driver_data['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)
- vector_type = type_vector(formdata)
- vector = name_vector(formdata)
-
- # 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:
- boundary = name_boundary_function(expr)
- osm = name_onestepmethod()
- apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
-
- # Setup visualization
- visualize_initial_condition()
- vtk_sequence_writer = name_vtk_sequence_writer()
-
- return ["",
- "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
- "double dt = {}.get<double>(\"instat.dt\", 0.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),
- "",
- " // Do time step",
- " {} {}new({});".format(vector_type, vector, vector),
- " {}".format(apply_call),
- "",
- " // Accept new time step",
- " {} = {}new;".format(vector, vector),
- " time += dt;",
- "",
- " // Output to VTK File",
- " {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
- "}",
- ""]
-
-
-def solve_instationary():
- # Test if form is linear in ansatzfunction
- linear = is_linear(_driver_data['formdata'].original_form)
-
- # Test wether we want to do matrix free operator evaluation
- matrix_free = get_option('matrix_free')
-
- # Create time loop
- if linear and matrix_free:
- assert False
- elif linear and not matrix_free:
- time_loop()
- if not linear and matrix_free:
- assert False
- elif not linear and not matrix_free:
- assert False
-
- 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
-def return_statement():
- from dune.perftool.pdelab.driver.error import name_test_fail_variable
- fail = name_test_fail_variable()
- return "return {};".format(fail)
-
-
def _blockstructured_degree_helper(element, op):
element._degree = op(element._degree, get_option("number_of_blocks"))
from ufl import MixedElement, VectorElement, TensorElement
@@ -1474,9 +259,12 @@ def generate_driver(formdatas, data):
evaluate_residual_timer()
apply_jacobian_timer()
elif is_stationary():
- # We could also use solve if we are not interested in visualization
+ from dune.perftool.pdelab.driver.solve import dune_solve
+ vec = dune_solve()
+ from dune.perftool.pdelab.driver.vtk import vtkoutput
vtkoutput()
else:
+ from dune.perftool.pdelab.driver.instationary import solve_instationary
solve_instationary()
# Make sure that timestream is declared before retrieving chache items
@@ -1485,6 +273,7 @@ def generate_driver(formdatas, data):
setup_timer()
timestream = name_timing_stream()
+ from dune.perftool.pdelab.driver.error import return_statement
return_statement()
from dune.perftool.generation import retrieve_cache_items
diff --git a/python/dune/perftool/pdelab/driver/constraints.py b/python/dune/perftool/pdelab/driver/constraints.py
new file mode 100644
index 00000000..a75016c2
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/constraints.py
@@ -0,0 +1,144 @@
+from dune.perftool.generation import (include_file,
+ preamble,
+ )
+from dune.perftool.pdelab.driver import (FEM_name_mangling,
+ get_formdata,
+ get_trial_element,
+ )
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
+ name_trial_gfs,
+ type_range,
+ type_trial_gfs,
+ preprocess_leaf_data,
+ )
+
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+
+
+def name_assembled_constraints():
+ name = name_constraintscontainer()
+ define_constraintscontainer(name)
+ assemble_constraints(name)
+ return name
+
+
+@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,
+ name,
+ )
+
+
+def name_bctype_function(element, is_dirichlet):
+ if isinstance(element, (VectorElement, TensorElement)):
+ subel = element.sub_element()[0]
+ subgfs = name_bctype_function(subel, is_dirichlet[:subel.value_size()])
+ name = "{}_pow{}bctype".format(subgfs, element.num_sub_elements())
+ define_power_bctype_function(element, 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()]))
+ k = k + subel.value_size()
+ name = "_".join(subgfs)
+ define_composite_bctype_function(element, is_dirichlet, name, subgfs)
+ return name
+ else:
+ assert isinstance(element, FiniteElement)
+ name = "{}_bctype".format(FEM_name_mangling(element).lower())
+ define_bctype_function(element, is_dirichlet[0], name)
+ return name
+
+
+@preamble
+def define_bctype_function(element, is_dirichlet, name):
+ gv = name_leafview()
+ bctype_lambda = name_bctype_lambda(name, is_dirichlet)
+ include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
+ return "auto {} = Dune::PDELab::makeBoundaryConditionFromCallable({}, {});".format(name,
+ gv,
+ bctype_lambda,
+ )
+
+
+@preamble
+def define_power_bctype_function(element, name, subgfs):
+ include_file('dune/pdelab/constraints/common/constraintsparameters.hh', filetag='driver')
+ return "Dune::PDELab::PowerConstraintsParameters<decltype({}), {}> {}({});".format(subgfs, element.num_sub_elements(), name, subgfs)
+
+
+@preamble
+def define_composite_bctype_function(element, is_dirichlet, name, subgfs):
+ include_file('dune/pdelab/constraints/common/constraintsparameters.hh', filetag='driver')
+ return "Dune::PDELab::CompositeConstraintsParameters<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in subgfs),
+ name,
+ ', '.join(c for c in subgfs)
+ )
+
+
+def name_bctype_lambda(name, func):
+ name = name + "_lambda"
+ define_intersection_lambda(name, func)
+ return name
+
+
+@preamble
+def define_intersection_lambda(name, func):
+ from dune.perftool.ufl.execution import Expression
+ from ufl.classes import Expr
+ if func is None:
+ func = 0.
+ if isinstance(func, (int, float)):
+ return "auto {} = [&](const auto& x){{ return {}; }};".format(name, float(func))
+ elif isinstance(func, Expression):
+ if expression.is_global:
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, func.c_expr[0])
+ else:
+ return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, func.c_expr[0])
+ elif isinstance(func, Expr):
+ # Set up a visitor
+ with global_context(integral_type="exterior_facet", formdata=get_formdata(), driver=True):
+ from dune.perftool.ufl.visitor import UFL2LoopyVisitor
+ from dune.perftool.pdelab import PDELabInterface
+ visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
+ from pymbolic.mapper.c_code import CCodeMapper
+ ccm = CCodeMapper()
+ expr = visitor(func)
+ return "auto {} = [&](const auto& x){{ return {}; }};".format(name, ccm(expr))
+
+ raise ValueError("Expression not understood")
+
+
+@preamble
+def typedef_constraintscontainer(name):
+ gfs = type_trial_gfs()
+ r = type_range()
+ return "using {} = {}::ConstraintsContainer<{}>::Type;".format(name, gfs, r)
+
+
+def type_constraintscontainer():
+ name = "{}_CC".format(type_trial_gfs())
+ typedef_constraintscontainer(name)
+ return name
+
+
+@preamble
+def define_constraintscontainer(name):
+ cctype = type_constraintscontainer()
+ return ["{} {};".format(cctype, name), "{}.clear();".format(name)]
+
+
+def name_constraintscontainer():
+ gfs = name_trial_gfs()
+ name = "{}_cc".format(gfs)
+ element = get_trial_element()
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ define_constraintscontainer(name)
+ return name
diff --git a/python/dune/perftool/pdelab/driver/error.py b/python/dune/perftool/pdelab/driver/error.py
index c02147a9..dbbb0b9c 100644
--- a/python/dune/perftool/pdelab/driver/error.py
+++ b/python/dune/perftool/pdelab/driver/error.py
@@ -5,6 +5,17 @@ from dune.perftool.generation import (cached,
preamble,
)
from dune.perftool.options import get_option
+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.interpolate import (interpolate_vector,
+ name_boundary_function,
+ )
+from dune.perftool.pdelab.driver.solve import (define_vector,
+ name_vector,
+ )
@preamble
@@ -18,37 +29,26 @@ def name_test_fail_variable():
return name
-def name_vector_from_solution_expression():
- interpolate_solution_expression("solution")
- return "solution"
+def name_exact_solution():
+ name = "solution"
+ define_vector(name, get_formdata())
+ interpolate_exact_solution(name)
+ return name
-@preamble
-def interpolate_solution_expression(name):
- from dune.perftool.pdelab.driver import (define_vector,
- get_formdata,
- get_trial_element,
- name_gfs,
- )
- formdata = get_formdata()
- define_vector(name, formdata)
+def interpolate_exact_solution(name):
element = get_trial_element()
- sol = name_solution_function()
- gfs = name_gfs(element)
- return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
- gfs,
- name,
- )
+ 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():
- from dune.perftool.pdelab.driver import (get_formdata,
- name_vector,
- )
-
v = name_vector(get_formdata())
- solution = name_vector_from_solution_expression()
+ solution = name_exact_solution()
fail = name_test_fail_variable()
return ["",
@@ -64,7 +64,6 @@ def compare_dofs():
" {} = true;".format(fail)]
-
def type_discrete_grid_function(gfs):
return "DGF_{}".format(gfs.upper())
@@ -77,149 +76,54 @@ def define_discrete_grid_function(gfs, vector_name, dgf_name):
def name_discrete_grid_function(gfs, vector_name):
- dgf_name = gfs + "_dgf"
+ dgf_name = "{}_dgf".format(gfs)
define_discrete_grid_function(gfs, vector_name, dgf_name)
return dgf_name
-def type_subgfs(element, tree_path):
- from dune.perftool.pdelab.driver import type_gfs
- include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
- gfs = type_gfs(element)
- return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(tree_path))
-
-
@preamble
-def define_subgfs(name, element, tree_path):
- t = type_subgfs(element, tree_path)
-
- from dune.perftool.pdelab.driver import name_gfs
- gfs = name_gfs(element)
-
- return '{} {}({});'.format(t, name, gfs)
-
-
-def name_subgfs(element, tree_path):
- from dune.perftool.pdelab.driver import name_gfs
- gfs = name_gfs(element)
- if len(tree_path) == 0:
- return gfs
- else:
- name = "subgfs_{}_{}".format(gfs, '_'.join(tree_path))
- define_subgfs(name, element, tree_path)
- return name
-
-
-@cached
-def name_solution_function(tree_path=()):
- from dune.perftool.generation import get_global_context_value
- expr = get_global_context_value("data").object_by_name[get_option("exact_solution_expression")]
-
- from dune.perftool.ufl.execution import Expression, split_expression
- from ufl.classes import ListTensor
- for i in tree_path:
- if isinstance(expr, Expression):
- expr = split_expression(expr)[int(i)]
- elif isinstance(expr, tuple):
- expr = expr[int(i)]
- elif isinstance(expr, ListTensor):
- expr = expr.ufl_operands[int(i)]
- else:
- raise TypeError("No idea how to split {}".format(type(expr)))
-
- from dune.perftool.generation import get_global_context_value
- try:
- name = get_global_context_value("data").object_names[id(expr)]
- except KeyError:
- from dune.perftool.generation import get_counter
- name = 'generic_{}'.format(get_counter('__generic'))
-
- from dune.perftool.pdelab.driver import define_boundary_function
- define_boundary_function(expr, name)
-
- return name
-
-
-@preamble
-def typedef_difference_squared_adapter(name, tree_path):
- from dune.perftool.pdelab.driver import (get_formdata,
- get_trial_element,
- name_vector,
- )
- element = get_trial_element()
- formdata = get_formdata()
-
- solution_function = name_solution_function(tree_path)
- gfs = name_subgfs(element, tree_path)
- vector = name_vector(formdata)
+def typedef_difference_squared_adapter(name):
+ sol = name_exact_solution()
+ vector = name_vector(get_formdata())
dgf = name_discrete_grid_function(gfs, vector)
- return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, solution_function, dgf)
+ return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, sol, dgf)
-def type_difference_squared_adapter(tree_path):
- t = 'DifferenceSquared_{}'.format('_'.join(tree_path))
- typedef_difference_squared_adapter(t, tree_path)
- return t
+def type_difference_squared_adapter():
+ name = 'DifferenceSquaredAdapter'
+ typedef_difference_squared_adapter(name)
+ return name
@preamble
-def define_difference_squared_adapter(name, tree_path):
- from dune.perftool.pdelab.driver import (get_formdata,
- get_trial_element,
- name_vector,
- )
- element = get_trial_element()
- formdata = get_formdata()
-
- t = type_difference_squared_adapter(tree_path)
- solution_function = name_solution_function(tree_path)
- gfs = name_subgfs(element, tree_path)
- vector = name_vector(formdata)
+def define_difference_squared_adapter(name):
+ t = type_difference_squared_adapter()
+ sol = name_exact_solution()
+ vector = name_vector(get_formdata())
dgf = name_discrete_grid_function(gfs, vector)
- return '{} {}({}, {});'.format(t, name, solution_function, dgf)
+ return '{} {}({}, {});'.format(t, name, sol, dgf)
-def name_difference_squared_adapter(tree_path):
- name = 'differencesquared_{}'.format('_'.join(tree_path))
- define_difference_squared_adapter(name, tree_path)
+def name_difference_squared_adapter():
+ name = 'differencesquared_adapter'
+ define_difference_squared_adapter(name)
return name
@preamble
-def _accumulate_L2_squared(tree_path):
- dsa = name_difference_squared_adapter(tree_path)
- t_dsa = type_difference_squared_adapter(tree_path)
+def accumulate_L2_squared():
+ dsa = name_difference_squared_adapter()
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",
- " // only subgfs with treepath: {}".format(', '.join(tree_path)),
- " typename {}::Traits::RangeType err(0.0);".format(t_dsa),
- " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
- " {} += err;".format(accum_error),
- " std::cout << \"Accumlated L2 Error for tree path {}: \" << err << std::endl;".format(tree_path),
- "}"]
-
-
-def accumulate_L2_squared(tree_path=()):
- from dune.perftool.pdelab.driver import get_trial_element
- element = get_trial_element()
- from ufl.functionview import select_subelement
- from ufl.classes import MultiIndex, FixedIndex
- element = select_subelement(element, MultiIndex(tuple(FixedIndex(int(i)) for i in tree_path)))
-
- from ufl import MixedElement
- if isinstance(element, MixedElement):
- for i, subel in enumerate(element.sub_elements()):
- accumulate_L2_squared(tree_path + (str(i),))
- else:
- _accumulate_L2_squared(tree_path)
+ return ["// L2 error squared of difference between numerical",
+ "// solution and the interpolation of exact solution",
+ " Dune::PDELab::integrateGridFunction({}, {}, 10);".format(dsa, accum_error),
+ ]
@preamble
@@ -242,3 +146,10 @@ def compare_L2_squared():
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")),
" {} = true;".format(fail)]
+
+
+@preamble
+def return_statement():
+ from dune.perftool.pdelab.driver.error import name_test_fail_variable
+ fail = name_test_fail_variable()
+ return "return {};".format(fail)
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
new file mode 100644
index 00000000..b439ded1
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
@@ -0,0 +1,345 @@
+from dune.perftool.generation import (include_file,
+ preamble,
+ )
+from dune.perftool.options import get_option, set_option
+from dune.perftool.pdelab.driver import (FEM_name_mangling,
+ get_cell,
+ get_dimension,
+ get_test_element,
+ get_trial_element,
+ isDG,
+ isPk,
+ isQk,
+ isQuadrilateral,
+ isSimplical,
+ name_initree,
+ preprocess_leaf_data,
+ )
+
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+
+
+@preamble
+def typedef_domainfield(name):
+ gridt = type_grid()
+ return "using {} = {}::ctype;".format(name, gridt)
+
+
+def type_domainfield():
+ typedef_domainfield("DF")
+ return "DF"
+
+
+@preamble
+def typedef_range(name):
+ if get_option('opcounter'):
+ from dune.perftool.pdelab.driver.timings import setup_timer
+ setup_timer()
+ return "using {} = oc::OpCounter<double>;".format(name)
+ else:
+ return "using {} = double;".format(name)
+
+
+def type_range():
+ name = "RangeType"
+ typedef_range(name)
+ return name
+
+
+@preamble
+def typedef_grid(name):
+ dim = get_dimension()
+ if isQuadrilateral(get_trial_element().cell()):
+ # For Yasp Grids the jacobi of the transformation is diagonal and constant on each cell
+ set_option('diagonal_transformation_matrix', True)
+ set_option('constant_transformation_matrix', True)
+
+ range_type = type_range()
+ gridt = "Dune::YaspGrid<{0}, Dune::EquidistantCoordinates<{1}, {0}>>".format(dim, range_type)
+ include_file("dune/grid/yaspgrid.hh", filetag="driver")
+ else:
+ if isSimplical(get_trial_element().cell()):
+ # gridt = "Dune::UGGrid<{}>".format(dim)
+ # include_file("dune/grid/uggrid.hh", filetag="driver")
+ gridt = "Dune::ALUGrid<{}, {}, Dune::simplex, Dune::conforming>".format(dim, dim)
+ include_file("dune/alugrid/grid.hh", filetag="driver")
+ else:
+ raise PerftoolCodegenError("Cant match your geometry with a DUNE grid. Please report bug.")
+ return "using {} = {};".format(name, gridt)
+
+
+def type_grid():
+ name = "Grid"
+ typedef_grid(name)
+ return name
+
+
+@preamble
+def define_grid(name):
+ include_file("dune/testtools/gridconstruction.hh", filetag="driver")
+ ini = name_initree()
+ _type = type_grid()
+ return ["IniGridFactory<{}> factory({});".format(_type, ini),
+ "std::shared_ptr<{}> grid = factory.getGrid();".format(_type)]
+
+
+def name_grid():
+ name = "grid"
+ define_grid(name)
+ return name
+
+
+@preamble
+def typedef_leafview(name):
+ grid = type_grid()
+ return "using {} = {}::LeafGridView;".format(name, grid)
+
+
+def type_leafview():
+ name = "GV"
+ typedef_leafview(name)
+ return name
+
+
+@preamble
+def define_leafview(name):
+ _type = type_leafview()
+ grid = name_grid()
+ return "{} {} = {}->leafGridView();".format(_type, name, grid)
+
+
+def name_leafview():
+ name = "gv"
+ define_leafview(name)
+ return name
+
+
+@preamble
+def typedef_fem(element, name):
+ gv = type_leafview()
+ df = type_domainfield()
+ r = type_range()
+ dim = get_dimension()
+ if isQk(element):
+ include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
+ return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
+ if isPk(element):
+ include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
+ return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
+ if isDG(element):
+ if isQuadrilateral(element.cell()):
+ include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
+ # TODO allow switching the basis here!
+ return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, df, r, element.degree(), dim)
+ if isSimplical(element.cell()):
+ include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
+ return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;".format(name, df, r, element.degree(), dim)
+ raise NotImplementedError("Geometry type not known in code generation")
+
+ raise NotImplementedError("FEM not implemented in dune-perftool")
+
+
+def type_fem(element):
+ name = "{}_FEM".format(FEM_name_mangling(element).upper())
+ typedef_fem(element, name)
+ return name
+
+
+@preamble
+def define_fem(element, name):
+ femtype = type_fem(element)
+ from dune.perftool.pdelab.driver import isDG
+ if isDG(element):
+ return "{} {};".format(femtype, name)
+ else:
+ gv = name_leafview()
+ return "{} {}({});".format(femtype, name, gv)
+
+
+def name_fem(element):
+ assert isinstance(element, FiniteElement)
+ name = "{}_fem".format(FEM_name_mangling(element).lower())
+ define_fem(element, name)
+ return name
+
+
+def name_trial_gfs():
+ element = get_trial_element()
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ return name_gfs(element, is_dirichlet)
+
+
+def name_test_gfs():
+ element = get_test_element()
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ return name_gfs(element, is_dirichlet)
+
+
+def name_gfs(element, is_dirichlet):
+ 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())
+ 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()]))
+ k = k + subel.value_size()
+ name = "_".join(subgfs)
+ define_composite_gfs(element, is_dirichlet, name, subgfs)
+ return name
+ else:
+ assert isinstance(element, FiniteElement)
+ name = "{}{}_gfs".format(FEM_name_mangling(element).lower(),
+ "_dirichlet" if is_dirichlet else "",
+ )
+ define_gfs(element, is_dirichlet, name)
+ return name
+
+
+def type_test_gfs():
+ element = get_test_element()
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ return type_gfs(element, is_dirichlet)
+
+
+def type_trial_gfs():
+ element = get_trial_element()
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ return type_gfs(element, is_dirichlet)
+
+
+def type_gfs(element, is_dirichlet):
+ if isinstance(element, (VectorElement, TensorElement)):
+ subel = element.sub_element()[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)
+ return name
+ if isinstance(element, MixedElement):
+ k = 0
+ subgfs = []
+ for subel in element.sub_elements():
+ 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)
+ return name
+ else:
+ assert isinstance(element, FiniteElement)
+ name = "{}{}_GFS".format(FEM_name_mangling(element).upper(),
+ "_dirichlet" if is_dirichlet else "",
+ )
+ typedef_gfs(element, is_dirichlet, name)
+ return name
+
+
+@preamble
+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)
+
+
+@preamble
+def define_power_gfs(element, is_dirichlet, name, subgfs):
+ gv = name_leafview()
+ fem = name_fem(element.sub_elements()[0])
+ return "{} {}({}, {});".format(gfstype, name, gv, fem)
+
+
+@preamble
+def define_composite_gfs(element, is_dirichlet, name, subgfs):
+ gfstype = type_gfs(element, is_dirichlet)
+ return "{} {}({});".format(gfstype, name, ", ".join(subgfs))
+
+
+@preamble
+def typedef_gfs(element, is_dirichlet, name):
+ vb = type_vectorbackend(element)
+ gv = type_leafview()
+ fem = type_fem(element)
+ cass = type_constraintsassembler(bool(is_dirichlet))
+ return "using {} = Dune::PDELab::GridFunctionSpace<{}, {}, {}, {}>;".format(name, gv, fem, cass, vb)
+
+
+@preamble
+def typedef_power_gfs(element, is_dirichlet, name, subgfs):
+ include_file("dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh", filetag="driver")
+ gv = type_leafview()
+ vb = type_vectorbackend(element)
+ fem = type_fem(element.sub_elements()[0])
+ dim = element.num_sub_elements()
+ cass = type_constraintsassembler(bool(is_dirichlet[0]))
+ return "using {} = Dune::PDELab::VectorGridFunctionSpace<{}, {}, {}, {}, {}, {}>;".format(name, gv, fem, dim, vb, vb, cass)
+
+
+@preamble
+def typedef_composite_gfs(element, name, subgfs):
+ vb = type_vectorbackend(element)
+ ot = type_orderingtag()
+ args = ", ".join(subgfs)
+ return "using {} = Dune::PDELab::CompositeGridFunctionSpace<{}, {}, {}>;".format(name, vb, ot, args)
+
+
+@preamble
+def define_blocksize(name, element):
+ from dune.perftool.pdelab.driver import isDG, isQuadrilateral
+ assert isDG(element)
+ assert isQuadrilateral(element.cell())
+ dimension = get_dimension()
+ degree = element.degree()
+ return "static const int {} = Dune::QkStuff::QkSize<{}, {}>::value;".format(name, degree, dimension)
+
+
+def name_blocksize(element):
+ name = "blocksize"
+ define_blocksize(name, element)
+ return name
+
+
+@preamble
+def typedef_vectorbackend(name, element):
+ include_file("dune/pdelab/backend/istl.hh", filetag="driver")
+ if get_option("fastdg"):
+ blocksize = name_blocksize(element)
+ return "using {} = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed, {}>;".format(name, blocksize)
+ else:
+ return "using {} = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none, 1>;".format(name)
+
+
+def type_vectorbackend(element):
+ name = "VectorBackend"
+ typedef_vectorbackend(name, element)
+ return name
+
+
+def type_orderingtag():
+ return "Dune::PDELab::LexicographicOrderingTag"
+
+
+@preamble
+def typedef_dirichlet_constraintsassembler(name):
+ include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
+ return "using {} = Dune::PDELab::ConformingDirichletConstraints;".format(name)
+
+
+@preamble
+def typedef_no_constraintsassembler(name):
+ return "using {} = Dune::PDELab::NoConstraints;".format(name)
+
+
+def type_constraintsassembler(is_dirichlet):
+ assert isinstance(is_dirichlet, bool)
+ if is_dirichlet:
+ name = "DirichletConstraintsAssember"
+ typedef_dirichlet_constraintsassembler(name)
+ return name
+ else:
+ name = "NoConstraintsAssembler"
+ typedef_no_constraintsassembler(name)
+ return name
diff --git a/python/dune/perftool/pdelab/driver/gridoperator.py b/python/dune/perftool/pdelab/driver/gridoperator.py
new file mode 100644
index 00000000..500d8bc7
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/gridoperator.py
@@ -0,0 +1,173 @@
+from dune.perftool.generation import (get_global_context_value,
+ include_file,
+ preamble,
+ )
+from dune.perftool.pdelab.driver import (form_name_suffix,
+ get_cell,
+ get_dimension,
+ get_test_element,
+ get_trial_element,
+ isQuadrilateral,
+ name_initree,
+ )
+from dune.perftool.pdelab.driver.constraints import (name_assembled_constraints,
+ type_constraintscontainer,
+ )
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_test_gfs,
+ name_trial_gfs,
+ preprocess_leaf_data,
+ type_domainfield,
+ type_range,
+ type_test_gfs,
+ type_trial_gfs,
+ )
+from dune.perftool.pdelab.localoperator import localoperator_basename
+from dune.perftool.pdelab.parameter import parameterclass_basename
+from dune.perftool.options import get_option
+
+
+@preamble
+def typedef_gridoperator(name, formdata):
+ ugfs = type_trial_gfs()
+ vgfs = type_test_gfs()
+ lop = type_localoperator(formdata)
+ cc = type_constraintscontainer()
+ mb = type_matrixbackend()
+ df = type_domainfield()
+ r = type_range()
+ if get_option("fastdg"):
+ if not get_option("sumfact"):
+ raise PerftoolCodegenError("FastDGGridOperator is only implemented for sumfactorization.")
+ include_file("dune/pdelab/gridoperator/fastdg.hh", filetag="driver")
+ return "using {} = Dune::PDELab::FastDGGridOperator<{}, {}, {}, {}, {}, {}, {}, {}, {}>;".format(name, ugfs, vgfs, lop, mb, df, r, r, cc, cc)
+ else:
+ include_file("dune/pdelab/gridoperator/gridoperator.hh", filetag="driver")
+ return "using {} = Dune::PDELab::GridOperator<{}, {}, {}, {}, {}, {}, {}, {}, {}>;".format(name, ugfs, vgfs, lop, mb, df, r, r, cc, cc)
+
+
+def type_gridoperator(formdata):
+ name = form_name_suffix("GO", formdata).upper()
+ typedef_gridoperator(name, formdata)
+ return name
+
+
+@preamble
+def define_gridoperator(name, formdata):
+ gotype = type_gridoperator(formdata)
+ ugfs = name_trial_gfs()
+ vgfs = name_test_gfs()
+ if ugfs != vgfs:
+ raise NotImplementedError("Non-Galerkin methods currently not supported!")
+ cc = name_assembled_constraints()
+ lop = name_localoperator(formdata)
+ mb = name_matrixbackend()
+ return ["{} {}({}, {}, {}, {}, {}, {});".format(gotype, name, ugfs, cc, vgfs, cc, lop, mb),
+ "std::cout << \"gfs with \" << {}.size() << \" dofs generated \"<< std::endl;".format(ugfs),
+ "std::cout << \"cc with \" << {}.size() << \" dofs generated \"<< std::endl;".format(cc)]
+
+
+def name_gridoperator(formdata):
+ name = form_name_suffix("go", formdata)
+ define_gridoperator(name, formdata)
+ return name
+
+
+@preamble
+def typedef_localoperator(name, formdata):
+ ugfs = type_trial_gfs()
+ vgfs = type_test_gfs()
+ data = get_global_context_value("data")
+ filename = get_option("operator_file")
+ include_file(filename, filetag="driver")
+ lopname = localoperator_basename(formdata, data)
+ range_type = type_range()
+ return "using {} = {}<{}, {}, {}>;".format(name, lopname, ugfs, vgfs, range_type)
+
+
+def type_localoperator(formdata):
+ name = form_name_suffix("LOP", formdata).upper()
+ typedef_localoperator(name, formdata)
+ return name
+
+
+@preamble
+def define_localoperator(name, formdata):
+ loptype = type_localoperator(formdata)
+ ini = name_initree()
+ params = name_parameters(formdata)
+ return "{} {}({}, {});".format(loptype, name, ini, params)
+
+
+def name_localoperator(formdata):
+ name = form_name_suffix("lop", formdata)
+ define_localoperator(name, formdata)
+ return name
+
+
+@preamble
+def define_dofestimate(name):
+ # Provide a worstcase estimate for the number of entries per row based
+ # on the given gridfunction space and cell geometry
+ if isQuadrilateral(get_cell()):
+ geo_factor = 2**get_dimension()
+ else:
+ if get_dimension() < 3:
+ geo_factor = 3*get_dimension()
+ else:
+ # TODO no idea how a generic estimate for 3D simplices looks like
+ geo_factor = 12
+
+ gfs = name_trial_gfs()
+ ini = name_initree()
+
+ return ["int generic_dof_estimate = {} * {}.maxLocalSize();".format(geo_factor, gfs),
+ "int {} = {}.get<int>(\"istl.number_of_nnz\", generic_dof_estimate);".format(name, ini)]
+
+
+def name_dofestimate():
+ name = "dofestimate"
+ define_dofestimate(name)
+ return name
+
+
+@preamble
+def typedef_matrixbackend(name):
+ include_file("dune/pdelab/backend/istl.hh", filetag="driver")
+ return "using {} = Dune::PDELab::istl::BCRSMatrixBackend<>;".format(name)
+
+
+def type_matrixbackend():
+ name = "MatrixBackend"
+ typedef_matrixbackend(name)
+ return name
+
+
+@preamble
+def define_matrixbackend(name):
+ mbtype = type_matrixbackend()
+ dof = name_dofestimate()
+ return "{} {}({});".format(mbtype, name, dof)
+
+
+def name_matrixbackend():
+ name = "mb"
+ define_matrixbackend(name)
+ return name
+
+
+def type_parameters(formdata):
+ data = get_global_context_value("data")
+ name = parameterclass_basename(formdata, data)
+ return name
+
+
+@preamble
+def define_parameters(name, formdata):
+ partype = type_parameters(formdata)
+ return "{} {};".format(partype, name)
+
+
+def name_parameters(formdata):
+ name = form_name_suffix("params", formdata)
+ define_parameters(name, formdata)
+ return name
diff --git a/python/dune/perftool/pdelab/driver/instationary.py b/python/dune/perftool/pdelab/driver/instationary.py
new file mode 100644
index 00000000..bf8759e7
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/instationary.py
@@ -0,0 +1,213 @@
+from dune.perftool.generation import preamble
+from dune.perftool.pdelab.driver import (get_formdata,
+ is_linear,
+ )
+from dune.perftool.pdelab.driver.gridoperator import (type_gridoperator,)
+from dune.perftool.pdelab.driver.solve import (print_matrix,
+ print_residual,
+ )
+from dune.perftool.options import get_option
+
+
+def solve_instationary():
+ # Test if form is linear in ansatzfunction
+ linear = is_linear()
+
+ # Test wether we want to do matrix free operator evaluation
+ matrix_free = get_option('matrix_free')
+
+ # Create time loop
+ if linear and matrix_free:
+ assert False
+ elif linear and not matrix_free:
+ time_loop()
+ if not linear and matrix_free:
+ assert False
+ elif not linear and not matrix_free:
+ assert False
+
+ 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
+def time_loop():
+ ini = name_initree()
+ 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)
+ vector_type = type_vector(formdata)
+ vector = name_vector(formdata)
+
+ # 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:
+ boundary = name_boundary_function(expr)
+ osm = name_onestepmethod()
+ apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
+
+ # Setup visualization
+ visualize_initial_condition()
+ vtk_sequence_writer = name_vtk_sequence_writer()
+
+ return ["",
+ "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
+ "double dt = {}.get<double>(\"instat.dt\", 0.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),
+ "",
+ " // Do time step",
+ " {} {}new({});".format(vector_type, vector, vector),
+ " {}".format(apply_call),
+ "",
+ " // Accept new time step",
+ " {} = {}new;".format(vector, vector),
+ " time += dt;",
+ "",
+ " // Output to VTK File",
+ " {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
+ "}",
+ ""]
+
+
+@preamble
+def define_time(name):
+ return "double {} = 0.0;".format(name)
+
+
+def name_time():
+ define_time("time")
+ return "time"
+
+
+
+@preamble
+def typedef_timesteppingmethod(name):
+ r_type = type_range()
+ explicit = get_option('explicit_time_stepping')
+ if explicit:
+ return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
+ else:
+ return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
+
+
+def type_timesteppingmethod():
+ typedef_timesteppingmethod("TSM")
+ return "TSM"
+
+
+@preamble
+def define_timesteppingmethod(name):
+ tsm_type = type_timesteppingmethod()
+ explicit = get_option('explicit_time_stepping')
+ if explicit:
+ return "{} {};".format(tsm_type, name)
+ else:
+ return "{} {}(1.0);".format(tsm_type, name)
+
+
+def name_timesteppingmethod():
+ define_timesteppingmethod("tsm")
+ return "tsm"
+
+
+@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'])
+ explicit = get_option('explicit_time_stepping')
+ if explicit:
+ return "using {} = Dune::PDELab::OneStepGridOperator<{},{},false>;".format(name, go_type, mass_go_type)
+ else:
+ return "using {} = Dune::PDELab::OneStepGridOperator<{},{}>;".format(name, go_type, mass_go_type)
+
+
+def type_instationarygridoperator():
+ typedef_instationarygridoperator("IGO")
+ return "IGO"
+
+
+@preamble
+def define_instationarygridoperator(name):
+ igo_type = type_instationarygridoperator()
+ go = name_gridoperator(_driver_data['formdata'])
+ mass_go = name_gridoperator(_driver_data['mass_formdata'])
+ return "{} {}({}, {});".format(igo_type, name, go, mass_go)
+
+
+def name_instationarygridoperator():
+ define_instationarygridoperator("igo")
+ return "igo"
+
+
+@preamble
+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'])
+ return "using {} = Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}>;".format(name, r_type, igo_type, snp_type, vector_type, vector_type)
+
+
+def type_onestepmethod():
+ typedef_onestepmethod("OSM")
+ return "OSM"
+
+
+@preamble
+def define_onestepmethod(name):
+ ilptype = type_onestepmethod()
+ tsm = name_timesteppingmethod()
+ igo_type = type_instationarygridoperator()
+ igo = name_instationarygridoperator()
+ snp = name_stationarynonlinearproblemsolver(igo_type, igo)
+ return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)
+
+
+def name_onestepmethod():
+ define_onestepmethod("osm")
+ return "osm"
+
+
+@preamble
+def typedef_explicitonestepmethod(name):
+ r_type = type_range()
+ igo_type = type_instationarygridoperator()
+ ls_type = type_linearsolver()
+ vector_type = type_vector(_driver_data['formdata'])
+ return "using {} = Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}>;".format(name, r_type, igo_type, ls_type, vector_type)
+
+
+def type_explicitonestepmethod():
+ typedef_explicitonestepmethod("EOSM")
+ return "EOSM"
+
+
+@preamble
+def define_explicitonestepmethod(name):
+ eosm_type = type_explicitonestepmethod()
+ tsm = name_timesteppingmethod()
+ igo = name_instationarygridoperator()
+ ls = name_linearsolver()
+ return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)
+
+
+def name_explicitonestepmethod():
+ define_explicitonestepmethod("eosm")
+ return "eosm"
diff --git a/python/dune/perftool/pdelab/driver/interpolate.py b/python/dune/perftool/pdelab/driver/interpolate.py
new file mode 100644
index 00000000..424388a0
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/interpolate.py
@@ -0,0 +1,115 @@
+from dune.perftool.generation import (cached,
+ get_counted_variable,
+ global_context,
+ include_file,
+ preamble,
+ )
+from dune.perftool.pdelab.driver import (FEM_name_mangling,
+ get_formdata,
+ get_trial_element,
+ is_stationary,
+ preprocess_leaf_data,
+ )
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
+ name_leafview,
+ )
+from dune.perftool.pdelab.driver.gridoperator import (name_parameters,)
+
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+
+
+def _do_interpolate(dirichlet):
+ if isinstance(dirichlet, (list, tuple)):
+ return any(bool(d) for d in dirichlet)
+ else:
+ return bool(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)
+ dirichlet = preprocess_leaf_data(element, "dirichlet_expression")
+ bf = name_boundary_function(element, dirichlet)
+ interpolate_vector(bf, gfs, name)
+
+
+@preamble
+def interpolate_vector(func, gfs, name):
+ return "Dune::PDELab::interpolate({}, {}, {});".format(func,
+ gfs,
+ name,
+ )
+
+
+def name_boundary_function(element, dirichlet):
+ if isinstance(element, MixedElement):
+ k = 0
+ childs = []
+ for subel in element.sub_elements():
+ childs.append(name_boundary_function(subel, dirichlet[k:k + subel.value_size()]))
+ k = k + subel.value_size()
+ name = "_".join(childs)
+ define_composite_boundary_function(name, childs)
+ return name
+ else:
+ assert isinstance(element, FiniteElement)
+ name = "{}_boundary".format(FEM_name_mangling(element).lower())
+ define_boundary_function(name, dirichlet[0])
+ return name
+
+
+@preamble
+def define_compositegfs_parameterfunction(name, children):
+ return "Dune::PDELab::CompositeGridFunction<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in children),
+ name,
+ ', '.join(children))
+
+
+@preamble
+def define_boundary_function(name, dirichlet):
+ gv = name_leafview()
+ lambdaname = name_boundary_lambda(dirichlet)
+ include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
+ if is_stationary():
+ return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
+ gv,
+ lambdaname,
+ )
+ else:
+ params = name_parameters(get_formdata())
+ return "auto {} = Dune::PDELab::makeInstationaryGridFunctionFromCallable({}, {}, {});".format(name,
+ gv,
+ lambdaname,
+ params)
+
+
+@cached
+def name_boundary_lambda(boundary):
+ name = get_counted_variable("lambda")
+ define_boundary_lambda(name, boundary)
+ return name
+
+
+@preamble
+def define_boundary_lambda(name, boundary):
+ from dune.perftool.ufl.execution import Expression
+ from ufl.classes import Expr
+ if boundary is None:
+ return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
+ elif isinstance(boundary, Expression):
+ if boundary.is_global:
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr[0])
+ else:
+ return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr[0])
+ elif isinstance(boundary, Expr):
+ # Set up a visitor
+ with global_context(integral_type="exterior_facet", formdata=get_formdata(), driver=True):
+ from dune.perftool.ufl.visitor import UFL2LoopyVisitor
+ from dune.perftool.pdelab import PDELabInterface
+ visitor = UFL2LoopyVisitor(PDELabInterface(), "exterior_facet", {})
+ from dune.perftool.loopy.target import numpy_to_cpp_dtype
+ return "auto {} = [&](const auto& x){{ return ({}){}; }};".format(name,
+ numpy_to_cpp_dtype("float64"),
+ visitor(boundary))
diff --git a/python/dune/perftool/pdelab/driver/solve.py b/python/dune/perftool/pdelab/driver/solve.py
new file mode 100644
index 00000000..5e602827
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/solve.py
@@ -0,0 +1,245 @@
+from dune.perftool.generation import (include_file,
+ preamble,
+ )
+from dune.perftool.options import get_option
+from dune.perftool.pdelab.driver import (form_name_suffix,
+ get_formdata,
+ is_linear,
+ name_initree,
+ )
+from dune.perftool.pdelab.driver.gridfunctionspace import name_trial_gfs
+from dune.perftool.pdelab.driver.gridoperator import (name_gridoperator,
+ type_gridoperator,
+ )
+from dune.perftool.pdelab.driver.interpolate import interpolate_dirichlet_data
+
+
+@preamble
+def dune_solve():
+ # Test if form is linear in ansatzfunction
+ linear = is_linear()
+
+ # Test wether we want to do matrix free operator evaluation
+ matrix_free = get_option('matrix_free')
+
+ # Get right solve command
+ if linear and matrix_free:
+ formdata = get_formdata()
+ go = name_gridoperator(formdata)
+ x = name_vector(formdata)
+ include_file("dune/perftool/matrixfree.hh", filetag="driver")
+ solve = "solveMatrixFree({},{});".format(go, x)
+ elif linear and not matrix_free:
+ slp = name_stationarylinearproblemsolver()
+ solve = "{}.apply();".format(slp)
+ elif not linear and matrix_free:
+ # TODO copy of linear case and obviously broken, used to generate something ;)
+ formdata = get_formdata()
+ go = name_gridoperator(formdata)
+ x = name_vector(formdata)
+ include_file("dune/perftool/matrixfree.hh", filetag="driver")
+ solve = "solveNonlinearMatrixFree({},{});".format(go, x)
+ elif not linear and not matrix_free:
+ go_type = type_gridoperator(_driver_data['formdata'])
+ go = name_gridoperator(_driver_data['formdata'])
+ snp = name_stationarynonlinearproblemsolver(go_type, go)
+ solve = "{}.apply();".format(snp)
+
+ 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()
+
+ if get_option('instrumentation_level') >= 2:
+ from dune.perftool.pdelab.driver.timings import setup_timer, name_timing_stream
+ setup_timer()
+ from dune.perftool.generation import post_include
+ post_include("HP_DECLARE_TIMER(solve);", filetag="driver")
+
+ # Print times after solving
+ from dune.perftool.generation import get_global_context_value
+ formdatas = get_global_context_value("formdatas")
+ print_times = []
+ for formdata in formdatas:
+ from dune.perftool.pdelab.driver.gridoperator import name_localoperator
+ lop_name = name_localoperator(formdata)
+ timestream = name_timing_stream()
+ print_times.append("{}.dump_timers({}, {}, true);".format(lop_name, timestream, name_timing_identifier()))
+
+ solve = ["HP_TIMER_START(solve);",
+ "{}".format(solve),
+ "HP_TIMER_STOP(solve);",
+ "DUMP_TIMER(solve, {}, true);".format(timestream),
+ ]
+ if get_option('instrumentation_level') >= 3:
+ solve.extend(print_times)
+
+ return solve
+
+
+def name_vector(formdata):
+ name = form_name_suffix("x", formdata)
+ define_vector(name, formdata)
+ interpolate_dirichlet_data(name)
+ return name
+
+
+@preamble
+def typedef_vector(name, formdata):
+ go_type = type_gridoperator(formdata)
+ return "using {} = {}::Traits::Domain;".format(name, go_type)
+
+
+def type_vector(formdata):
+ name = form_name_suffix("V", formdata).upper()
+ typedef_vector(name, formdata)
+ return name
+
+
+@preamble
+def define_vector(name, formdata):
+ vtype = type_vector(formdata)
+ gfs = name_trial_gfs()
+ return ["{} {}({});".format(vtype, name, gfs), "{} = 0.0;".format(name)]
+
+
+
+@preamble
+def typedef_linearsolver(name):
+ include_file("dune/pdelab/backend/istl.hh", filetag="driver")
+ return "using {} = Dune::PDELab::ISTLBackend_SEQ_UMFPack;".format(name)
+
+
+def type_linearsolver():
+ name = "LinearSolver"
+ typedef_linearsolver(name)
+ return name
+
+
+@preamble
+def define_linearsolver(name):
+ lstype = type_linearsolver()
+ return "{} {}(false);".format(lstype, name)
+
+
+def name_linearsolver():
+ name = "ls"
+ define_linearsolver(name)
+ return name
+
+
+@preamble
+def define_reduction(name):
+ ini = name_initree()
+ return "double {} = {}.get<double>(\"reduction\", 1e-12);".format(name, ini)
+
+
+def name_reduction():
+ name = "reduction"
+ define_reduction(name)
+ return name
+
+
+@preamble
+def typedef_stationarylinearproblemsolver(name):
+ include_file("dune/pdelab/stationary/linearproblem.hh", filetag="driver")
+ gotype = type_gridoperator(get_formdata())
+ lstype = type_linearsolver()
+ xtype = type_vector(get_formdata())
+ return "using {} = Dune::PDELab::StationaryLinearProblemSolver<{}, {}, {}>;".format(name, gotype, lstype, xtype)
+
+
+def type_stationarylinearproblemsolver():
+ typedef_stationarylinearproblemsolver("SLP")
+ return "SLP"
+
+
+@preamble
+def define_stationarylinearproblemsolver(name):
+ slptype = type_stationarylinearproblemsolver()
+ formdata = get_formdata()
+ go = name_gridoperator(formdata)
+ ls = name_linearsolver()
+ x = name_vector(formdata)
+ red = name_reduction()
+ return "{} {}({}, {}, {}, {});".format(slptype, name, go, ls, x, red)
+
+
+def name_stationarylinearproblemsolver():
+ define_stationarylinearproblemsolver("slp")
+ return "slp"
+
+
+@preamble
+def typedef_stationarynonlinearproblemsolver(name, go_type):
+ include_file("dune/pdelab/newton/newton.hh", filetag="driver")
+ ls_type = type_linearsolver()
+ x_type = type_vector(_driver_data['formdata'])
+ return "using {} = Dune::PDELab::Newton<{}, {}, {}>;".format(name, go_type, ls_type, x_type)
+
+
+def type_stationarynonlinearproblemssolver(go_type):
+ name = "SNP"
+ typedef_stationarynonlinearproblemsolver(name, go_type)
+ return name
+
+
+@preamble
+def define_stationarynonlinearproblemsolver(name, go_type, go):
+ snptype = type_stationarynonlinearproblemssolver(go_type)
+ x = name_vector(get_formdata())
+ ls = name_linearsolver()
+ return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
+
+
+def name_stationarynonlinearproblemsolver(go_type, go):
+ name = "snp"
+ define_stationarynonlinearproblemsolver(name, go_type, go)
+ return name
+
+
+@preamble
+def print_residual():
+ ini = name_initree()
+ formdata = get_formdata()
+ n_go = name_gridoperator(formdata)
+ v = name_vector(formdata)
+ t_v = type_vector(formdata)
+ include_file("random", system=True, filetag="driver")
+
+ return ["if ({}.get<bool>(\"printresidual\", false)) {{".format(ini),
+ " using Dune::PDELab::Backend::native;",
+ " {} r({});".format(t_v, v),
+ " // Setup random input",
+ " std::size_t seed = 0;",
+ " auto rng = std::mt19937_64(seed);",
+ " auto dist = std::uniform_real_distribution<>(-1., 1.);",
+ " for (auto& v : {})".format(v),
+ " v = dist(rng);",
+ " r=0.0;",
+ " {}.residual({}, r);".format(n_go, v),
+ " Dune::printvector(std::cout, native(r), \"residual vector\", \"row\");",
+ "}"]
+
+
+@preamble
+def print_matrix():
+ formdata = get_formdata()
+ ini = name_initree()
+ t_go = type_gridoperator(formdata)
+ n_go = name_gridoperator(formdata)
+ v = name_vector(formdata)
+ t_v = type_vector(formdata)
+
+ return ["if ({}.get<bool>(\"printmatrix\", false)) {{".format(ini),
+ " using M = typename {}::Traits::Jacobian;".format(t_go),
+ " M m({});".format(n_go),
+ " {}.jacobian({},m);".format(n_go, v),
+ " using Dune::PDELab::Backend::native;",
+ " Dune::printmatrix(std::cout, native(m),\"global stiffness matrix\",\"row\",9,1);",
+ "}"]
diff --git a/python/dune/perftool/pdelab/driver/vtk.py b/python/dune/perftool/pdelab/driver/vtk.py
new file mode 100644
index 00000000..2f68140e
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/vtk.py
@@ -0,0 +1,151 @@
+from dune.perftool.generation import (include_file,
+ preamble,
+ )
+from dune.perftool.pdelab.driver import (get_formdata,
+ get_trial_element,
+ name_initree,
+ )
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
+ name_trial_gfs,
+ type_leafview,
+ )
+from dune.perftool.pdelab.driver.solve import name_vector
+
+
+@preamble
+def define_vtkfile(name):
+ ini = name_initree()
+ include_file("string", filetag="driver")
+ return "std::string {} = {}.get<std::string>(\"wrapper.vtkcompare.name\", \"output\");".format(name, ini)
+
+
+def name_vtkfile():
+ define_vtkfile("vtkfile")
+ return "vtkfile"
+
+
+@preamble
+def typedef_vtkwriter(name):
+ include_file("dune/grid/io/file/vtk/subsamplingvtkwriter.hh", filetag="driver")
+ gv = type_leafview()
+ return "using {} = Dune::SubsamplingVTKWriter<{}>;".format(name, gv)
+
+
+def type_vtkwriter():
+ typedef_vtkwriter("VTKWriter")
+ return "VTKWriter"
+
+
+@preamble
+def define_subsamplinglevel(name):
+ ini = name_initree()
+ return "int {} = {}.get<int>(\"vtk.subsamplinglevel\", 0);".format(name, ini)
+
+
+def name_subsamplinglevel():
+ define_subsamplinglevel("sublevel")
+ return "sublevel"
+
+
+@preamble
+def define_vtkwriter(name):
+ _type = type_vtkwriter()
+ gv = name_leafview()
+ subsamp = name_subsamplinglevel()
+ return "{} {}({}, {});".format(_type, name, gv, subsamp)
+
+
+def name_vtkwriter():
+ define_vtkwriter("vtkwriter")
+ return "vtkwriter"
+
+
+@preamble
+def vtkoutput():
+ include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
+ vtkwriter = name_vtkwriter()
+ gfs = name_trial_gfs()
+ vtkfile = name_vtkfile()
+ predicate = name_predicate()
+ vec = name_vector(get_formdata())
+
+ return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate),
+ "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]
+
+
+def type_predicate():
+ include_file("dune/perftool/vtkpredicate.hh", filetag="driver")
+ return "CuttingPredicate"
+
+
+@preamble
+def define_predicate(name):
+ t = type_predicate()
+ return "{} {};".format(t, name)
+
+
+def name_predicate():
+ define_predicate("predicate")
+ return "predicate"
+
+
+
+@preamble
+def typedef_vtk_sequence_writer(name):
+ include_file("dune/grid/io/file/vtk/vtksequencewriter.hh", filetag="driver")
+ gv_type = type_leafview()
+ return "using {} = Dune::VTKSequenceWriter<{}>;".format(name, gv_type)
+
+
+def type_vtk_sequence_writer():
+ typedef_vtk_sequence_writer("VTKSW")
+ return "VTKSW"
+
+
+@preamble
+def define_vtk_sequence_writer(name):
+ vtksw_type = type_vtk_sequence_writer()
+ vtkw_type = type_vtkwriter()
+ vtkw = name_vtkwriter()
+ vtkfile = name_vtkfile()
+ return "{} {}(std::make_shared<{}>({}), {});".format(vtksw_type, name, vtkw_type, vtkw, vtkfile)
+
+
+def name_vtk_sequence_writer():
+ define_vtk_sequence_writer("vtkSequenceWriter")
+ return "vtkSequenceWriter"
+
+
+@preamble
+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_gfs(element)
+ vector = name_vector(get_formdata())
+ predicate = name_predicate()
+ from dune.perftool.pdelab.driver.instationary import name_time
+ 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
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 9e857d97..0f753451 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -9,16 +9,6 @@ import ufl
from ufl import *
-# Monkey patch to the FiniteElementBase
-# We want to subclass finite elements in order to allow meta data
-# connected to constraints. This is an intermediate solution until
-# we have a proper notion of (global) function spaces in UFL. As is,
-# two finite elements with different meta data would not compare equally,
-# but they should in order to allow hashing generation magic based on
-# them (as done everywhere)
-ufl.FiniteElementBase.__eq__ = lambda s, o: repr(s) == repr(o)
-
-
class TrialFunction(ufl.Coefficient):
""" A coefficient that always takes the reserved index 0 """
def __init__(self, element, count=None):
@@ -136,51 +126,3 @@ def split_expression(expr):
element=expr.element.sub_elements()[i])
for i in range(expr.element.num_sub_elements())
)
-
-
-class FiniteElement(ufl.FiniteElement):
-
- __slots = ufl.FiniteElement.__slots__ + ["dirichlet_expression", "dirichlet_constraints"]
-
- def __new__(cls, *args, **kwargs):
- kwargs.pop('dirichlet_constraints', None)
- kwargs.pop('dirichlet_expression', None)
- return ufl.FiniteElement.__new__(cls, *args, **kwargs)
-
- def __init__(self, *args, **kwargs):
- if ('dirichlet_constraints' in kwargs) or ('dirichlet_expression' in kwargs):
- # Get dirichlet_constraints and convert it to Expression if necessary!
- self.dirichlet_constraints = kwargs.pop('dirichlet_constraints', 'return true;')
- if isinstance(self.dirichlet_constraints, str):
- self.dirichlet_constraints = Expression(self.dirichlet_constraints)
- from ufl.classes import Expr
- assert isinstance(self.dirichlet_constraints, (Expression, Expr))
-
- # Get dirichlet_constraints and convert it to Expression if necessary!
- self.dirichlet_expression = kwargs.pop('dirichlet_expression', 'return 0.0;')
- if isinstance(self.dirichlet_expression, str):
- self.dirichlet_expression = Expression(self.dirichlet_expression)
- assert isinstance(self.dirichlet_expression, (Expression, Expr))
-
- # Initialize the original finite element from ufl
- ufl.FiniteElement.__init__(self, *args, **kwargs)
-
-
-class VectorElement(ufl.VectorElement):
- def __init__(self, *args, **kwargs):
- if ('dirichlet_constraints' in kwargs) or ('dirichlet_expression' in kwargs):
- # Get dirichlet_constraints and convert it to Expression if necessary!
- self.dirichlet_constraints = kwargs.pop('dirichlet_constraints', 'return true;')
- if isinstance(self.dirichlet_constraints, str):
- self.dirichlet_constraints = Expression(self.dirichlet_constraints)
- from ufl.classes import Expr
- assert isinstance(self.dirichlet_constraints, (Expression, Expr))
-
- # Get dirichlet_constraints and convert it to Expression if necessary!
- self.dirichlet_expression = kwargs.pop('dirichlet_expression', 'return 0.0;')
- if isinstance(self.dirichlet_expression, str):
- self.dirichlet_expression = Expression(self.dirichlet_expression)
- assert isinstance(self.dirichlet_expression, (Expression, Expr))
-
- # Initialize the original finite element from ufl
- ufl.VectorElement.__init__(self, *args, **kwargs)
diff --git a/test/poisson/poisson.mini b/test/poisson/poisson.mini
index 063be404..3597ac05 100644
--- a/test/poisson/poisson.mini
+++ b/test/poisson/poisson.mini
@@ -13,5 +13,4 @@ extension = vtu
[formcompiler]
numerical_jacobian = 1, 0 | expand num
-exact_solution_expression = g
compare_l2errorsquared = 1e-7
diff --git a/test/poisson/poisson.ufl b/test/poisson/poisson.ufl
index 73515675..2bfe3313 100644
--- a/test/poisson/poisson.ufl
+++ b/test/poisson/poisson.ufl
@@ -6,9 +6,12 @@ c = (0.5-x[0])**2 + (0.5-x[1])**2
g = exp(-1.*c)
f = 4*(1.-c)*g
-V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
+V = FiniteElement("CG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
forms = [(inner(grad(u), grad(v)) - f*v)*dx]
+exact_solution = g
+dirichlet_expression = g
+is_dirichlet = 1
\ No newline at end of file
--
GitLab