Add better LaplacianJumpIndicator tests.

These new tests are based on the Biharmonic equation (and its
transient variant) so the Laplacian jump indicator is appropriate for
estimating the error.

In the process of developing these tests, it was discovered that MOOSE
is not currently imposing Dirichlet BCs for Hermite elements
correctly. Specifically, MOOSE constrains only the *value* degrees of
freedom when Dirichlet BCs are imposed, but setting BCs for Hermite
elements involves the gradient degrees of freedom as well, according
to @roystgnr. The current workaround is to use a penalty approach for
enforcing the boundary conditions, but this comes with the obvious
drawbacks of dramatically increasing the initial residual and changing
the conditioning of the system.

Refs #2190.

test/tests/indicators/laplacian_jump_indicator/biharmonic.i

@@ -0,0 +1,123 @@
+[GlobalParams]
+  # Parameters used by Functions.
+  vars   = 'c'
+  vals   = '50'
+[]
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = -.5
+  xmax = .5
+  ymin = -.5
+  ymax = .5
+  nx = 10
+  ny = 10
+[]
+
+[Variables]
+  [./u]
+    order = THIRD
+    family = HERMITE
+  [../]
+[]
+
+[Kernels]
+  [./biharmonic]
+    type = Biharmonic
+    variable = u
+  [../]
+  [./body_force]
+    type = BodyForce
+    variable = u
+    function = forcing_func
+  [../]
+[]
+
+[BCs]
+  active = 'all_value all_flux'
+  [./all_value]
+    type = FunctionPenaltyDirichletBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+  [./all_flux]
+    type = FunctionPenaltyFluxBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+  [./all_laplacian]
+    type = BiharmonicLapBC
+    variable = u
+    boundary = 'left right top bottom'
+    laplacian_function = lapu_func
+  [../]
+[]
+
+[Adaptivity]
+  [./Indicators]
+    [./error]
+      type = LaplacianJumpIndicator
+      variable = u
+      scale_by_flux_faces = true
+    [../]
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+
+  # Note: the unusually tight tolerances here are due to the penalty
+  # BCs (currently the only way of accurately Dirichlet boundary
+  # conditions on Hermite elements in MOOSE).
+  nl_rel_tol = 1.e-15
+  l_tol = 1.e-15
+
+  # We have exact Jacobians
+  solve_type = 'NEWTON'
+
+  # Use 6x6 quadrature to ensure the forcing function is integrated
+  # accurately.
+  [./Quadrature]
+    type = GAUSS
+    order = ELEVENTH
+  [../]
+[]
+
+[Functions]
+  [./u_func]
+    type   = ParsedGradFunction
+    value  = 'exp(-c*(x^2+y^2))'
+    grad_x = '-2*c*exp(-c*(x^2+y^2))*x'
+    grad_y = '-2*c*exp(-c*(x^2+y^2))*y'
+  [../]
+  [./lapu_func]
+    type   = ParsedFunction
+    value  = '4*c*(c*(x^2+y^2) - 1)*exp(-c*(x^2+y^2))'
+  [../]
+  [./forcing_func]
+    type   = ParsedFunction
+    value  = '16*c^2*(c^2*(x^2+y^2)^2 - 4*c*(x^2+y^2) + 2)*exp(-c*(x^2+y^2))'
+  [../]
+[]
+
+[Postprocessors]
+  [./l2_error]
+    type = ElementL2Error
+    variable = u
+    function = u_func
+  [../]
+  [./h1_error]
+    type = ElementH1Error
+    variable = u
+    function = u_func
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/indicators/laplacian_jump_indicator/biharmonic.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+from sympy import *
+from fractions import Fraction
+
+r, theta, t, c = sympify('r, theta, t, c')
+
+# A simple Gaussian function (all derivatives wrt theta are 0).
+# This script should print the following:
+#
+# Steady version: f = du/dt + Biharmonic(u)
+# lapu = 4*c*(c*r**2 - 1)*exp(-c*r**2)
+# bihu = 16*c**2*(c**2*r**4 - 4*c*r**2 + 2)*exp(-c*r**2)
+#
+# Time dependent version: f = du/dt + Biharmonic(u)
+u = exp(-t) * exp(-c * r**2)
+
+ut = diff(u,t)
+ur = diff(u,r)
+urr = diff(ur,r)
+u_theta = diff(u,theta)
+u_theta_theta = diff(u_theta,theta)
+
+# Gradient component
+gradu_r=ur
+gradu_t=(1/r)*u_theta
+# print('gradu_r = {}'.format(gradu_r))
+# print('gradu_t = {}'.format(gradu_t))
+
+lapu = simplify(urr + ur/r + u_theta_theta/r**2)
+print('lapu = {}'.format(lapu))
+
+# The biharmonic operator in polar coodinates has a bunch of nested
+# derivatives in r. Here we are ignoring the theta derivatives because
+# our Gaussian does not have any theta dependence.
+# https://en.wikipedia.org/wiki/Biharmonic_equation
+bihu = simplify((1/r) * diff(r * diff((1/r) * diff(r*diff(u,r),r),r),r))
+print('bihu = {}'.format(bihu))
+
+# 16*c**2*(c**2*r**4 - 4*c*r**2 + 2)*exp(-c*r**2 - t) - exp(-t)*exp(-c*r**2)
+f = ut + bihu
+print('Time-dependent forcing function = {}'.format(f))

