Add upper bound test

Refs #2999

framework/include/nodalkernels/UpperBoundNodalKernel.h

@@ -0,0 +1,45 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "NodalKernel.h"
+
+// Forward Declarations
+class UpperBoundNodalKernel;
+
+template <>
+InputParameters validParams<UpperBoundNodalKernel>();
+
+/**
+ * Class used to enforce a upper bound on a coupled variable
+ */
+class UpperBoundNodalKernel : public NodalKernel
+{
+public:
+  UpperBoundNodalKernel(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+private:
+  /// The number of the coupled variable
+  const unsigned int _v_var;
+
+  /// The value of the coupled variable
+  const VariableValue & _v;
+
+  /// Boundaries on which we should not execute this object
+  std::set<BoundaryID> _bnd_ids;
+
+  /// The upper bound on the coupled variable
+  const Real _upper_bound;
+};

framework/src/nodalkernels/UpperBoundNodalKernel.C

@@ -0,0 +1,80 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "UpperBoundNodalKernel.h"
+
+registerMooseObject("MooseApp", UpperBoundNodalKernel);
+
+template <>
+InputParameters
+validParams<UpperBoundNodalKernel>()
+{
+  InputParameters params = validParams<NodalKernel>();
+  params.addClassDescription("Used to prevent a coupled variable from going above a upper bound");
+  params.addRequiredCoupledVar(
+      "v", "The coupled variable we require to be greater than the upper bound");
+  params.addParam<Real>("upper_bound", 0, "The upper bound on the coupled variable");
+  params.addParam<std::vector<BoundaryName>>(
+      "exclude_boundaries",
+      "Boundaries on which not to execute the NodalKernel. This can be useful for avoiding "
+      "singuarility in the matrix in case a constraint is active in the same place that a "
+      "DirichletBC is set");
+  return params;
+}
+
+UpperBoundNodalKernel::UpperBoundNodalKernel(const InputParameters & parameters)
+  : NodalKernel(parameters),
+    _v_var(coupled("v")),
+    _v(coupledValue("v")),
+    _upper_bound(getParam<Real>("upper_bound"))
+{
+  if (_var.number() == _v_var)
+    mooseError("Coupled variable 'v' needs to be different from 'variable' with "
+               "UpperBoundNodalKernel");
+
+  const auto & bnd_names = getParam<std::vector<BoundaryName>>("exclude_boundaries");
+  for (const auto & bnd_name : bnd_names)
+    _bnd_ids.insert(_mesh.getBoundaryID(bnd_name));
+}
+
+Real
+UpperBoundNodalKernel::computeQpResidual()
+{
+  for (auto bnd_id : _bnd_ids)
+    if (_mesh.isBoundaryNode(_current_node->id(), bnd_id))
+      return _u[_qp];
+
+  return std::min(_u[_qp], _upper_bound - _v[_qp]);
+}
+
+Real
+UpperBoundNodalKernel::computeQpJacobian()
+{
+  for (auto bnd_id : _bnd_ids)
+    if (_mesh.isBoundaryNode(_current_node->id(), bnd_id))
+      return 1;
+
+  if (_u[_qp] <= _upper_bound - _v[_qp])
+    return 1;
+  return 0;
+}
+
+Real
+UpperBoundNodalKernel::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  for (auto bnd_id : _bnd_ids)
+    if (_mesh.isBoundaryNode(_current_node->id(), bnd_id))
+      return 0;
+
+  if (jvar == _v_var)
+    if (_upper_bound - _v[_qp] < _u[_qp])
+      return -1;
+
+  return 0.0;
+}

test/tests/nodalkernels/constraint_enforcement/lower-bound.i

