From ac2c1a15292f832aa3c917e281748aeb1c4c10ec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Tue, 23 Aug 2016 15:51:13 +0200
Subject: [PATCH] [WIP] Large steps toward solving instationary problems
---
python/dune/perftool/pdelab/driver.py | 459 ++++++++++++++-----
python/dune/perftool/pdelab/localoperator.py | 8 +
python/dune/perftool/pdelab/parameter.py | 26 ++
3 files changed, 374 insertions(+), 119 deletions(-)
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 1ae59a3c..9baebbc2 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -13,6 +13,40 @@ from dune.perftool.generation import (generator_factory,
)
from dune.perftool.options import get_option
+# Have a global variable with the entire form data. This allows functions that depend
+# deterministically on the entire data set to directly access it instead of passing it
+# through the entire generator chain.
+_driver_data = {}
+
+
+# Have a function access this global data structure
+def set_driver_data(formdatas, data):
+ assert (len(formdatas) <= 2)
+ if len(formdatas) == 1:
+ _driver_data['form'] = formdatas[0].preprocessed_form
+ _driver_data['formdata'] = formdatas[0]
+ else:
+ mass_index = mass_form_index(formdatas, data)
+ if mass_index is None:
+ raise NotImplementedError("Form for mass matrix needs to have name 'mass' in ufl file.")
+ _driver_data['mass_form'] = formdatas[mass_index].preprocessed_form
+ _driver_data['mass_formdata'] = formdatas[mass_index]
+ _driver_data['form'] = formdatas[1 - mass_index].preprocessed_form
+ _driver_data['formdata'] = formdatas[1 - mass_index]
+
+ _driver_data['data'] = data
+
+
+def is_stationary():
+ return 'mass_form' not in _driver_data
+
+
+def form_name_suffix(name, formdata):
+ from dune.perftool.pdelab.localoperator import name_form
+ data = _driver_data['data']
+ form_name = name_form(formdata, data)
+ return name + '_' + form_name
+
def has_constraints(element):
from ufl import MixedElement, VectorElement
@@ -47,26 +81,6 @@ def mass_form_index(formdatas, data):
except:
continue
-# Have a global variable with the entire form data. This allows functions that depend
-# deterministically on the entire data set to directly access it instead of passing it
-# through the entire generator chain.
-_driver_data = {}
-
-
-# Have a function access this global data structure
-def set_driver_data(formdatas, data):
- assert (len(formdatas) <= 2)
- if len(formdatas) == 1:
- _driver_data['form'] = formdatas[0].preprocessed_form
- _driver_data['formdata'] = formdatas[0]
- else:
- mass_index = mass_form_index(formdatas, data)
- if mass_index == None:
- raise NotImplementedError("Form for mass matrix needs to have name 'mass' in ufl file.")
- _driver_data['mass_form'] = formdatas[mass_index].preprocessed_form
- _driver_data['form'] = formdatas[1-mass_index].preprocessed_form
- _driver_data['formdata'] = formdatas[1-mass_index]
-
def is_linear(form):
'''Test if form is linear in test function'''
@@ -596,29 +610,29 @@ def name_matrixbackend():
@symbol
-def type_parameters():
+def type_parameters(formdata):
from dune.perftool.generation import get_global_context_value
data = get_global_context_value("data")
- formdata = _driver_data['formdata']
from dune.perftool.pdelab.parameter import parameterclass_basename
name = parameterclass_basename(formdata, data)
return name
@preamble
-def define_parameters(name):
- partype = type_parameters()
+def define_parameters(name, formdata):
+ partype = type_parameters(formdata)
return "{} {};".format(partype, name)
@symbol
-def name_parameters():
- define_parameters("params")
- return "params"
+def name_parameters(formdata):
+ name = form_name_suffix("params", formdata)
+ define_parameters(name, formdata)
+ return name
@preamble
-def typedef_localoperator(name):
+def typedef_localoperator(name, formdata):
# No Parameter class here, yet
# params = type_parameters()
# return "typedef LocalOperator<{}> {};".format(params, name)
@@ -626,7 +640,6 @@ def typedef_localoperator(name):
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")
- formdata = _driver_data['formdata']
filename = get_option("operator_file")
include_file(filename, filetag="driver")
from dune.perftool.pdelab.localoperator import localoperator_basename
@@ -635,30 +648,32 @@ def typedef_localoperator(name):
@symbol
-def type_localoperator():
- typedef_localoperator("LOP")
- return "LOP"
+def type_localoperator(formdata):
+ name = form_name_suffix("LOP", formdata).upper()
+ typedef_localoperator(name, formdata)
+ return name
@preamble
-def define_localoperator(name):
- loptype = type_localoperator()
+def define_localoperator(name, formdata):
+ loptype = type_localoperator(formdata)
ini = name_initree()
- params = name_parameters()
+ params = name_parameters(formdata)
return "{} {}({}, {});".format(loptype, name, ini, params)
@symbol
-def name_localoperator():
- define_localoperator("lop")
- return "lop"
+def name_localoperator(formdata):
+ name = form_name_suffix("lop", formdata)
+ define_localoperator(name, formdata)
+ return name
@preamble
-def typedef_gridoperator(name):
+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()
+ 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()
@@ -669,19 +684,21 @@ def typedef_gridoperator(name):
@symbol
-def type_gridoperator():
- typedef_gridoperator("GO")
- return "GO"
+def type_gridoperator(formdata):
+ name = form_name_suffix("GO", formdata).upper()
+ typedef_gridoperator(name, formdata)
+ return name
@preamble
-def define_gridoperator(name):
- gotype = type_gridoperator()
+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()
+ 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),
@@ -689,26 +706,28 @@ def define_gridoperator(name):
@symbol
-def name_gridoperator():
- define_gridoperator("go")
- return "go"
+def name_gridoperator(formdata):
+ name = form_name_suffix("go", formdata)
+ define_gridoperator(name, formdata)
+ return name
@preamble
-def typedef_vector(name):
- gotype = type_gridoperator()
- return "typedef {}::Traits::Domain {};".format(gotype, name)
+def typedef_vector(name, formdata):
+ go_type = type_gridoperator(formdata)
+ return "typedef {}::Traits::Domain {};".format(go_type, name)
@symbol
-def type_vector():
- typedef_vector("V")
- return "V"
+def type_vector(formdata):
+ name = form_name_suffix("V", formdata).upper()
+ typedef_vector(name, formdata)
+ return name
@preamble
-def define_vector(name):
- vtype = type_vector()
+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)]
@@ -734,10 +753,17 @@ def define_boundary_function(boundary, name):
gv = name_leafview()
lambdaname = name_boundary_lambda(boundary, name)
include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
- return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
- gv,
- lambdaname,
- )
+ 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
@@ -787,8 +813,8 @@ def name_solution_function(expr):
@preamble
-def interpolate_vector(name):
- define_vector(name)
+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)
@@ -810,18 +836,19 @@ def interpolate_solution_expression(name):
)
-def maybe_interpolate_vector(name):
+def maybe_interpolate_vector(name, formdata):
element = _driver_data['form'].coefficients()[0].ufl_element()
if has_constraints(element):
- interpolate_vector(name)
+ interpolate_vector(name, formdata)
else:
- define_vector(name)
+ define_vector(name, formdata)
@symbol
-def name_vector():
- maybe_interpolate_vector("x")
- return "x"
+def name_vector(formdata):
+ name = form_name_suffix("x", formdata)
+ maybe_interpolate_vector(name, formdata)
+ return name
@symbol
@@ -869,82 +896,274 @@ def name_reduction():
@preamble
def typedef_stationarylinearproblemsolver(name):
include_file("dune/pdelab/stationary/linearproblem.hh", filetag="driver")
- gotype = type_gridoperator()
+ gotype = type_gridoperator(_driver_data['formdata'])
lstype = type_linearsolver()
- xtype = type_vector()
+ xtype = type_vector(_driver_data['formdata'])
return "typedef Dune::PDELab::StationaryLinearProblemSolver<{}, {}, {}> {};".format(gotype, lstype, xtype, name)
-@preamble
-def typedef_stationarynonlinearproblemsolver(name):
- include_file("dune/pdelab/newton/newton.hh", filetag="driver")
- gotype = type_gridoperator()
- lstype = type_linearsolver()
- xtype = type_vector()
- return "typedef Dune::PDELab::Newton<{}, {}, {}> {};".format(gotype, lstype, xtype, name)
-
-
@symbol
def type_stationarylinearproblemsolver():
typedef_stationarylinearproblemsolver("SLP")
return "SLP"
-@symbol
-def type_stationarynonlinearproblemssolver():
- typedef_stationarynonlinearproblemsolver("SNP")
- return "SNP"
-
-
@preamble
def define_stationarylinearproblemsolver(name):
slptype = type_stationarylinearproblemsolver()
- go = name_gridoperator()
+ formdata = _driver_data['formdata']
+ go = name_gridoperator(formdata)
ls = name_linearsolver()
- x = name_vector()
+ x = name_vector(formdata)
red = name_reduction()
return "{} {}({}, {}, {}, {});".format(slptype, name, go, ls, x, red)
-@preamble
-def define_stationarynonlinearproblemsolver(name):
- snptype = type_stationarynonlinearproblemssolver()
- go = name_gridoperator()
- x = name_vector()
- ls = name_linearsolver()
- return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
-
-
@symbol
def name_stationarylinearproblemsolver():
define_stationarylinearproblemsolver("slp")
return "slp"
+@preamble
+def typedef_stationarynonlinearproblemsolver(name):
+ include_file("dune/pdelab/newton/newton.hh", filetag="driver")
+ go_type = type_gridoperator(_driver_data['formdata'])
+ ls_type = type_linearsolver()
+ x_type = type_vector(_driver_data['formdata'])
+ return "typedef Dune::PDELab::Newton<{}, {}, {}> {};".format(go_type, ls_type, x_type, name)
+
+
+@symbol
+def type_stationarynonlinearproblemssolver():
+ typedef_stationarynonlinearproblemsolver("SNP")
+ return "SNP"
+
+
+@preamble
+def define_stationarynonlinearproblemsolver(name):
+ snptype = type_stationarynonlinearproblemssolver()
+ go = name_gridoperator(_driver_data['formdata'])
+ x = name_vector(_driver_data['formdata'])
+ ls = name_linearsolver()
+ return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
+
+
@symbol
def name_stationarynonlinearproblemsolver():
define_stationarynonlinearproblemsolver("snp")
return "snp"
+@preamble
+def typedef_timesteppingmethod(name):
+ r_type = type_range()
+ return "typedef Dune::PDELab::OneStepThetaParameter<{}> {};".format(r_type, name)
+
+
+@symbol
+def type_timesteppingmethod():
+ typedef_timesteppingmethod("TSM")
+ return "TSM"
+
+
+@preamble
+def define_timesteppingmethod(name, explicit):
+ tsm_type = type_timesteppingmethod()
+ # TODO enable theta from ini file
+ if explicit:
+ return "{} {}(0.0);".format(tsm_type, name)
+ else:
+ return "{} {}(1.0);".format(tsm_type, name)
+
+
+@symbol
+def name_timesteppingmethod(explicit):
+ define_timesteppingmethod("tsm", explicit)
+ 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'])
+ return "typedef Dune::PDELab::OneStepGridOperator<{},{}> {};".format(go_type, mass_go_type, name)
+
+
+@symbol
+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)
+
+
+@symbol
+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()
+ vector_type = type_vector(_driver_data['formdata'])
+ return "typedef Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}> {};".format(r_type, igo_type, snp_type, vector_type, vector_type, name)
+
+
+@symbol
+def type_onestepmethod():
+ typedef_onestepmethod("OSM")
+ return "OSM"
+
+
+@preamble
+def define_onestepmethod(name):
+ ilptype = type_onestepmethod()
+ explicit = False
+ tsm = name_timesteppingmethod(explicit)
+ igo = name_instationarygridoperator()
+ snp = name_stationarynonlinearproblemsolver()
+ return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)
+
+
+@symbol
+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 "typedef Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}> {};".format(r_type, igo_type, ls_type, vector_type, name)
+
+
+@symbol
+def type_explicitonestepmethod():
+ typedef_explicitonestepmethod("EOSM")
+ return "EOSM"
+
+
+@preamble
+def define_explicitonestepmethod(name):
+ eosm_type = type_explicitonestepmethod()
+ explicit = True
+ tsm = name_timesteppingmethod(explicit)
+ igo = name_instationarygridoperator()
+ ls = name_linearsolver()
+ return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)
+
+
+@symbol
+def name_explicitonestepmethod():
+ define_explicitonestepmethod("eosm")
+ return "eosm"
+
+
+def time_loop(osm):
+ ini = name_initree()
+ params = name_parameters()
+ formdata = _driver_data['formdata']
+ vector_type = type_vector(formdata)
+ vector_name = name_vector(formdata)
+ expr = _driver_data['form'].coefficients()[0].ufl_element()
+ boundary_name = name_boundary_function(expr)
+ bctype_name = name_bctype_function(expr)
+ gfs_name = name_gfs(expr)
+ cc_name = name_constraintscontainer(expr)
+
+ return ["",
+ "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
+ "double dt = {}.get<double>(\"instat.dt\", 0.01);".format(ini),
+ "double time = 0;",
+ "",
+ "while (time<T-1e-8){",
+ " {}.setTime(time+dt);".format(params),
+ " Dune::PDELab::constraints({}, {}, {});".format(bctype_name, gfs_name, cc_name),
+ "",
+ " {} {}new({});".format(vector_type, vector_name, vector_name),
+ " {}.apply(time,dt, {}, {}, {}new);".format(osm, vector_name, boundary_name, vector_name),
+ "",
+ " {} = {}new;".format(vector_name, vector_name),
+ " time += dt;",
+ "}",
+ ""]
+
+
+def explicit_time_loop(eosm):
+ ini = name_initree()
+ formdata = _driver_data['formdata']
+ params = name_parameters(formdata)
+ vector_type = type_vector(formdata)
+ vector_name = name_vector(formdata)
+ expr = _driver_data['form'].coefficients()[0].ufl_element()
+ bctype_name = name_bctype_function(expr)
+ gfs_name = name_gfs(expr)
+ cc_name = name_constraintscontainer(expr)
+
+ return ["",
+ "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
+ "double dt = {}.get<double>(\"instat.dt\", 0.01);".format(ini),
+ "double time = 0;",
+ "",
+ "while (time<T-1e-8){",
+ " {}.setTime(time+dt);".format(params),
+ " Dune::PDELab::constraints({}, {}, {});".format(bctype_name, gfs_name, cc_name),
+ "",
+ " {} {}new({});".format(vector_type, vector_name, vector_name),
+ " {}.apply(time,dt, {}, {}new);".format(eosm, vector_name, vector_name),
+ "",
+ " {} = {}new;".format(vector_name, vector_name),
+ " time += dt;",
+ "}",
+ ""]
+
+
@preamble
def dune_solve():
# Test if form is linear in ansatzfunction
- if is_linear(_driver_data['form']):
- if get_option("matrix_free"):
- go = name_gridoperator()
- x = name_vector()
- include_file("dune/perftool/matrixfree.hh", filetag="driver")
- return "solveMatrixFree({},{});".format(go, x)
- else:
- slp = name_stationarylinearproblemsolver()
- return "{}.apply();".format(slp)
+ linear = is_linear(_driver_data['form'])
+
+ # Test if problem is stationary
+ stationary = is_stationary()
+
+ # Test wether we want to do matrix free operator evaluation
+ matrix_free = get_option('matrix_free')
+
+ if linear and stationary and matrix_free:
+ go = name_gridoperator(_driver_data['formdata'])
+ x = name_vector()
+ include_file("dune/perftool/matrixfree.hh", filetag="driver")
+ return "solveMatrixFree({},{});".format(go, x)
+ elif linear and stationary and not matrix_free:
+ slp = name_stationarylinearproblemsolver()
+ return "{}.apply();".format(slp)
+ elif linear and not stationary and not matrix_free:
+ # TODO -> form compiler argument for switching explicit/implicit
+ # osm = name_onestepmethod()
+ # return time_loop(osm)
+ eosm = name_explicitonestepmethod()
+ return explicit_time_loop(eosm)
+ elif not linear and stationary and matrix_free:
+ raise NotImplementedError("Nonlinear and matrix free is not yet implemented")
+ elif not linear and stationary and not matrix_free:
+ snp = name_stationarynonlinearproblemsolver()
+ return "{}.apply();".format(snp)
else:
- if get_option("matrix_free"):
- raise NotImplementedError("Nonlinear and matrix free is not yet implemented")
- else:
- snp = name_stationarynonlinearproblemsolver()
- return "{}.apply();".format(snp)
+ assert False
@preamble
@@ -1016,9 +1235,10 @@ def compare_L2_squared():
@preamble
def print_residual():
ini = name_initree()
- n_go = name_gridoperator()
- v = name_vector()
- t_v = type_vector()
+ formdata = _driver_data['formdata']
+ n_go = name_gridoperator(formdata)
+ v = name_vector(formdata)
+ t_v = type_vector(formdata)
return ["if ({}.get<bool>(\"printresidual\", false)) {{".format(ini),
" using Dune::PDELab::Backend::native;",
@@ -1031,11 +1251,12 @@ def print_residual():
@preamble
def print_matrix():
+ formdata = _driver_data['formdata']
ini = name_initree()
- t_go = type_gridoperator()
- n_go = name_gridoperator()
- v = name_vector()
- t_v = type_vector()
+ 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),
" typedef typename {}::Traits::Jacobian M;".format(t_go),
@@ -1086,7 +1307,7 @@ def vtkoutput():
include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
vtkwriter = name_vtkwriter()
gfs = name_gfs(element)
- vec = name_vector()
+ vec = name_vector(_driver_data['formdata'])
vtkfile = name_vtkfile()
predicate = name_predicate()
dune_solve()
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index b4990c0f..965efb39 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -339,6 +339,14 @@ def generate_localoperator_kernels(formdata, data):
name_paramclass()
base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
+ from dune.perftool.pdelab.driver import is_stationary
+ if not is_stationary():
+ # TODO replace double with clever typename stuff
+ base_class('Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>', classtag="operator")
+
+ # Create set time method in parameter class
+ from dune.perftool.pdelab.parameter import define_set_time_method
+ define_set_time_method()
# Have a data structure collect the generated kernels
operator_kernels = {}
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index cda24945..05c21e82 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -43,6 +43,32 @@ def name_paramclass():
return "param"
+@class_member(classtag="parameterclass", access=AccessModifier.PRIVATE)
+def define_time(name):
+ initializer_list(name, ["0.0"], classtag="parameterclass")
+ return "double {};".format(name)
+
+
+@symbol
+def name_time():
+ define_time("t")
+ return "t"
+
+
+@class_member("parameterclass", access=AccessModifier.PUBLIC)
+def define_set_time_method():
+ time_name = name_time()
+ # TODO double?
+ result = ["// Set time in instationary case",
+ "void setTime (double t_)",
+ "{",
+ " {} = t_;".format(time_name),
+ "}"
+ ]
+
+ return result
+
+
@class_member("parameterclass", access=AccessModifier.PUBLIC)
def define_parameter_function_class_member(name, expr, t, cell):
geot = "E" if cell else "I"
--
GitLab