From 10bd5d43e7a0ff9212ffb24e3f7a9dd5c177e82b Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Tue, 20 Jun 2017 13:42:31 +0200
Subject: [PATCH] Introduce subpackage driver in pdelab
---
.../dune/perftool/pdelab/driver/__init__.py | 1745 +++++++++++++++--
1 file changed, 1626 insertions(+), 119 deletions(-)
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index b5efebf3..48fa4574 100644
--- a/python/dune/perftool/pdelab/driver/__init__.py
+++ b/python/dune/perftool/pdelab/driver/__init__.py
@@ -43,30 +43,6 @@ def set_driver_data(formdatas, data):
_driver_data['data'] = data
-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()
-
-
-def get_trial_element():
- return _driver_data['form'].coefficients()[0].ufl_element()
-
-
-def get_formdata():
- return _driver_data['formdata']
-
-
-def get_mass_formdata():
- return _driver_data["mass_formdata"]
-
-
def is_stationary():
return 'mass_form' not in _driver_data
@@ -78,8 +54,29 @@ def form_name_suffix(name, formdata):
return name + '_' + form_name
-def get_object(name):
- return _driver_data['data'].object_by_name.get(name, None)
+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 mass_form_index(formdatas, data):
@@ -91,10 +88,8 @@ def mass_form_index(formdatas, data):
continue
-def is_linear(form=None):
- '''Test if form is linear in trial function'''
- if form is None:
- form = get_formdata().original_form
+def is_linear(form):
+ '''Test if form is linear in test function'''
from ufl import derivative
from ufl.algorithms import expand_derivatives
jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
@@ -109,26 +104,20 @@ def isLagrange(fem):
return fem._short_name is 'CG'
-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 isSimplical(fem):
+ return any(fem._cell._cellname in x for x in ["triangle", "tetrahedron"])
-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 isQuadrilateral(fem):
+ return any(fem._cell._cellname in x for x in ["vertex", "interval", "quadrilateral", "hexahedron"])
def isPk(fem):
- return isLagrange(fem) and isSimplical(fem.cell())
+ return isLagrange(fem) and isSimplical(fem)
def isQk(fem):
- return isLagrange(fem) and isQuadrilateral(fem.cell())
+ return isLagrange(fem) and isQuadrilateral(fem)
def isDG(fem):
@@ -159,72 +148,1000 @@ 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
+def name_inifile():
+ # TODO pass some other option here.
+ return "argv[1]"
+
+
+@preamble
+def parse_initree(varname):
+ include_file("dune/common/parametertree.hh", filetag="driver")
+ include_file("dune/common/parametertreeparser.hh", filetag="driver")
+ filename = name_inifile()
+ return ["Dune::ParameterTree initree;", "Dune::ParameterTreeParser::readINITree({}, {});".format(filename, varname)]
+
+
+def name_initree():
+ parse_initree("initree")
+ # TODO we can get some other ini file here.
+ return "initree"
+
+
+@preamble
+def define_dimension(name):
+ return "static const int {} = {};".format(name, _driver_data['form'].ufl_cell().geometric_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:
- yield l
+ 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 _unroll_list_tensors(expr):
- from ufl.classes import ListTensor
- if isinstance(expr, ListTensor):
- for op in expr.ufl_operands:
- yield op
+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"
+
+
+@preamble
+def typedef_range(name):
+ if get_option('opcounter'):
+ setup_timer()
+ return "using {} = oc::OpCounter<double>;".format(name)
else:
- yield expr
+ return "using {} = double;".format(name)
-def unroll_list_tensors(data):
- for expr in data:
- for e in _unroll_list_tensors(expr):
- yield e
+def type_range():
+ typedef_range("RangeType")
+ return "RangeType"
-def preprocess_leaf_data(element, data):
- data = get_object(data)
- from ufl import MixedElement
+@preamble
+def typedef_fem(expr, name):
+ gv = type_leafview()
+ df = type_domainfield()
+ r = type_range()
+ dim = name_dimension()
+ 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):
- # data is None -> use 0 default
- if data is None:
- data = (0,) * element.value_size()
+ 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
+
+
+@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(_driver_data['form'].coefficients()[0].ufl_element()):
+ geo_factor = "4"
+ else:
+ geo_factor = "6"
+ gfs = name_gfs(_driver_data['form'].coefficients()[0].ufl_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)
- # Flatten nested lists
- data = tuple(i for i in _flatten_list(data))
- # Expand any list tensors
- data = tuple(i for i in unroll_list_tensors(data))
+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(_driver_data['form'].coefficients()[0].ufl_element())
+ vgfs = type_gfs(_driver_data['form'].arguments()[0].ufl_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)
- assert len(data) == element.value_size()
- return data
+
+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(_driver_data['form'].coefficients()[0].ufl_element())
+ vgfs = type_gfs(_driver_data['form'].arguments()[0].ufl_element())
+ lop = type_localoperator(formdata)
+ ucc = type_constraintscontainer(_driver_data['form'].coefficients()[0].ufl_element())
+ vcc = type_constraintscontainer(_driver_data['form'].arguments()[0].ufl_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:
- # Do not return lists for non-MixedElement
- if not isinstance(data, (tuple, list)):
- return (data,)
+ 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(_driver_data['form'].coefficients()[0].ufl_element())
+ ucc = name_assembled_constraints(_driver_data['form'].coefficients()[0].ufl_element())
+ vgfs = name_gfs(_driver_data['form'].arguments()[0].ufl_element())
+ vcc = name_assembled_constraints(_driver_data['form'].arguments()[0].ufl_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(_driver_data['form'].coefficients()[0].ufl_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:
- assert len(data) == 1
- return data
+ 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_inifile():
- # TODO pass some other option here.
- return "argv[1]"
+
+def name_boundary_lambda(boundary, name):
+ define_boundary_lambda(boundary, name + "lambda")
+ return name + "lambda"
@preamble
-def parse_initree(varname):
- include_file("dune/common/parametertree.hh", filetag="driver")
- include_file("dune/common/parametertreeparser.hh", filetag="driver")
- filename = name_inifile()
- return ["Dune::ParameterTree initree;", "Dune::ParameterTreeParser::readINITree({}, {});".format(filename, varname)]
+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)
-def name_initree():
- parse_initree("initree")
- # TODO we can get some other ini file here.
- return "initree"
+@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)
+
+
+@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'))
+
+ define_boundary_function(expr, name)
+
+ return name
+
+
+@preamble
+def interpolate_vector(name, formdata):
+ define_vector(name, formdata)
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ bf = name_boundary_function(element)
+ gfs = name_gfs(element)
+ return "Dune::PDELab::interpolate({}, {}, {});".format(bf,
+ gfs,
+ name,
+ )
+
+
+@preamble
+def interpolate_solution_expression(name):
+ formdata = _driver_data['formdata']
+ define_vector(name, formdata)
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ sol = name_solution_function()
+ gfs = name_gfs(element)
+ return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
+ gfs,
+ name,
+ )
+
+
+def maybe_interpolate_vector(name, formdata):
+ element = _driver_data['form'].coefficients()[0].ufl_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
+
+
+def name_vector_from_solution_expression():
+ interpolate_solution_expression("solution")
+ return "solution"
+
+
+@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
@@ -239,34 +1156,638 @@ def name_mpihelper():
return name
-def _blockstructured_degree_helper(element, op):
- element._degree = op(element._degree, get_option("number_of_blocks"))
+@preamble
+def define_timing_stream(name):
+ include_file('fstream', filetag='driver', system=True)
+ include_file('sstream', filetag='driver', system=True)
+ include_file('sys/types.h', filetag='driver', system=True)
+ include_file('unistd.h', filetag='driver', system=True)
+
+ return ["std::stringstream ss;",
+ "ss << \"{}/timings-rank-\" << {}.rank() << \"-pid-\" << getpid() << \".csv\";".format(get_option('project_basedir'), name_mpihelper()),
+ "std::ofstream {};".format(name),
+ "{}.open(ss.str(), std::ios_base::app);".format(name),
+ ]
+
+
+@preamble
+def dump_dof_numbers(stream):
+ ident = name_timing_identifier()
+ return "{} << {} << \" dofs dofs \" << {}.size() << std::endl;".format(stream,
+ ident,
+ name_gfs(_driver_data['form'].coefficients()[0].ufl_element()))
+
+
+def name_timing_stream():
+ name = "timestream"
+ define_timing_stream(name)
+ dump_dof_numbers(name)
+ 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:
+ 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 compare_dofs():
+ v = name_vector(_driver_data['formdata'])
+ solution = name_vector_from_solution_expression()
+ fail = name_test_fail_variable()
+ return ["",
+ "// Maximum norm of difference between dof vectors of the",
+ "// numerical solution and the interpolation of the exact solution",
+ "using Dune::PDELab::Backend::native;",
+ "double maxerror(0.0);",
+ "for (std::size_t i=0; i<native({}).size(); ++i)".format(v),
+ " if (std::abs(native({})[i]-native({})[i]) > maxerror)".format(v, solution),
+ " maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
+ "std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
+ "if (maxerror>{})".format(get_option("compare_dofs")),
+ " {} = true;".format(fail)]
+
+
+def type_discrete_grid_function(gfs):
+ return "DGF_{}".format(gfs.upper())
+
+
+@preamble
+def define_discrete_grid_function(gfs, vector_name, dgf_name):
+ dgf_type = type_discrete_grid_function(gfs)
+ return ["using {} = Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})>;".format(dgf_type, gfs, vector_name),
+ "{} {}({},{});".format(dgf_type, dgf_name, gfs, vector_name)]
+
+
+def name_discrete_grid_function(gfs, vector_name):
+ dgf_name = gfs + "_dgf"
+ define_discrete_grid_function(gfs, vector_name, dgf_name)
+ return dgf_name
+
+
+def type_subgfs(element, tree_path):
+ 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)
+ gfs = name_gfs(element)
+
+ return '{} {}({});'.format(t, name, gfs)
+
+
+def name_subgfs(element, tree_path):
+ 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
+
+
+@preamble
+def typedef_difference_squared_adapter(name, tree_path):
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ formdata = _driver_data['formdata']
+
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, solution_function, dgf)
+
+
+def type_difference_squared_adapter(tree_path):
+ t = 'DifferenceSquared_{}'.format('_'.join(tree_path))
+ typedef_difference_squared_adapter(t, tree_path)
+ return t
+
+
+@preamble
+def define_difference_squared_adapter(name, tree_path):
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ formdata = _driver_data['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)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return '{} {}({}, {});'.format(t, name, solution_function, dgf)
+
+
+def name_difference_squared_adapter(tree_path):
+ name = 'differencesquared_{}'.format('_'.join(tree_path))
+ define_difference_squared_adapter(name, tree_path)
+ return name
+
+
+@preamble
+def _accumulate_L2_squared(tree_path):
+ dsa = name_difference_squared_adapter(tree_path)
+ t_dsa = type_difference_squared_adapter(tree_path)
+ 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=()):
+ element = _driver_data['form'].coefficients()[0].ufl_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)
+
+
+@preamble
+def define_accumulated_L2_error(name):
+ return "double {}(0.0);".format(name)
+
+
+def name_accumulated_L2_error():
+ name = 'l2error'
+ define_accumulated_L2_error(name)
+ return name
+
+
+@preamble
+def compare_L2_squared():
+ accumulate_L2_squared()
+
+ accum_error = name_accumulated_L2_error()
+ fail = name_test_fail_variable()
+ return ["std::cout << \"l2errorsquared: \" << {} << std::endl;".format(accum_error),
+ "if (std::isnan({0}) or std::abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
+ " {} = true;".format(fail)]
+
+
+@preamble
+def define_test_fail_variable(name):
+ return 'bool {}(false);'.format(name)
+
+
+def name_test_fail_variable():
+ name = "testfail"
+ define_test_fail_variable(name)
+ return name
+
+
+@cached
+def setup_timer():
+ # Necessary includes and defines
+ from dune.perftool.generation import pre_include
+
+ # TODO check that we are using YASP?
+ if get_option('opcounter'):
+ pre_include("#define ENABLE_COUNTER", filetag="driver")
+ pre_include("#define ENABLE_HP_TIMERS", filetag="driver")
+ include_file("dune/perftool/common/timer.hh", filetag="driver")
+
+
+@preamble
+def define_timing_identifier(name):
+ ini = name_initree()
+ return "auto {} = {}.get<std::string>(\"identifier\", std::string(argv[0]));".format(name, ini)
+
+
+def name_timing_identifier():
+ name = "ident"
+ define_timing_identifier(name)
+ return name
+
+
+@preamble
+def evaluate_residual_timer():
+ formdata = _driver_data['formdata']
+ n_go = name_gridoperator(formdata)
+ v = name_vector(formdata)
+ t_v = type_vector(formdata)
+ setup_timer()
+
+ if get_option('instrumentation_level') >= 2:
+ # Write back times
+ from dune.perftool.generation import post_include
+ post_include("HP_DECLARE_TIMER(residual_evaluation);", filetag="driver")
+ timestream = name_timing_stream()
+ print_times = []
+
+ from dune.perftool.generation import get_global_context_value
+ formdatas = get_global_context_value("formdatas")
+ for formdata in formdatas:
+ lop_name = name_localoperator(formdata)
+ if get_option('instrumentation_level') >= 3:
+ print_times.append("{}.dump_timers({}, {}, true);".format(lop_name, timestream, name_timing_identifier()))
+
+ if get_option('instrumentation_level') >= 2:
+ evaluation = ["HP_TIMER_START(residual_evaluation);",
+ "{}.residual({}, r);".format(n_go, v),
+ "HP_TIMER_STOP(residual_evaluation);",
+ "DUMP_TIMER(residual_evaluation, {}, true);".format(timestream)]
+ evaluation.extend(print_times)
+ else:
+ evaluation = ["{}.residual({}, r);".format(n_go, v)]
+
+ evaluation = ["{} r({});".format(t_v, v), "r=0.0;"] + evaluation
+
+ return evaluation
+
+
+@preamble
+def apply_jacobian_timer():
+ # Set the matrix_free option to True!
+ from dune.perftool.options import set_option
+ set_option("matrix_free", True)
+
+ formdata = _driver_data['formdata']
+ n_go = name_gridoperator(formdata)
+ v = name_vector(formdata)
+ t_v = type_vector(formdata)
+ setup_timer()
+
+ if get_option('instrumentation_level') >= 2:
+ # Write back times
+ from dune.perftool.generation import post_include
+ post_include("HP_DECLARE_TIMER(apply_jacobian);", filetag="driver")
+ timestream = name_timing_stream()
+ print_times = []
+
+ from dune.perftool.generation import get_global_context_value
+ formdatas = get_global_context_value("formdatas")
+ for formdata in formdatas:
+ lop_name = name_localoperator(formdata)
+ if get_option('instrumentation_level') >= 3:
+ print_times.append("{}.dump_timers({}, {}, true);".format(lop_name, timestream, name_timing_identifier()))
+
+ if get_option('instrumentation_level') >= 2:
+ evaluation = ["HP_TIMER_START(apply_jacobian);",
+ "{}.jacobian_apply({}, j);".format(n_go, v),
+ "HP_TIMER_STOP(apply_jacobian);",
+ "DUMP_TIMER(apply_jacobian, {}, true);".format(timestream)]
+ evaluation.extend(print_times)
+ else:
+ evaluation = ["{}.jacobian_apply({}, j);".format(n_go, v)]
+
+ evaluation = ["{} j({});".format(t_v, v), "j=0.0;"] + evaluation
+
+ return evaluation
+
+
+@preamble
+def assemble_matrix_timer():
+ formdata = _driver_data['formdata']
+ t_go = type_gridoperator(formdata)
+ n_go = name_gridoperator(formdata)
+ v = name_vector(formdata)
+ t_v = type_vector(formdata)
+ setup_timer()
+
+ if get_option('instrumentation_level') >= 2:
+ # Write back times
+ from dune.perftool.generation import post_include
+ post_include("HP_DECLARE_TIMER(matrix_assembly);", filetag="driver")
+ timestream = name_timing_stream()
+ print_times = []
+
+ from dune.perftool.generation import get_global_context_value
+ formdatas = get_global_context_value("formdatas")
+ for formdata in formdatas:
+ lop_name = name_localoperator(formdata)
+ if get_option('instrumentation_level') >= 3:
+ print_times.append("{}.dump_timers({}, {}, true);".format(lop_name, timestream, name_timing_identifier()))
+
+ if get_option('instrumentation_level') >= 2:
+ assembly = ["HP_TIMER_START(matrix_assembly);",
+ "{}.jacobian({},m);".format(n_go, v),
+ "HP_TIMER_STOP(matrix_assembly);",
+ "DUMP_TIMER(matrix_assembly, {}, true);".format(timestream)]
+ assembly.extend(print_times)
+ else:
+ assembly = ["{}.jacobian({},m);".format(n_go, v)]
+
+ assembly = ["using M = typename {}::Traits::Jacobian;".format(t_go),
+ "M m({});".format(n_go)] + assembly
+
+ return assembly
+
+
+@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) or isinstance(element, TensorElement):
- _blockstructured_degree_helper(element.sub_elements()[0], op)
- elif isinstance(element, MixedElement):
- for e in element.sub_elements():
- _blockstructured_degree_helper(e, op)
+ 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"
-def set_blockstructured_degree():
- from operator import mul
- _blockstructured_degree_helper(_driver_data["form"].coefficients()[0].ufl_element(), mul)
+
+@preamble
+def define_predicate(name):
+ t = type_predicate()
+ return "{} {};".format(t, name)
-def unset_blockstructured_degree():
- from operator import floordiv
- _blockstructured_degree_helper(_driver_data["form"].coefficients()[0].ufl_element(), floordiv)
+def name_predicate():
+ define_predicate("predicate")
+ return "predicate"
+
+
+@preamble
+def vtkoutput():
+ element = _driver_data['form'].coefficients()[0].ufl_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()
+
+ 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 = _driver_data['form'].coefficients()[0].ufl_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 = _driver_data['form'].coefficients()[0].ufl_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()
+ 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():
+ fail = name_test_fail_variable()
+ return "return {};".format(fail)
def generate_driver(formdatas, data):
# The driver module uses a global dictionary for storing necessary data
set_driver_data(formdatas, data)
- # TODO find a better solution ...
- if get_option("blockstructured"):
- set_blockstructured_degree()
-
# Entrypoint for driver generation
if get_option("opcounter") or get_option("time_opcounter"):
if get_option("time_opcounter"):
@@ -275,29 +1796,19 @@ def generate_driver(formdatas, data):
["vertex", "interval", "quadrilateral", "hexahedron"]))
# In case of operator conunting we only assemble the matrix and evaluate the residual
# assemble_matrix_timer()
- from dune.perftool.pdelab.driver.timings import apply_jacobian_timer, evaluate_residual_timer
evaluate_residual_timer()
apply_jacobian_timer()
elif is_stationary():
- from dune.perftool.pdelab.driver.solve import dune_solve
- vec = dune_solve()
- from dune.perftool.pdelab.driver.vtk import vtkoutput
+ # We could also use solve if we are not interested in visualization
vtkoutput()
else:
- from dune.perftool.pdelab.driver.instationary import solve_instationary
solve_instationary()
- from dune.perftool.pdelab.driver.error import compare_L2_squared
- if get_option("compare_l2errorsquared"):
- compare_L2_squared()
-
# Make sure that timestream is declared before retrieving chache items
if get_option("instrumentation_level") >= 1:
- from dune.perftool.pdelab.driver.timings import setup_timer, name_timing_stream
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
@@ -324,7 +1835,3 @@ def generate_driver(formdatas, data):
# Reset the caching data structure
from dune.perftool.generation import delete_cache_items
delete_cache_items()
-
- # TODO find a better solution ...
- if get_option("blockstructured"):
- unset_blockstructured_degree()
--
GitLab