@@ -0,0 +1,122 @@
+l=10
+nx=100
+num_steps=10
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 1
+  xmax = ${l}
+  nx = ${nx}
+[]
+
+[Variables]
+  [u]
+  []
+  [lm]
+  []
+[]
+
+[ICs]
+  [u]
+    type = FunctionIC
+    variable = u
+    function = '${l} - x'
+  []
+[]
+
+[Kernels]
+  [time]
+    type = TimeDerivative
+    variable = u
+  []
+  [diff]
+    type = Diffusion
+    variable = u
+  []
+  [ffn]
+    type = BodyForce
+    variable = u
+    function = '-1'
+  []
+[]
+
+[NodalKernels]
+  [positive_constraint]
+    type = LowerBoundNodalKernel
+    variable = lm
+    v = u
+    exclude_boundaries = 'left right'
+  []
+  [forces]
+    type = CoupledForceNodalKernel
+    variable = u
+    v = lm
+  []
+[]
+
+
+[BCs]
+  [left]
+    type = DirichletBC
+    boundary = left
+    value = ${l}
+    variable = u
+  []
+  [right]
+    type = DirichletBC
+    boundary = right
+    value = 0
+    variable = u
+  []
+[]
+
+[NodalKernels]
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = ${num_steps}
+  solve_type = NEWTON
+  dtmin = 1
+  petsc_options_iname = '-snes_max_linear_solve_fail -ksp_max_it -pc_type -sub_pc_factor_levels -snes_linesearch_type'
+  petsc_options_value = '0                           30          asm      16                    basic'
+[]
+
+[Outputs]
+  exodus = true
+  [csv]
+    type = CSV
+    execute_on = 'nonlinear timestep_end'
+  []
+  [dof]
+    type = DOFMap
+    execute_on = 'initial'
+  []
+[]
+
+[Debug]
+  show_var_residual_norms = true
+[]
+
+[Postprocessors]
+  [active_lm]
+    type = LMActiveSetSize
+    variable = lm
+    execute_on = 'nonlinear timestep_end'
+    value = 1e-8
+  []
+  [violations]
+    type = LMActiveSetSize
+    variable = u
+    execute_on = 'nonlinear timestep_end'
+    value = -1e-8
+    comparator = 'less'
+  []
+[]

test/tests/nodalkernels/constraint_enforcement/tests

@@ -1,17 +1,39 @@
 [Tests]
-  [exo]
-    input = interpolated-ncp-lm-nodal-enforcement-nodal-forces.i
-    type = Exodiff
-    exodiff = interpolated-ncp-lm-nodal-enforcement-nodal-forces_out.e
-    requirement = 'The system shall be able to solve a positively constrained PDE using nodal NCP, nodal application of resultant forces, and have no oscillations in the solution'
+  issues = '#2999'
+  [lower_bound]
+    design = LowerBoundNodalKernel.md
+    requirement = 'The system shall be able to enforce a lower bound on a variable using nodal NCP, nodal application of resultant forces,'
+    [exo]
+      input = lower-bound.i
+      type = Exodiff
+      exodiff = lower-bound_out.e
+      detail = 'have no oscillations in the solution, and'
+    []
+    [non_singular]
+      input = lower-bound.i
+      type = RunApp
+      absent_out = '[1-9]+[0-9]* of 22 singular values'
+      expect_out = '0 of 22 singular values'
+      cli_args = "nx=10 num_steps=5 Outputs/exodus=false Outputs/active='' -pc_type svd -pc_svd_monitor"
+      detail = 'have a non-singular matrix'
+    []
   []
-  [non_singular]
-    input = interpolated-ncp-lm-nodal-enforcement-nodal-forces.i
-    type = RunApp
-    prereq = exo
-    requirement = 'The system shall be able to solve a positively constrained PDE using nodal NCP, nodal application of resultant forces, and have a non-singular matrix'
-    absent_out = '[1-9]+[0-9]* of 22 singular values'
-    expect_out = '0 of 22 singular values'
-    cli_args = 'nx=10 num_steps=5 -pc_type svd -pc_svd_monitor'
+  [upper_bound]
+    design = UpperBoundNodalKernel.md
+    requirement = 'The system shall be able to enforce an upper bound on a variable using nodal NCP, nodal application of resultant forces,'
+    [exo]
+      input = upper-bound.i
+      type = Exodiff
+      exodiff = upper-bound_out.e
+      detail = 'have no oscillations in the solution, and'
+    []
+    [non_singular]
+      input = upper-bound.i
+      type = RunApp
+      absent_out = '[1-9]+[0-9]* of 22 singular values'
+      expect_out = '0 of 22 singular values'
+      cli_args = "nx=10 num_steps=5 Outputs/exodus=false Outputs/active='' -pc_type svd -pc_svd_monitor"
+      detail = 'have a non-singular matrix'
+    []
   []
 []