diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py index e4b51c3a15f191f86691b1b0816e1a3057f23073..ff563cfb2bd1a50af135571fc43a248009b26018 100644 --- a/python/dune/perftool/options.py +++ b/python/dune/perftool/options.py @@ -15,8 +15,18 @@ def get_form_compiler_arguments(): parser.add_argument("--numerical-jacobian", action="store_true", help="use numerical jacobians (only makes sense, if uflpdelab for some reason fails to generate analytic jacobians)") parser.add_argument( "--exact-solution-expression", - nargs = 2, - help="requires two arguments: the name of the expression corresponding to solution in the ufl file and the maximal allowed error (between difference of the dof vectors of num. solution and interpolated exact solution") + type=str, + help="name of the exact solution expression in the ufl file") + parser.add_argument( + "--compare-dofs", + type=str, + help="maximal allowed maximum error of difference between degrees of freedom vectors of numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)" + ) + parser.add_argument( + "--compare-l2errorsquared", + type=str, + help="maximal allowed l2 error squared of difference between numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)" + ) parser.add_argument("--interactive", action="store_true", help="whether the optimization process should be guided interactively (also useful for debugging)") # These are the options that I deemed necessary in uflpdelab diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py index 12c6b67486d5c1be803da21aa2b9d1651415e1a6..f246e0d95b68235eb7c467ce01dbc3fae3e4c2a6 100644 --- a/python/dune/perftool/pdelab/driver.py +++ b/python/dune/perftool/pdelab/driver.py @@ -749,7 +749,7 @@ def name_boundary_function(expr): def name_solution_function(expr): # Get expression representing solution from dune.perftool.generation import get_global_context_value - solution_expression_name = get_option("exact_solution_expression")[0] + solution_expression_name = get_option("exact_solution_expression") solution_expression = get_global_context_value("data").object_by_name[solution_expression_name] # Create function just like it is done for boundary_function @@ -924,19 +924,55 @@ def name_vtkfile(): @preamble -def exact_solution(): +def compare_dofs(): v = name_vector() solution = name_vector_from_solution_expression() - - return ["using Dune::PDELab::Backend::native;", + 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("exact_solution_expression")[1]), + "if (maxerror>{})".format(get_option("compare_dofs")), + " throw;"] + + +@preamble +def define_discrete_grid_function(gfs, vector_name, dgf_name): + return ["typedef Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})> DGF;".format(gfs, vector_name), + "DGF {}({},{});".format(dgf_name, gfs, vector_name)] + + +@symbol +def name_discrete_grid_function(gfs, vector_name): + dgf_name = vector_name + "_dgf" + define_discrete_grid_function(gfs, vector_name, dgf_name) + return dgf_name + + +@preamble +def compare_L2_squared(): + element = _form.coefficients()[0].ufl_element() + v = name_vector() + gfs = name_gfs(element) + vdgf = name_discrete_grid_function(gfs, v) + solution_function = name_solution_function(element) + + 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", + "typedef DifferenceSquaredAdapter<decltype({}),decltype({})> DifferenceSquared;".format(solution_function, vdgf), + "DifferenceSquared differencesquared({},{});".format(solution_function, vdgf), + "typename DifferenceSquared::Traits::RangeType l2errorsquared(0.0);", + "Dune::PDELab::integrateGridFunction(differencesquared,l2errorsquared,10);", + "std::cout << \"l2errorsquared: \" << l2errorsquared << std::endl;", + "if (l2errorsquared>{})".format(get_option("compare_l2errorsquared")), " throw;"] @@ -1024,7 +1060,12 @@ def vtkoutput(): from ufl import MixedElement, VectorElement, TensorElement if isinstance(_form.coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)): raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.") - exact_solution() + + 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)] diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt index 5ebedaa84ecbcecfa73d5c17ed4d24597aa7038c..a78cccd05c560d4ed69edf860d2cbbc25f3928b4 100644 --- a/test/nonlinear/CMakeLists.txt +++ b/test/nonlinear/CMakeLists.txt @@ -11,7 +11,7 @@ dune_add_system_test(TARGET nonlinear_symdiff add_generated_executable(UFLFILE nonlinear_dg.ufl TARGET nonlinear_dg - FORM_COMPILER_ARGS --exact-solution-expression g 1e-2 + FORM_COMPILER_ARGS --exact-solution-expression g --compare-dofs 1e-1 --compare-l2errorsquared 1e-3 ) dune_add_system_test(TARGET nonlinear_dg diff --git a/test/nonlinear/nonlinear_dg.mini b/test/nonlinear/nonlinear_dg.mini new file mode 100644 index 0000000000000000000000000000000000000000..d2b26e27c4575901a0d094a8da176f2dd13d73d7 --- /dev/null +++ b/test/nonlinear/nonlinear_dg.mini @@ -0,0 +1,11 @@ +__name = nonlinear_dg + +lowerleft = 0.0 0.0 +upperright = 1.0 1.0 +elements = 32 32 +elementType = simplical + +[wrapper.vtkcompare] +name = {__name} +reference = nonlinear_dg +extension = vtu