Time dependent parameters/boundary condition

patches/apply_patches.sh

@@ -10,6 +10,7 @@ popd
 
 pushd python/ufl
 git apply ../../patches/ufl/conditional-uflid.patch
+git apply ../../patches/ufl/0001-Remove-special-case-for-variable-in-ufl2dot.patch
 popd
 
 pushd python/ufl

patches/ufl/0001-Remove-special-case-for-variable-in-ufl2dot.patch

+From 6f2931706f28cd29e3ed72851a7712815a23f474 Mon Sep 17 00:00:00 2001
+From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
+Date: Thu, 9 Nov 2017 14:13:34 +0100
+Subject: [PATCH] Remove special case for variable in ufl2dot
+
+---
+ ufl/formatting/ufl2dot.py | 9 ++-------
+ 1 file changed, 2 insertions(+), 7 deletions(-)
+
+diff --git a/ufl/formatting/ufl2dot.py b/ufl/formatting/ufl2dot.py
+index 5fdca148..e0387a9b 100644
+--- a/ufl/formatting/ufl2dot.py
++++ b/ufl/formatting/ufl2dot.py
+@@ -176,13 +176,8 @@ def build_entities(e, nodes, edges, nodeoffset, prefix="", labeller=None):
+     if labeller is None:
+         labeller = ReprLabeller()
+
+-    # Special-case Variable instances
+-    if isinstance(e, Variable):  # FIXME: Is this really necessary?
+-        ops = (e._expression,)
+-        label = "variable %d" % e._label._count
+-    else:
+-        ops = e.ufl_operands
+-        label = labeller(e)
++    ops = e.ufl_operands
++    label = labeller(e)
+
+     # Create node for parent e
+     nodename = "%sn%04d" % (prefix, len(nodes) + nodeoffset)
+--
+2.11.0

python/dune/perftool/pdelab/argument.py

@@ -9,7 +9,6 @@ from dune.perftool.options import get_option
 from dune.perftool.generation import (domain,
                                       function_mangler,
                                       iname,
-                                      globalarg,
                                       valuearg,
                                       get_global_context_value,
                                       kernel_cached,

python/dune/perftool/pdelab/driver/instationary.py

@@ -79,6 +79,8 @@ def time_loop():
     return ["",
             "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
             "double dt = {}.get<double>(\"instat.dt\", 0.1);".format(ini),
+            "int step_number(0);"
+            "int nth = {}.get<int>(\"instat.nth\", 1);".format(ini),
             "while (time<T-1e-8){",
             "  // Assemble constraints for new time step",
             "  {}.setTime({}+dt);".format(params, time),
@@ -92,8 +94,11 @@ def time_loop():
             "  {} = {}new;".format(vector, vector),
             "  time += dt;",
             "",
-            "  // Output to VTK File",
-            "  {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
+            "  step_number += 1;",
+            "  if (step_number%nth == 0){",
+            "    // Output to VTK File",
+            "    {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
+            "  }",
             "}",
             ""]
 

python/dune/perftool/pdelab/parameter.py

@@ -5,6 +5,7 @@ from dune.perftool.generation import (class_basename,
                                       constructor_parameter,
                                       generator_factory,
                                       get_backend,
+                                      get_global_context_value,
                                       initializer_list,
                                       kernel_cached,
                                       preamble,
@@ -38,8 +39,11 @@ def define_parameterclass(name):
 
 
 def name_paramclass():
-    define_parameterclass("param")
-    return "param"
+    formdata = get_global_context_value("formdata")
+    from dune.perftool.pdelab.driver.gridoperator import name_parameters
+    name = name_parameters(formdata)
+    define_parameterclass(name)
+    return name
 
 
 @class_member(classtag="parameterclass")

python/dune/perftool/ufl/visitor.py

@@ -3,13 +3,18 @@ This module defines the main visitor algorithm transforming ufl expressions
 to pymbolic and loopy.
 """
 from dune.perftool.error import PerftoolUFLError
-from dune.perftool.generation import get_global_context_value, domain
+from dune.perftool.generation import (get_global_context_value,
+                                      domain,
+                                      globalarg,
+                                      valuearg,
+                                      )
 from dune.perftool.ufl.flatoperators import get_operands
 from dune.perftool.ufl.modified_terminals import (ModifiedTerminalTracker,
                                                   Restriction,
                                                   )
 from dune.perftool.tools import maybe_wrap_subscript
 from dune.perftool.options import get_option
+from dune.perftool.pdelab.parameter import name_paramclass, name_time
 from loopy import Reduction
 
 from pymbolic.primitives import (Call,
@@ -27,7 +32,8 @@ from ufl import (VectorElement,
                  TensorElement,
                  TensorProductElement,
                  )
-from ufl.classes import (FixedIndex,
+from ufl.classes import (Coefficient,
+                         FixedIndex,
                          IndexSum,
                          JacobianDeterminant,
                          )
@@ -172,6 +178,21 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             # And return a symbol
             return Variable(name)
 
+    def variable(self, o):
+        # Find Coefficient node below the variable node and check that
+        # is has the id reserved for the time variable
+        coeff = o.ufl_operands[0]
+        while not isinstance(coeff, Coefficient):
+            coeff = coeff.ufl_operands[0]
+        # TODO: Find a better way to create the time Variable in the ufl file
+        assert coeff.count() == 2
+
+        param = name_paramclass()
+        time = name_time()
+        name = param + "." + time
+        valuearg(name, dtype=np.float64)
+        return Variable(name)
+
     #
     # Handlers for all indexing related stuff
     #

test/heatequation/CMakeLists.txt

@@ -7,6 +7,11 @@ dune_add_formcompiler_system_test(UFLFILE heatequation.ufl
                                   INIFILE heatequation.mini
                                   )
 
+dune_add_formcompiler_system_test(UFLFILE heatequation_time_dependent_bc.ufl
+                                  BASENAME heatequation_time_dependent_bc
+                                  INIFILE heatequation_time_dependent_bc.mini
+                                  )
+
 #=====
 # DG
 #=====

test/heatequation/heatequation_time_dependent_bc.mini

+__name = heatequation_time_dependent_bc_{__exec_suffix}
+__exec_suffix = implicit, explicit | expand scheme
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 16 16
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = heatequation_ref
+extension = vtu
+
+[formcompiler]
+explicit_time_stepping = 0, 1 | expand scheme
+compare_l2errorsquared = 2e-4
+
+[instat]
+T = 1.0
+dt = 1e-1
+nth = 1
+
+
+# Disable explicit tests for now
+{__exec_suffix} == explicit | exclude

test/heatequation/heatequation_time_dependent_bc.ufl

+cell = triangle
+
+x = SpatialCoordinate(cell)
+time = variable(Constant(cell, count=2))
+
+nu = 1.0/10
+
+c = (0.5-x[0])**2 + (0.5-x[1])**2
+g = exp(-nu*time) * exp(-1.*c)
+f = 4*(1.-c)*g - nu*g
+
+V = FiniteElement("CG", "triangle", 1)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+mass = (u*v)*dx
+poisson = (inner(grad(u), grad(v)) - f*v)*dx
+
+forms = [mass, poisson]
+dirichlet_expression = g
+is_dirichlet = 1
+exact_solution = g