Compare with exact solution

python/dune/perftool/pdelab/driver.py

@@ -1387,6 +1387,17 @@ def time_loop():
     elif not linear and not matrix_free:
         assert False
 
+    print_residual()
+    print_matrix()
+    if get_option("exact_solution_expression"):
+        from ufl import MixedElement, VectorElement, TensorElement
+        if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
+            raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
+
+        if get_option("compare_dofs"):
+            compare_dofs()
+        if get_option("compare_l2errorsquared"):
+            compare_L2_squared()
 
 def generate_driver(formdatas, data):
     # The driver module uses a global dictionary for storing necessary data

test/heatequation/CMakeLists.txt

 add_generated_executable(UFLFILE heatequation.ufl
                          TARGET heatequation_implicit
-                         FORM_COMPILER_ARGS
+                         FORM_COMPILER_ARGS --exact-solution-expression g --compare-l2errorsquared 1e-7
                          )
 
 
@@ -10,7 +10,7 @@ dune_add_system_test(TARGET heatequation_implicit
 
 add_generated_executable(UFLFILE heatequation.ufl
                          TARGET heatequation_explicit
-                         FORM_COMPILER_ARGS --explicit-time-stepping
+                         FORM_COMPILER_ARGS --explicit-time-stepping --exact-solution-expression g --compare-l2errorsquared 1e-7
                          )
 
 
@@ -18,14 +18,10 @@ dune_add_system_test(TARGET heatequation_explicit
                      INIFILE heatequation_explicit.mini
                      )
 
+# Hand written reference solutions
 add_executable(heatequation_implicit_ref heatequation_implicit_ref.cc)
 dune_symlink_to_source_files(FILES heatequation_implicit_ref.ini)
 set_target_properties(heatequation_implicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
 
 add_executable(heatequation_explicit_ref heatequation_explicit_ref.cc)
 set_target_properties(heatequation_explicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
-
-# # Add the reference code with the PDELab localoperator that produced
-# # the reference vtk file
-# add_executable(heatequation_ref heatequation_main.cc)
-# set_target_properties(heatequation_ref PROPERTIES EXCLUDE_FROM_ALL 1)

test/heatequation/heatequation.ufl

+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
 g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
 
 V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
@@ -5,6 +6,6 @@ u = TrialFunction(V)
 v = TestFunction(V)
 
 mass = (u*v)*dx
-poisson = (inner(grad(u), grad(v)))*dx
+poisson = (inner(grad(u), grad(v)) - f*v)*dx
 
 forms = [mass, poisson]
\ No newline at end of file

test/heatequation/heatequation_explicit_ref.cc

@@ -53,9 +53,13 @@ public:
 
   //! source term
   typename Traits::RangeFieldType
-  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  f (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
   {
-    return 0.0;
+    typename Traits::DomainType x = e.geometry().global(xlocal);
+    Dune::FieldVector<double,2> c(0.5);
+    c-= x;
+    return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+    // return 0.0;
   }
 
   //! boundary condition type function

test/heatequation/problem.hh

@@ -29,9 +29,9 @@ public:
   template<typename E, typename X>
   Number f (const E& e, const X& local) const
   {
-    return 0.0;
-    // auto x = e.geometry().global(local);
-    // Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+    // return 0.0;
+    auto x = e.geometry().global(local);
+    Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
   }
 
   //! boundary condition type function (true = Dirichlet)