test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_dirichlet.csv

@@ -0,0 +1,4 @@
+.1,8.411260e-04,3.207318e-02
+.05,5.885662e-05,4.515238e-03
+.025,3.721949e-06,5.738347e-04
+.0125,2.332737e-07,7.202340e-05

test/tests/indicators/laplacian_jump_indicator/biharmonic_centered_weak_bc.csv

@@ -0,0 +1,4 @@
+.1,8.411280e-04,3.207327e-02
+.05,5.886178e-05,4.515250e-03
+.025,3.722349e-06,5.738352e-04
+.0125,2.332996e-07,7.202341e-05

test/tests/indicators/laplacian_jump_indicator/biharmonic_transient.i

@@ -0,0 +1,122 @@
+[GlobalParams]
+  # Parameters used by Functions.
+  vars   = 'c'
+  vals   = '50'
+[]
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = -.5
+  xmax = .5
+  ymin = -.5
+  ymax = .5
+  nx = 10
+  ny = 10
+[]
+
+[Variables]
+  [./u]
+    order = THIRD
+    family = HERMITE
+  [../]
+[]
+
+[Kernels]
+  [./biharmonic]
+    type = Biharmonic
+    variable = u
+  [../]
+  [./body_force]
+    type = BodyForce
+    variable = u
+    function = forcing_func
+  [../]
+[]
+
+[BCs]
+  [./all_value]
+    type = FunctionPenaltyDirichletBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+  [./all_flux]
+    type = FunctionPenaltyFluxBC
+    variable = u
+    boundary = 'left right top bottom'
+    function = u_func
+    penalty = 1e10
+  [../]
+[]
+
+[Adaptivity]
+  [./Indicators]
+    [./error]
+      type = LaplacianJumpIndicator
+      variable = u
+      scale_by_flux_faces = true
+    [../]
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 4
+  dt = 0.1
+
+  # Note: the unusually tight tolerances here are due to the penalty
+  # BCs (currently the only way of accurately Dirichlet boundary
+  # conditions on Hermite elements in MOOSE).
+  nl_rel_tol = 1.e-15
+  l_tol = 1.e-15
+
+  # We have exact Jacobians
+  solve_type = 'NEWTON'
+
+  # Use 6x6 quadrature to ensure the forcing function is integrated
+  # accurately.
+  [./Quadrature]
+    type = GAUSS
+    order = ELEVENTH
+  [../]
+[]
+
+[Functions]
+  [./u_func]
+    type   = ParsedGradFunction
+    value  = 'exp(-c*(x^2+y^2))*exp(-t)'
+    grad_x = '-2*c*exp(-c*(x^2+y^2))*x*exp(-t)'
+    grad_y = '-2*c*exp(-c*(x^2+y^2))*y*exp(-t)'
+  [../]
+  [./forcing_func]
+    type   = ParsedFunction
+    value  = '16*c^2*(c^2*(x^2+y^2)^2 - 4*c*(x^2+y^2) + 2)*exp(-c*(x^2+y^2))*exp(-t)'
+  [../]
+[]
+
+[ICs]
+  [./u_ic]
+    type = FunctionIC
+    function = u_func
+    variable = u
+  [../]
+[]
+
+[Postprocessors]
+  [./l2_error]
+    type = ElementL2Error
+    variable = u
+    function = u_func
+  [../]
+  [./h1_error]
+    type = ElementH1Error
+    variable = u
+    function = u_func
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]