From 6e0836d14b81f5e2e284967f823e58bc8d46d66f 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 | 1532 +----------------
.../perftool/pdelab/driver/constraints.py | 41 +-
python/dune/perftool/pdelab/driver/error.py | 129 +-
.../pdelab/driver/gridfunctionspace.py | 57 +-
.../perftool/pdelab/driver/gridoperator.py | 9 +-
.../perftool/pdelab/driver/instationary.py | 60 +-
.../perftool/pdelab/driver/interpolate.py | 24 +-
python/dune/perftool/pdelab/driver/solve.py | 14 +-
python/dune/perftool/pdelab/driver/vtk.py | 19 +-
9 files changed, 197 insertions(+), 1688 deletions(-)
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index 47603d74..e246e55e 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,984 +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"
-
-
-@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():
- typedef_range("RangeType")
- return "RangeType"
-
-
-@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):
- 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(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)
-
-
-@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 = get_trial_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 = get_trial_element()
- sol = name_solution_function()
- gfs = name_gfs(element)
- return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
- 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
-
-
-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
def define_mpihelper(name):
include_file("dune/common/parallel/mpihelper.hh", filetag="driver")
@@ -1173,468 +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 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 = get_trial_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 = get_trial_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 = 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)
-
-
-@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
-
-
-@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()
-
- 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()
- 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)
@@ -1651,9 +235,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
@@ -1662,6 +249,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
index 2a093a9c..a75016c2 100644
--- a/python/dune/perftool/pdelab/driver/constraints.py
+++ b/python/dune/perftool/pdelab/driver/constraints.py
@@ -1,13 +1,11 @@
-from dune.perftool.generation import (global_context,
- include_file,
+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_gfs,
- name_leafview,
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
name_trial_gfs,
type_range,
type_trial_gfs,
@@ -29,31 +27,28 @@ def assemble_constraints(name):
element = get_trial_element()
gfs = name_trial_gfs()
is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
- if bool(is_dirichlet[0]):
- bctype_function = name_bctype_function(element, is_dirichlet)
- return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
- gfs,
- name,
- )
- else:
- return []
+ 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_elements()[0]
- child = name_bctype_function(subel, is_dirichlet[:subel.value_size()])
- name = "{}_pow{}bctype".format(child, element.num_sub_elements())
- define_power_bctype_function(element, name, child)
+ 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
- childs = []
+ subgfs = []
for subel in element.sub_elements():
- childs.append(name_bctype_function(subel, is_dirichlet[k:k + subel.value_size()]))
+ subgfs.append(name_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
k = k + subel.value_size()
- name = "{}_bctype".format("_".join(childs))
- define_composite_bctype_function(element, is_dirichlet, name, tuple(childs))
+ name = "_".join(subgfs)
+ define_composite_bctype_function(element, is_dirichlet, name, subgfs)
return name
else:
assert isinstance(element, FiniteElement)
@@ -96,11 +91,17 @@ def name_bctype_lambda(name, func):
@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):
diff --git a/python/dune/perftool/pdelab/driver/error.py b/python/dune/perftool/pdelab/driver/error.py
index 2a8f35b4..dbbb0b9c 100644
--- a/python/dune/perftool/pdelab/driver/error.py
+++ b/python/dune/perftool/pdelab/driver/error.py
@@ -9,19 +9,13 @@ from dune.perftool.pdelab.driver import (get_formdata,
get_trial_element,
preprocess_leaf_data,
)
-from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
- name_trial_gfs,
- name_trial_subgfs,
- type_range,
- )
+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,
- dune_solve,
name_vector,
)
-from ufl import MixedElement, TensorElement, VectorElement
@preamble
@@ -35,18 +29,43 @@ def name_test_fail_variable():
return name
-def name_exact_solution_gridfunction(treepath):
+def name_exact_solution():
+ name = "solution"
+ define_vector(name, get_formdata())
+ interpolate_exact_solution(name)
+ return name
+
+
+def interpolate_exact_solution(name):
element = get_trial_element()
- func = preprocess_leaf_data(element, "exact_solution")
- if isinstance(element, MixedElement):
- index = treepath_to_index(element, treepath)
- func = (func[index],)
- element = element.extract_component(index)[1]
- return name_boundary_function(element, func)
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ gfs = name_gfs(element, is_dirichlet)
+ sol = preprocess_leaf_data(element, "exact_solution")
+ func = name_boundary_function(element, sol)
+ interpolate_vector(func, gfs, name)
+
+
+@preamble
+def compare_dofs():
+ v = name_vector(get_formdata())
+ solution = name_exact_solution()
+ fail = name_test_fail_variable()
+
+ return ["",
+ "// Maximum norm of difference between dof vectors of the",
+ "// numerical solution and the interpolation of the exact solution",
+ "using Dune::PDELab::Backend::native;",
+ "double maxerror(0.0);",
+ "for (std::size_t i=0; i<native({}).size(); ++i)".format(v),
+ " if (std::abs(native({})[i]-native({})[i]) > maxerror)".format(v, solution),
+ " maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
+ "std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
+ "if (maxerror>{})".format(get_option("compare_dofs")),
+ " {} = true;".format(fail)]
def type_discrete_grid_function(gfs):
- return "{}_DGF".format(gfs.upper())
+ return "DGF_{}".format(gfs.upper())
@preamble
@@ -63,93 +82,53 @@ def name_discrete_grid_function(gfs, vector_name):
@preamble
-def typedef_difference_squared_adapter(name, treepath):
- sol = name_exact_solution_gridfunction(treepath)
+def typedef_difference_squared_adapter(name):
+ sol = name_exact_solution()
vector = name_vector(get_formdata())
- gfs = name_trial_subgfs(treepath)
dgf = name_discrete_grid_function(gfs, vector)
return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, sol, dgf)
-def type_difference_squared_adapter(treepath):
- name = 'DifferenceSquaredAdapter_{}'.format("_".join(str(t) for t in treepath))
- typedef_difference_squared_adapter(name, treepath)
+def type_difference_squared_adapter():
+ name = 'DifferenceSquaredAdapter'
+ typedef_difference_squared_adapter(name)
return name
@preamble
-def define_difference_squared_adapter(name, treepath):
- t = type_difference_squared_adapter(treepath)
- sol = name_exact_solution_gridfunction(treepath)
+def define_difference_squared_adapter(name):
+ t = type_difference_squared_adapter()
+ sol = name_exact_solution()
vector = name_vector(get_formdata())
- gfs = name_trial_subgfs(treepath)
dgf = name_discrete_grid_function(gfs, vector)
return '{} {}({}, {});'.format(t, name, sol, dgf)
-def name_difference_squared_adapter(treepath):
- name = 'dsa_{}'.format("_".join(str(t) for t in treepath))
- define_difference_squared_adapter(name, treepath)
+def name_difference_squared_adapter():
+ name = 'differencesquared_adapter'
+ define_difference_squared_adapter(name)
return name
@preamble
-def _accumulate_L2_squared(treepath):
- dsa = name_difference_squared_adapter(treepath)
+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")
- strtp = ", ".join(str(t) for t in treepath)
-
- return ["{",
- " // L2 error squared of difference between numerical",
- " // solution and the interpolation of exact solution",
- " // for treepath ({})".format(strtp),
- " typename decltype({})::Traits::RangeType err(0.0);".format(dsa),
- " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
- " {} += err;".format(accum_error),
- " std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
- "}",
+ return ["// L2 error squared of difference between numerical",
+ "// solution and the interpolation of exact solution",
+ " Dune::PDELab::integrateGridFunction({}, {}, 10);".format(dsa, accum_error),
]
-def get_treepath(element, index):
- if isinstance(element, (VectorElement, TensorElement)):
- return (index,)
- if isinstance(element, MixedElement):
- pos, rest = element.extract_subelement_component(index)
- offset = sum(element.sub_elements()[i].value_size() for i in range(pos))
- return (pos,) + get_treepath(element.sub_elements()[pos], index - offset)
- else:
- return ()
-
-
-def treepath_to_index(element, treepath, offset=0):
- if len(treepath) == 0:
- return offset
- index = treepath[0]
- offset = offset + sum(element.sub_elements()[i].value_size() for i in range(index))
- subel = element.sub_elements()[index]
- return treepath_to_index(subel, treepath[1:], offset)
-
-
-def accumulate_L2_squared():
- element = get_trial_element()
- if isinstance(element, MixedElement):
- for i in range(element.value_size()):
- _accumulate_L2_squared(get_treepath(element, i))
- else:
- _accumulate_L2_squared(())
-
-
@preamble
def define_accumulated_L2_error(name):
- t = type_range()
- return "{} {}(0.0);".format(t, name)
+ return "double {}(0.0);".format(name)
def name_accumulated_L2_error():
@@ -164,10 +143,8 @@ def compare_L2_squared():
accum_error = name_accumulated_L2_error()
fail = name_test_fail_variable()
- return ["using std::abs;",
- "using std::isnan;",
- "std::cout << \"\\nl2errorsquared: \" << {} << std::endl << std::endl;".format(accum_error),
- "if (isnan({0}) or abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
+ 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)]
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
index 8509b30f..b439ded1 100644
--- a/python/dune/perftool/pdelab/driver/gridfunctionspace.py
+++ b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
@@ -30,18 +30,14 @@ def type_domainfield():
return "DF"
-def basetype_range():
+@preamble
+def typedef_range(name):
if get_option('opcounter'):
from dune.perftool.pdelab.driver.timings import setup_timer
setup_timer()
- return "oc::OpCounter<double>"
+ return "using {} = oc::OpCounter<double>;".format(name)
else:
- return "double"
-
-
-@preamble
-def typedef_range(name):
- return "using {} = {};".format(name, basetype_range())
+ return "using {} = double;".format(name)
def type_range():
@@ -124,9 +120,6 @@ def typedef_fem(element, name):
df = type_domainfield()
r = type_range()
dim = get_dimension()
- if get_option("blockstructured"):
- include_file("dune/perftool/blockstructured/blockstructuredqkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
if isQk(element):
include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
@@ -184,7 +177,7 @@ def name_test_gfs():
def name_gfs(element, is_dirichlet):
if isinstance(element, (VectorElement, TensorElement)):
- subel = element.sub_elements()[0]
+ 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)
@@ -196,12 +189,12 @@ def name_gfs(element, is_dirichlet):
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, tuple(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[0] else "",
+ "_dirichlet" if is_dirichlet else "",
)
define_gfs(element, is_dirichlet, name)
return name
@@ -221,7 +214,7 @@ def type_trial_gfs():
def type_gfs(element, is_dirichlet):
if isinstance(element, (VectorElement, TensorElement)):
- subel = element.sub_elements()[0]
+ 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)
@@ -233,12 +226,12 @@ def type_gfs(element, is_dirichlet):
subgfs.append(type_gfs(subel, is_dirichlet[k:k + subel.value_size()]))
k = k + subel.value_size()
name = "_".join(subgfs)
- typedef_composite_gfs(element, name, tuple(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[0] else "",
+ "_dirichlet" if is_dirichlet else "",
)
typedef_gfs(element, is_dirichlet, name)
return name
@@ -256,7 +249,6 @@ def define_gfs(element, is_dirichlet, name):
def define_power_gfs(element, is_dirichlet, name, subgfs):
gv = name_leafview()
fem = name_fem(element.sub_elements()[0])
- gfstype = type_gfs(element, is_dirichlet)
return "{} {}({}, {});".format(gfstype, name, gv, fem)
@@ -271,7 +263,7 @@ def typedef_gfs(element, is_dirichlet, name):
vb = type_vectorbackend(element)
gv = type_leafview()
fem = type_fem(element)
- cass = type_constraintsassembler(bool(is_dirichlet[0]))
+ cass = type_constraintsassembler(bool(is_dirichlet))
return "using {} = Dune::PDELab::GridFunctionSpace<{}, {}, {}, {}>;".format(name, gv, fem, cass, vb)
@@ -351,30 +343,3 @@ def type_constraintsassembler(is_dirichlet):
name = "NoConstraintsAssembler"
typedef_no_constraintsassembler(name)
return name
-
-
-def name_trial_subgfs(treepath):
- if len(treepath) == 0:
- return name_trial_gfs()
- else:
- return name_subgfs(treepath)
-
-
-def name_subgfs(treepath):
- gfs = name_trial_gfs()
- name = "{}_{}".format(gfs, "_".join(str(t) for t in treepath))
- define_subgfs(name, treepath)
- return name
-
-
-@preamble
-def define_subgfs(name, treepath):
- t = type_subgfs(treepath)
- gfs = name_trial_gfs()
- return "{} {}({});".format(t, name, gfs)
-
-
-def type_subgfs(tree_path):
- include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
- gfs = type_trial_gfs()
- return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(str(t) for t in tree_path))
diff --git a/python/dune/perftool/pdelab/driver/gridoperator.py b/python/dune/perftool/pdelab/driver/gridoperator.py
index b701ceca..500d8bc7 100644
--- a/python/dune/perftool/pdelab/driver/gridoperator.py
+++ b/python/dune/perftool/pdelab/driver/gridoperator.py
@@ -104,12 +104,6 @@ def name_localoperator(formdata):
return name
-@preamble
-def update_gfs(gfs):
- return ["// Update the GFS to trigger GFS initialization",
- "{}.update();".format(gfs)]
-
-
@preamble
def define_dofestimate(name):
# Provide a worstcase estimate for the number of entries per row based
@@ -118,14 +112,13 @@ def define_dofestimate(name):
geo_factor = 2**get_dimension()
else:
if get_dimension() < 3:
- geo_factor = 3 * get_dimension()
+ 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()
- update_gfs(gfs)
return ["int generic_dof_estimate = {} * {}.maxLocalSize();".format(geo_factor, gfs),
"int {} = {}.get<int>(\"istl.number_of_nnz\", generic_dof_estimate);".format(name, ini)]
diff --git a/python/dune/perftool/pdelab/driver/instationary.py b/python/dune/perftool/pdelab/driver/instationary.py
index 9ab49c9c..bf8759e7 100644
--- a/python/dune/perftool/pdelab/driver/instationary.py
+++ b/python/dune/perftool/pdelab/driver/instationary.py
@@ -1,35 +1,11 @@
-from dune.perftool.generation import (include_file,
- preamble,
- )
+from dune.perftool.generation import preamble
from dune.perftool.pdelab.driver import (get_formdata,
- get_mass_formdata,
- get_trial_element,
is_linear,
- name_initree,
- preprocess_leaf_data,
)
-from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
- type_range,
- )
-from dune.perftool.pdelab.driver.gridoperator import (name_gridoperator,
- name_parameters,
- type_gridoperator,)
-from dune.perftool.pdelab.driver.constraints import (name_bctype_function,
- name_constraintscontainer,
- )
-from dune.perftool.pdelab.driver.interpolate import (interpolate_dirichlet_data,
- name_boundary_function,
- )
+from dune.perftool.pdelab.driver.gridoperator import (type_gridoperator,)
from dune.perftool.pdelab.driver.solve import (print_matrix,
print_residual,
- name_stationarynonlinearproblemsolver,
- name_vector,
- type_stationarynonlinearproblemssolver,
- type_vector,
)
-from dune.perftool.pdelab.driver.vtk import (name_vtk_sequence_writer,
- visualize_initial_condition,
- )
from dune.perftool.options import get_option
@@ -52,6 +28,12 @@ def solve_instationary():
print_residual()
print_matrix()
+ from dune.perftool.pdelab.driver.error import compare_dofs, compare_L2_squared
+ if get_option("exact_solution_expression"):
+ if get_option("compare_dofs"):
+ compare_dofs()
+ if get_option("compare_l2errorsquared"):
+ compare_L2_squared()
@preamble
@@ -60,14 +42,12 @@ def time_loop():
formdata = get_formdata()
params = name_parameters(formdata)
time = name_time()
- element = get_trial_element()
- is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
- bctype = name_bctype_function(element, is_dirichlet)
- gfs = name_gfs(element, is_dirichlet)
- cc = name_constraintscontainer()
+ 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)
- interpolate_dirichlet_data(vector)
# Choose between explicit and implicit time stepping
explicit = get_option('explicit_time_stepping')
@@ -75,8 +55,7 @@ def time_loop():
osm = name_explicitonestepmethod()
apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
else:
- dirichlet = preprocess_leaf_data(element, "dirichlet_expression")
- boundary = name_boundary_function(element, dirichlet)
+ boundary = name_boundary_function(expr)
osm = name_onestepmethod()
apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
@@ -116,6 +95,7 @@ def name_time():
return "time"
+
@preamble
def typedef_timesteppingmethod(name):
r_type = type_range()
@@ -149,8 +129,8 @@ def name_timesteppingmethod():
@preamble
def typedef_instationarygridoperator(name):
include_file("dune/pdelab/gridoperator/onestep.hh", filetag="driver")
- go_type = type_gridoperator(get_formdata())
- mass_go_type = type_gridoperator(get_mass_formdata())
+ 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)
@@ -166,8 +146,8 @@ def type_instationarygridoperator():
@preamble
def define_instationarygridoperator(name):
igo_type = type_instationarygridoperator()
- go = name_gridoperator(get_formdata())
- mass_go = name_gridoperator(get_mass_formdata())
+ go = name_gridoperator(_driver_data['formdata'])
+ mass_go = name_gridoperator(_driver_data['mass_formdata'])
return "{} {}({}, {});".format(igo_type, name, go, mass_go)
@@ -181,7 +161,7 @@ def typedef_onestepmethod(name):
r_type = type_range()
igo_type = type_instationarygridoperator()
snp_type = type_stationarynonlinearproblemssolver(igo_type)
- vector_type = type_vector(get_formdata())
+ vector_type = type_vector(_driver_data['formdata'])
return "using {} = Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}>;".format(name, r_type, igo_type, snp_type, vector_type, vector_type)
@@ -210,7 +190,7 @@ def typedef_explicitonestepmethod(name):
r_type = type_range()
igo_type = type_instationarygridoperator()
ls_type = type_linearsolver()
- vector_type = type_vector(get_formdata())
+ vector_type = type_vector(_driver_data['formdata'])
return "using {} = Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}>;".format(name, r_type, igo_type, ls_type, vector_type)
diff --git a/python/dune/perftool/pdelab/driver/interpolate.py b/python/dune/perftool/pdelab/driver/interpolate.py
index 763283c1..424388a0 100644
--- a/python/dune/perftool/pdelab/driver/interpolate.py
+++ b/python/dune/perftool/pdelab/driver/interpolate.py
@@ -43,21 +43,20 @@ def interpolate_vector(func, gfs, name):
)
-@cached
-def name_boundary_function(element, func):
+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, func[k:k + subel.value_size()]))
+ childs.append(name_boundary_function(subel, dirichlet[k:k + subel.value_size()]))
k = k + subel.value_size()
name = "_".join(childs)
- define_compositegfs_parameterfunction(name, tuple(childs))
+ define_composite_boundary_function(name, childs)
return name
else:
assert isinstance(element, FiniteElement)
- name = get_counted_variable("func")
- define_boundary_function(name, func[0])
+ name = "{}_boundary".format(FEM_name_mangling(element).lower())
+ define_boundary_function(name, dirichlet[0])
return name
@@ -95,13 +94,16 @@ def name_boundary_lambda(boundary):
@preamble
def define_boundary_lambda(name, boundary):
+ from dune.perftool.ufl.execution import Expression
from ufl.classes import Expr
if boundary is None:
- boundary = 0.0
- if isinstance(boundary, (int, float)):
- return "auto {} = [&](const auto& x){{ return {}; }};".format(name, boundary)
- else:
- assert isinstance(boundary, Expr)
+ 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
diff --git a/python/dune/perftool/pdelab/driver/solve.py b/python/dune/perftool/pdelab/driver/solve.py
index 59ef32eb..5e602827 100644
--- a/python/dune/perftool/pdelab/driver/solve.py
+++ b/python/dune/perftool/pdelab/driver/solve.py
@@ -40,14 +40,21 @@ def dune_solve():
include_file("dune/perftool/matrixfree.hh", filetag="driver")
solve = "solveNonlinearMatrixFree({},{});".format(go, x)
elif not linear and not matrix_free:
- go_type = type_gridoperator(get_formdata())
- go = name_gridoperator(get_formdata())
+ 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()
@@ -101,6 +108,7 @@ def define_vector(name, formdata):
return ["{} {}({});".format(vtype, name, gfs), "{} = 0.0;".format(name)]
+
@preamble
def typedef_linearsolver(name):
include_file("dune/pdelab/backend/istl.hh", filetag="driver")
@@ -171,7 +179,7 @@ def name_stationarylinearproblemsolver():
def typedef_stationarynonlinearproblemsolver(name, go_type):
include_file("dune/pdelab/newton/newton.hh", filetag="driver")
ls_type = type_linearsolver()
- x_type = type_vector(get_formdata())
+ x_type = type_vector(_driver_data['formdata'])
return "using {} = Dune::PDELab::Newton<{}, {}, {}>;".format(name, go_type, ls_type, x_type)
diff --git a/python/dune/perftool/pdelab/driver/vtk.py b/python/dune/perftool/pdelab/driver/vtk.py
index 276570de..2f68140e 100644
--- a/python/dune/perftool/pdelab/driver/vtk.py
+++ b/python/dune/perftool/pdelab/driver/vtk.py
@@ -4,10 +4,8 @@ from dune.perftool.generation import (include_file,
from dune.perftool.pdelab.driver import (get_formdata,
get_trial_element,
name_initree,
- preprocess_leaf_data,
)
from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
- name_gfs,
name_trial_gfs,
type_leafview,
)
@@ -67,7 +65,6 @@ def vtkoutput():
include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
vtkwriter = name_vtkwriter()
gfs = name_trial_gfs()
- define_gfs_names(get_trial_element())
vtkfile = name_vtkfile()
predicate = name_predicate()
vec = name_vector(get_formdata())
@@ -92,6 +89,7 @@ def name_predicate():
return "predicate"
+
@preamble
def typedef_vtk_sequence_writer(name):
include_file("dune/grid/io/file/vtk/vtksequencewriter.hh", filetag="driver")
@@ -124,7 +122,7 @@ def visualize_initial_condition():
vtkwriter = name_vtk_sequence_writer()
element = get_trial_element()
define_gfs_names(element)
- gfs = name_trial_gfs()
+ gfs = name_gfs(element)
vector = name_vector(get_formdata())
predicate = name_predicate()
from dune.perftool.pdelab.driver.instationary import name_time
@@ -134,15 +132,12 @@ def visualize_initial_condition():
@preamble
-def _define_gfs_name(element, is_dirichlet, offset):
- from ufl import MixedElement, VectorElement
- if isinstance(element, VectorElement):
- gfs = name_gfs(element, is_dirichlet)
- return "{}.name(\"data_{}to{}\");".format(gfs, offset, offset + element.value_size())
- elif isinstance(element, MixedElement):
+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 + k)
+ define_gfs_name(subel, is_dirichlet[k:k + subel.value_size()], offset=offset + k)
k = k + subel.value_size()
return []
else:
@@ -153,4 +148,4 @@ def _define_gfs_name(element, is_dirichlet, offset):
def define_gfs_names(element):
is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
- return _define_gfs_name(element, is_dirichlet, 0)
+ return _define_gfs_name(element, is_dirichlet)
\ No newline at end of file
--
GitLab