[WIP] Large steps toward solving instationary problems

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()

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 = {}

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"