diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..d988ede
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,64 @@
+# Adapted from the RACCOON project's test workflow:
+# https://github.com/hugary1995/raccoon/blob/devel/.github/workflows/tests.yml
+
+name: Tests
+
+on:
+  # Triggers the workflow on pushes to the master branch
+  push:
+    branches: [ master ]
+  pull_request:
+    branches: [ master ]
+
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
+
+env:
+  # Squirrel does not carry MOOSE as a submodule, so CI clones it next to the
+  # checkout.  Pin this to the MOOSE commit the downstream apps use whenever
+  # master drifts away from squirrel.
+  MOOSE_REF: master
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+
+    defaults:
+      run:
+        shell: bash -el {0}
+
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          path: squirrel
+
+      - name: Clone MOOSE
+        run: |
+          git clone --depth 1 --branch "${MOOSE_REF}" https://github.com/idaholab/moose.git moose
+
+      - name: Setup environment
+        uses: mamba-org/setup-micromamba@v2
+        with:
+          environment-name: moose
+          condarc: |
+            channels:
+              - https://conda.software.inl.gov/public
+              - conda-forge
+            channel_priority: flexible
+          init-shell: bash
+          create-args: moose-dev mpich
+        continue-on-error: false
+
+      - name: Compile squirrel
+        working-directory: squirrel
+        env:
+          MOOSE_DIR: ${{ github.workspace }}/moose
+        run: |
+          make -j 4
+
+      - name: Regression tests
+        working-directory: squirrel
+        env:
+          MOOSE_DIR: ${{ github.workspace }}/moose
+        run: |
+          ./run_tests -j 4
diff --git a/.gitignore b/.gitignore
index 9082198..1c4bc34 100644
--- a/.gitignore
+++ b/.gitignore
@@ -152,3 +152,6 @@ squirrel.yaml
 
 # VSCode dotfolder
 .vscode
+
+# JIT files
+**/.jitcache/*
diff --git a/include/bcs/PostprocessorVelocityFunctionTemperatureInflowBC.h b/include/bcs/PostprocessorVelocityFunctionTemperatureInflowBC.h
new file mode 100644
index 0000000..61a56d9
--- /dev/null
+++ b/include/bcs/PostprocessorVelocityFunctionTemperatureInflowBC.h
@@ -0,0 +1,24 @@
+#pragma once
+
+#include "PostprocessorVelocityFunctionInflowBC.h"
+#include "JvarMapInterface.h"
+#include "DerivativeMaterialInterface.h"
+
+class PostprocessorVelocityFunctionTemperatureInflowBC
+    : public DerivativeMaterialInterface<JvarMapIntegratedBCInterface<PostprocessorVelocityFunctionInflowBC>>
+{
+public:
+  PostprocessorVelocityFunctionTemperatureInflowBC(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  virtual void initialSetup() override;
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+
+  const MaterialProperty<Real> & _rho;
+  const MaterialProperty<Real> & _d_rho_d_u;
+  const MaterialProperty<Real> & _cp;
+  const MaterialProperty<Real> & _d_cp_d_u;
+};
diff --git a/src/bcs/DiffusiveFluxBC.C b/src/bcs/DiffusiveFluxBC.C
index c6426bc..fa93a9d 100644
--- a/src/bcs/DiffusiveFluxBC.C
+++ b/src/bcs/DiffusiveFluxBC.C
@@ -27,12 +27,12 @@ DiffusiveFluxBC::initialSetup()
 Real
 DiffusiveFluxBC::computeQpResidual()
 {
-  return _normals[_qp] * -_D[_qp] * _grad_u[_qp];
+  return -_test[_i][_qp] * _normals[_qp] * _D[_qp] * _grad_u[_qp];
 }
 
 Real
 DiffusiveFluxBC::computeQpJacobian()
 {
-  return -_normals[_qp] *
+  return -_test[_i][_qp] * _normals[_qp] *
          (_d_D_d_u[_qp] * _phi[_j][_qp] * _grad_u[_qp] + _D[_qp] * _grad_phi[_j][_qp]);
 }
diff --git a/src/bcs/PostprocessorTemperatureInflowBC.C b/src/bcs/PostprocessorTemperatureInflowBC.C
index e478639..377e1e9 100644
--- a/src/bcs/PostprocessorTemperatureInflowBC.C
+++ b/src/bcs/PostprocessorTemperatureInflowBC.C
@@ -6,6 +6,8 @@ InputParameters
 PostprocessorTemperatureInflowBC::validParams()
 {
   InputParameters params = PostprocessorInflowBC::validParams();
+  params.addParam<MaterialPropertyName>("rho", "rho", "Density material property name");
+  params.addParam<MaterialPropertyName>("cp", "cp", "Specific heat material property name");
   return params;
 }
 
diff --git a/src/bcs/PostprocessorVelocityFunctionTemperatureInflowBC.C b/src/bcs/PostprocessorVelocityFunctionTemperatureInflowBC.C
new file mode 100644
index 0000000..e87f5e2
--- /dev/null
+++ b/src/bcs/PostprocessorVelocityFunctionTemperatureInflowBC.C
@@ -0,0 +1,43 @@
+#include "PostprocessorVelocityFunctionTemperatureInflowBC.h"
+
+registerMooseObject("SquirrelApp", PostprocessorVelocityFunctionTemperatureInflowBC);
+
+InputParameters
+PostprocessorVelocityFunctionTemperatureInflowBC::validParams()
+{
+  InputParameters params = PostprocessorVelocityFunctionInflowBC::validParams();
+  params.addParam<MaterialPropertyName>("rho", "rho", "Density material property name");
+  params.addParam<MaterialPropertyName>("cp", "cp", "Specific heat material property name");
+  return params;
+}
+
+PostprocessorVelocityFunctionTemperatureInflowBC::PostprocessorVelocityFunctionTemperatureInflowBC(
+    const InputParameters & parameters)
+  : DerivativeMaterialInterface<JvarMapIntegratedBCInterface<PostprocessorVelocityFunctionInflowBC>>(parameters),
+    _rho(getMaterialProperty<Real>("rho")),
+    _d_rho_d_u(getMaterialPropertyDerivative<Real>("rho", _var.name())),
+    _cp(getMaterialProperty<Real>("cp")),
+    _d_cp_d_u(getMaterialPropertyDerivative<Real>("cp", _var.name()))
+{
+}
+
+void
+PostprocessorVelocityFunctionTemperatureInflowBC::initialSetup()
+{
+  validateNonlinearCoupling<Real>("rho");
+  validateNonlinearCoupling<Real>("cp");
+}
+
+Real
+PostprocessorVelocityFunctionTemperatureInflowBC::computeQpResidual()
+{
+  return _rho[_qp] * _cp[_qp] * PostprocessorVelocityFunctionInflowBC::computeQpResidual();
+}
+
+Real
+PostprocessorVelocityFunctionTemperatureInflowBC::computeQpJacobian()
+{
+  return _rho[_qp] * _cp[_qp] * PostprocessorVelocityFunctionInflowBC::computeQpJacobian() +
+         _d_rho_d_u[_qp] * _phi[_j][_qp] * _cp[_qp] * PostprocessorVelocityFunctionInflowBC::computeQpResidual() +
+         _rho[_qp] * _d_cp_d_u[_qp] * _phi[_j][_qp] * PostprocessorVelocityFunctionInflowBC::computeQpResidual();
+}
diff --git a/src/bcs/VelocityFunctionTemperatureOutflowBC.C b/src/bcs/VelocityFunctionTemperatureOutflowBC.C
index 8d5fdca..409f587 100644
--- a/src/bcs/VelocityFunctionTemperatureOutflowBC.C
+++ b/src/bcs/VelocityFunctionTemperatureOutflowBC.C
@@ -10,6 +10,8 @@ VelocityFunctionTemperatureOutflowBC::validParams()
   params.addRequiredParam<FunctionName>("vel_x_func", "The x velocity function");
   params.addRequiredParam<FunctionName>("vel_y_func", "The y velocity function");
   params.addRequiredParam<FunctionName>("vel_z_func", "The z velocity function");
+  params.addParam<MaterialPropertyName>("rho", "rho", "Density material property name");
+  params.addParam<MaterialPropertyName>("cp", "cp", "Specific heat material property name");
   return params;
 }
 
diff --git a/src/dgkernels/DGFunctionTemperatureAdvection.C b/src/dgkernels/DGFunctionTemperatureAdvection.C
index 63adbc1..b75eb09 100644
--- a/src/dgkernels/DGFunctionTemperatureAdvection.C
+++ b/src/dgkernels/DGFunctionTemperatureAdvection.C
@@ -7,6 +7,8 @@ DGFunctionTemperatureAdvection::validParams()
 {
   InputParameters params = DGFunctionConvection::validParams();
   params.addClassDescription("DG upwinding for func temp convection");
+  params.addParam<MaterialPropertyName>("rho", "rho", "Density material property name");
+  params.addParam<MaterialPropertyName>("cp", "cp", "Specific heat material property name");
   return params;
 }
 
diff --git a/src/interface_kernels/InterTemperatureAdvection.C b/src/interface_kernels/InterTemperatureAdvection.C
index 24a0945..75fd53f 100644
--- a/src/interface_kernels/InterTemperatureAdvection.C
+++ b/src/interface_kernels/InterTemperatureAdvection.C
@@ -42,12 +42,12 @@ InterTemperatureAdvection::computeQpResidual(Moose::DGResidualType type)
       if (vdotn >= 0)
         r += vdotn * _u[_qp] * _test[_i][_qp];
       else
-        r += vdotn * _neighbor_value[_qp] * _test[_i][_qp] - _heat_source;
+        r += (vdotn * _neighbor_value[_qp] - _heat_source) * _test[_i][_qp];
       break;
 
     case Moose::Neighbor:
       if (vdotn >= 0)
-        r -= vdotn * _u[_qp] * _test_neighbor[_i][_qp] - _heat_source;
+        r -= (vdotn * _u[_qp] - _heat_source) * _test_neighbor[_i][_qp];
       else
         r -= vdotn * _neighbor_value[_qp] * _test_neighbor[_i][_qp];
       break;
diff --git a/src/kernels/MatINSTemperatureTimeDerivative.C b/src/kernels/MatINSTemperatureTimeDerivative.C
index c7af700..d172ee4 100644
--- a/src/kernels/MatINSTemperatureTimeDerivative.C
+++ b/src/kernels/MatINSTemperatureTimeDerivative.C
@@ -34,7 +34,7 @@ MatINSTemperatureTimeDerivative::computeQpResidual()
 Real
 MatINSTemperatureTimeDerivative::computeQpJacobian()
 {
-  return _rho[_qp] * _cp[_qp] * TimeDerivative::computeQpJacobian(); // +
-        //  _d_rho_d_u[_qp] * _phi[_j][_qp] * _cp[_qp] * TimeDerivative::computeQpResidual() +
-        //  _rho[_qp] * _d_cp_d_u[_qp] * _phi[_j][_qp] * TimeDerivative::computeQpResidual();
+  return _rho[_qp] * _cp[_qp] * TimeDerivative::computeQpJacobian() +
+         _d_rho_d_u[_qp] * _phi[_j][_qp] * _cp[_qp] * TimeDerivative::computeQpResidual() +
+         _rho[_qp] * _d_cp_d_u[_qp] * _phi[_j][_qp] * TimeDerivative::computeQpResidual();
 }
diff --git a/src/kernels/VelocityFunctionTemperatureAdvection.C b/src/kernels/VelocityFunctionTemperatureAdvection.C
index b1a1d8d..39fe9d3 100644
--- a/src/kernels/VelocityFunctionTemperatureAdvection.C
+++ b/src/kernels/VelocityFunctionTemperatureAdvection.C
@@ -10,6 +10,8 @@ VelocityFunctionTemperatureAdvection::validParams()
   params.addRequiredParam<FunctionName>("vel_x_func", "The x velocity function");
   params.addRequiredParam<FunctionName>("vel_y_func", "The y velocity function");
   params.addRequiredParam<FunctionName>("vel_z_func", "The z velocity function");
+  params.addParam<MaterialPropertyName>("rho", "rho", "Density material property name");
+  params.addParam<MaterialPropertyName>("cp", "cp", "Specific heat material property name");
   return params;
 }
 
diff --git a/tests/auxkernels/density_from_log/density_from_log.i b/tests/auxkernels/density_from_log/density_from_log.i
new file mode 100644
index 0000000..5282248
--- /dev/null
+++ b/tests/auxkernels/density_from_log/density_from_log.i
@@ -0,0 +1,58 @@
+# DensityFromLog exponentiates a log-density field. Here the log-density falls
+# linearly with x, so the recovered density is rho = exp(-x/H).
+
+H = 0.5
+
+[Mesh]
+  [column]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 4
+    xmax = 1
+  []
+[]
+
+[Problem]
+  solve = false
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [dummy]
+  []
+[]
+
+[AuxVariables]
+  [log_rho]
+  []
+  [rho]
+  []
+[]
+
+[Functions]
+  [log_rho_func]
+    type = ParsedFunction
+    expression = '-x / ${H}'
+  []
+[]
+
+[AuxKernels]
+  [log_rho]
+    type = FunctionAux
+    variable = log_rho
+    function = log_rho_func
+  []
+  [rho]
+    type = DensityFromLog
+    variable = rho
+    density_log = log_rho
+  []
+[]
+
+[Executioner]
+  type = Steady
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/auxkernels/density_from_log/gold/density_from_log_out.e b/tests/auxkernels/density_from_log/gold/density_from_log_out.e
new file mode 100644
index 0000000..8634549
Binary files /dev/null and b/tests/auxkernels/density_from_log/gold/density_from_log_out.e differ
diff --git a/tests/auxkernels/density_from_log/gold/larger_scale_height_out.e b/tests/auxkernels/density_from_log/gold/larger_scale_height_out.e
new file mode 100644
index 0000000..48377a6
Binary files /dev/null and b/tests/auxkernels/density_from_log/gold/larger_scale_height_out.e differ
diff --git a/tests/auxkernels/density_from_log/tests b/tests/auxkernels/density_from_log/tests
new file mode 100644
index 0000000..d9786cf
--- /dev/null
+++ b/tests/auxkernels/density_from_log/tests
@@ -0,0 +1,16 @@
+[Tests]
+  [density_from_log]
+    type = Exodiff
+    input = 'density_from_log.i'
+    exodiff = 'density_from_log_out.e'
+    requirement = 'The system shall recover a density field by exponentiating a log-density variable.'
+  []
+  [larger_scale_height]
+    type = Exodiff
+    input = 'density_from_log.i'
+    cli_args = 'H=2
+                Outputs/file_base=larger_scale_height_out'
+    exodiff = 'larger_scale_height_out.e'
+    requirement = 'The system shall flatten the recovered density profile as the log-density scale height increases.'
+  []
+[]
diff --git a/tests/auxkernels/function_derivative_aux/function_derivative_aux.i b/tests/auxkernels/function_derivative_aux/function_derivative_aux.i
new file mode 100644
index 0000000..3bed7e6
--- /dev/null
+++ b/tests/auxkernels/function_derivative_aux/function_derivative_aux.i
@@ -0,0 +1,61 @@
+# FunctionDerivativeAux samples a component of a function gradient. For the
+# field u(x) = x^2 the x derivative is 2x, evaluated both at nodes and at
+# element quadrature points.
+
+[Mesh]
+  [rod]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 4
+    xmax = 1
+  []
+[]
+
+[Problem]
+  solve = false
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [dummy]
+  []
+[]
+
+[Functions]
+  [quadratic]
+    type = ParsedFunction
+    expression = 'x * x'
+  []
+[]
+
+[AuxVariables]
+  [dudx_nodal]
+  []
+  [dudx_elemental]
+    order = CONSTANT
+    family = MONOMIAL
+  []
+[]
+
+[AuxKernels]
+  [nodal]
+    type = FunctionDerivativeAux
+    variable = dudx_nodal
+    function = quadratic
+    component = 1
+  []
+  [elemental]
+    type = FunctionDerivativeAux
+    variable = dudx_elemental
+    function = quadratic
+    component = 1
+  []
+[]
+
+[Executioner]
+  type = Steady
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/auxkernels/function_derivative_aux/gold/function_derivative_aux_out.e b/tests/auxkernels/function_derivative_aux/gold/function_derivative_aux_out.e
new file mode 100644
index 0000000..9557c1e
Binary files /dev/null and b/tests/auxkernels/function_derivative_aux/gold/function_derivative_aux_out.e differ
diff --git a/tests/auxkernels/function_derivative_aux/gold/y_derivative_out.e b/tests/auxkernels/function_derivative_aux/gold/y_derivative_out.e
new file mode 100644
index 0000000..90dbc3f
Binary files /dev/null and b/tests/auxkernels/function_derivative_aux/gold/y_derivative_out.e differ
diff --git a/tests/auxkernels/function_derivative_aux/tests b/tests/auxkernels/function_derivative_aux/tests
new file mode 100644
index 0000000..2514883
--- /dev/null
+++ b/tests/auxkernels/function_derivative_aux/tests
@@ -0,0 +1,17 @@
+[Tests]
+  [x_derivative]
+    type = Exodiff
+    input = 'function_derivative_aux.i'
+    exodiff = 'function_derivative_aux_out.e'
+    requirement = 'The system shall evaluate a component of a function gradient at nodes and at element quadrature points.'
+  []
+  [y_derivative]
+    type = Exodiff
+    input = 'function_derivative_aux.i'
+    cli_args = 'AuxKernels/nodal/component=2
+                AuxKernels/elemental/component=2
+                Outputs/file_base=y_derivative_out'
+    exodiff = 'y_derivative_out.e'
+    requirement = 'The system shall return a zero derivative for the direction in which the sampled function is constant.'
+  []
+[]
diff --git a/tests/bcs/channel_gradient/channel_gradient.i b/tests/bcs/channel_gradient/channel_gradient.i
new file mode 100644
index 0000000..7c9a327
--- /dev/null
+++ b/tests/bcs/channel_gradient/channel_gradient.i
@@ -0,0 +1,145 @@
+# A heated channel T alongside a coolant channel at a prescribed temperature.
+# ChannelGradient forms the axial profile of the temperature difference, and
+# ChannelGradientBC applies wall cooling h*(T - T_coolant) from that profile.
+
+Nu = 4
+k = 1
+half_channel_width = 0.5
+h = ${fparse Nu * k / half_channel_width}
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 2
+    ny = 10
+    xmax = 1
+    ymax = 1
+  []
+[]
+
+[Variables]
+  [T]
+    initial_condition = 2
+  []
+[]
+
+[AuxVariables]
+  [T_coolant]
+  []
+[]
+
+[AuxKernels]
+  [T_coolant]
+    type = FunctionAux
+    variable = T_coolant
+    function = coolant_temperature
+  []
+[]
+
+[Functions]
+  [coolant_temperature]
+    type = ParsedFunction
+    expression = '1 + y'
+  []
+[]
+
+[Kernels]
+  [advection]
+    type = ConservativeAdvection
+    variable = T
+    velocity_variable = '0 1 0'
+  []
+  [conduction]
+    type = Diffusion
+    variable = T
+  []
+  [heat]
+    type = BodyForce
+    variable = T
+    value = 1
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = DirichletBC
+    variable = T
+    boundary = bottom
+    value = 2
+  []
+  [outlet]
+    type = OutflowBC
+    variable = T
+    boundary = top
+    velocity = '0 1 0'
+  []
+  [coolant_wall]
+    type = ChannelGradientBC
+    variable = T
+    boundary = right
+    channel_gradient_pps = channel_gradient
+    axis = y
+    h_name = h
+  []
+[]
+
+[Materials]
+  [heat_transfer_coefficient]
+    type = GenericConstantMaterial
+    prop_names = 'h'
+    prop_values = '${h}'
+  []
+[]
+
+[VectorPostprocessors]
+  [fuel_line]
+    type = LineValueSampler
+    variable = T
+    start_point = '0 0 0'
+    end_point = '0 1 0'
+    num_points = 5
+    sort_by = y
+    execute_on = 'initial timestep_begin linear nonlinear timestep_end'
+    outputs = none
+  []
+  [coolant_line]
+    type = LineValueSampler
+    variable = T_coolant
+    start_point = '1 0 0'
+    end_point = '1 1 0'
+    num_points = 5
+    sort_by = y
+    execute_on = 'initial timestep_begin linear nonlinear timestep_end'
+    outputs = none
+  []
+  [channel_gradient]
+    type = ChannelGradient
+    lv1 = fuel_line
+    lv2 = coolant_line
+    var1 = T
+    var2 = T_coolant
+    axis = y
+    execute_on = 'initial timestep_begin linear nonlinear timestep_end'
+    outputs = csv
+  []
+[]
+
+[Postprocessors]
+  [T_outlet]
+    type = SideAverageValue
+    variable = T
+    boundary = top
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  csv = true
+[]
diff --git a/tests/bcs/channel_gradient/gold/adiabatic_wall_out.csv b/tests/bcs/channel_gradient/gold/adiabatic_wall_out.csv
new file mode 100644
index 0000000..f581644
--- /dev/null
+++ b/tests/bcs/channel_gradient/gold/adiabatic_wall_out.csv
@@ -0,0 +1,3 @@
+time,T_outlet
+0,0
+1,2.3675725423829
diff --git a/tests/bcs/channel_gradient/gold/adiabatic_wall_out_channel_gradient_0001.csv b/tests/bcs/channel_gradient/gold/adiabatic_wall_out_channel_gradient_0001.csv
new file mode 100644
index 0000000..a5ce698
--- /dev/null
+++ b/tests/bcs/channel_gradient/gold/adiabatic_wall_out_channel_gradient_0001.csv
@@ -0,0 +1,6 @@
+gradient,y
+1,0
+0.89491037241394,0.25
+0.76129493073712,0.5
+0.58795910970821,0.75
+0.36757254238287,1
diff --git a/tests/bcs/channel_gradient/gold/channel_gradient_out.csv b/tests/bcs/channel_gradient/gold/channel_gradient_out.csv
new file mode 100644
index 0000000..b5c97d8
--- /dev/null
+++ b/tests/bcs/channel_gradient/gold/channel_gradient_out.csv
@@ -0,0 +1,3 @@
+time,T_outlet
+0,0
+1,1.5874816473487
diff --git a/tests/bcs/channel_gradient/gold/channel_gradient_out_channel_gradient_0001.csv b/tests/bcs/channel_gradient/gold/channel_gradient_out_channel_gradient_0001.csv
new file mode 100644
index 0000000..a72b216
--- /dev/null
+++ b/tests/bcs/channel_gradient/gold/channel_gradient_out_channel_gradient_0001.csv
@@ -0,0 +1,6 @@
+gradient,y
+1,0
+0.72851598013405,0.25
+0.4382658663882,0.5
+0.12745947296864,0.75
+-0.1485619113861,1
diff --git a/tests/bcs/channel_gradient/tests b/tests/bcs/channel_gradient/tests
new file mode 100644
index 0000000..97171c8
--- /dev/null
+++ b/tests/bcs/channel_gradient/tests
@@ -0,0 +1,18 @@
+[Tests]
+  [channel_gradient]
+    type = CSVDiff
+    input = 'channel_gradient.i'
+    csvdiff = 'channel_gradient_out.csv
+               channel_gradient_out_channel_gradient_0001.csv'
+    requirement = 'The system shall compute the axial profile of the difference between two line samplers and apply it as a wall heat transfer boundary condition.'
+  []
+  [adiabatic_wall]
+    type = CSVDiff
+    input = 'channel_gradient.i'
+    cli_args = 'Nu=0
+                Outputs/file_base=adiabatic_wall_out'
+    csvdiff = 'adiabatic_wall_out.csv
+               adiabatic_wall_out_channel_gradient_0001.csv'
+    requirement = 'The system shall recover an insulated wall when the heat transfer coefficient of the channel gradient boundary condition is zero.'
+  []
+[]
diff --git a/tests/bcs/diffusive_flux/diffusive_flux.i b/tests/bcs/diffusive_flux/diffusive_flux.i
new file mode 100644
index 0000000..29c3d44
--- /dev/null
+++ b/tests/bcs/diffusive_flux/diffusive_flux.i
@@ -0,0 +1,96 @@
+# Heated flow down a channel: -k T'' + v T' = q, T(0) = 1, with the enthalpy
+# leaving through an OutflowBC. The conductive part of the outlet flux is
+# supplied by one of two boundary conditions selected by the tests.
+
+k = 0.1
+v = 1
+q = 1
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [T]
+    initial_condition = 1
+  []
+[]
+
+[Kernels]
+  [conduction]
+    type = MatDiffusion
+    variable = T
+    diffusivity = k
+  []
+  [advection]
+    type = ConservativeAdvection
+    variable = T
+    velocity_variable = '${v} 0 0'
+  []
+  [heating]
+    type = BodyForce
+    variable = T
+    value = ${q}
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = DirichletBC
+    variable = T
+    boundary = left
+    value = 1
+  []
+  [outlet_advection]
+    type = OutflowBC
+    variable = T
+    boundary = right
+    velocity = '${v} 0 0'
+  []
+  [outlet_conduction]
+    type = MatINSTemperatureNoBCBC
+    variable = T
+    boundary = right
+    k = k
+  []
+  [outlet_diffusive_flux]
+    type = DiffusiveFluxBC
+    variable = T
+    boundary = right
+    D_name = k
+    enable = false
+  []
+[]
+
+[Materials]
+  [conductivity]
+    type = DerivativeParsedMaterial
+    property_name = k
+    coupled_variables = 'T'
+    expression = '${k}'
+    derivative_order = 1
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/bcs/diffusive_flux/gold/diffusive_flux_bc_out.e b/tests/bcs/diffusive_flux/gold/diffusive_flux_bc_out.e
new file mode 100644
index 0000000..de88fea
Binary files /dev/null and b/tests/bcs/diffusive_flux/gold/diffusive_flux_bc_out.e differ
diff --git a/tests/bcs/diffusive_flux/gold/diffusive_flux_out.e b/tests/bcs/diffusive_flux/gold/diffusive_flux_out.e
new file mode 100644
index 0000000..582e5bc
Binary files /dev/null and b/tests/bcs/diffusive_flux/gold/diffusive_flux_out.e differ
diff --git a/tests/bcs/diffusive_flux/tests b/tests/bcs/diffusive_flux/tests
new file mode 100644
index 0000000..e4b9e94
--- /dev/null
+++ b/tests/bcs/diffusive_flux/tests
@@ -0,0 +1,27 @@
+[Tests]
+  [mat_ins_temperature_no_bc]
+    type = Exodiff
+    input = 'diffusive_flux.i'
+    exodiff = 'diffusive_flux_out.e'
+    requirement = 'The system shall reintroduce the conductive surface term at an outflow boundary so that heat conducts out consistently with the interior temperature gradient.'
+  []
+  [diffusive_flux]
+    type = Exodiff
+    input = 'diffusive_flux.i'
+    cli_args = 'BCs/outlet_conduction/enable=false
+                BCs/outlet_diffusive_flux/enable=true
+                Outputs/file_base=diffusive_flux_bc_out'
+    exodiff = 'diffusive_flux_bc_out.e'
+    requirement = 'The system shall apply a diffusive flux at a boundary equal to the conductivity times the outward normal temperature gradient.'
+  []
+  [diffusive_flux_jacobian]
+    type = PetscJacobianTester
+    input = 'diffusive_flux.i'
+    run_sim = true
+    cli_args = 'BCs/outlet_conduction/enable=false
+                BCs/outlet_diffusive_flux/enable=true
+                k=0.1*(1+0.5*T)
+                Outputs/exodus=false'
+    requirement = 'The system shall provide an exact Jacobian for the diffusive flux boundary condition with a temperature-dependent conductivity.'
+  []
+[]
diff --git a/tests/bcs/example_shape_side/example_shape_side.i b/tests/bcs/example_shape_side/example_shape_side.i
new file mode 100644
index 0000000..1b5b0ec
--- /dev/null
+++ b/tests/bcs/example_shape_side/example_shape_side.i
@@ -0,0 +1,97 @@
+# Nonlocal electrode boundary condition. The species density v is pinned to 1
+# and 2 across the gap, so at the right electrode the user objects supply the
+# surface integrals int(grad(v).n) = 1 and int(v) = 2, floating the electrode
+# potential to (Vb - 1) / 3 while the potential stays linear across the gap.
+
+Vb = 3
+
+[Mesh]
+  [gap]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 4
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [potential]
+  []
+  [v]
+    initial_condition = 1
+  []
+[]
+
+[Kernels]
+  [potential_diffusion]
+    type = Diffusion
+    variable = potential
+  []
+  [v_diffusion]
+    type = Diffusion
+    variable = v
+  []
+[]
+
+[BCs]
+  [grounded_electrode]
+    type = DirichletBC
+    variable = potential
+    boundary = left
+    value = 0
+  []
+  [floating_electrode]
+    type = ExampleShapeSideIntegratedBC
+    variable = potential
+    boundary = right
+    v = v
+    Vb = ${Vb}
+    num_user_object = num_integral
+    denom_user_object = denom_integral
+  []
+  [v_left]
+    type = DirichletBC
+    variable = v
+    boundary = left
+    value = 1
+  []
+  [v_right]
+    type = DirichletBC
+    variable = v
+    boundary = right
+    value = 2
+  []
+[]
+
+[UserObjects]
+  [num_integral]
+    type = NumShapeSideUserObject
+    u = v
+    boundary = right
+    execute_on = 'linear nonlinear'
+  []
+  [denom_integral]
+    type = DenomShapeSideUserObject
+    u = v
+    boundary = right
+    execute_on = 'linear nonlinear'
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/bcs/example_shape_side/gold/example_shape_side_out.e b/tests/bcs/example_shape_side/gold/example_shape_side_out.e
new file mode 100644
index 0000000..49147e4
Binary files /dev/null and b/tests/bcs/example_shape_side/gold/example_shape_side_out.e differ
diff --git a/tests/bcs/example_shape_side/gold/larger_bias_out.e b/tests/bcs/example_shape_side/gold/larger_bias_out.e
new file mode 100644
index 0000000..3221d3a
Binary files /dev/null and b/tests/bcs/example_shape_side/gold/larger_bias_out.e differ
diff --git a/tests/bcs/example_shape_side/tests b/tests/bcs/example_shape_side/tests
new file mode 100644
index 0000000..ced6375
--- /dev/null
+++ b/tests/bcs/example_shape_side/tests
@@ -0,0 +1,23 @@
+[Tests]
+  [floating_electrode]
+    type = Exodiff
+    input = 'example_shape_side.i'
+    exodiff = 'example_shape_side_out.e'
+    requirement = 'The system shall float a boundary potential using nonlocal surface integrals of a coupled species supplied by shape side user objects.'
+  []
+  [larger_bias]
+    type = Exodiff
+    input = 'example_shape_side.i'
+    cli_args = 'Vb=5
+                Outputs/file_base=larger_bias_out'
+    exodiff = 'larger_bias_out.e'
+    requirement = 'The system shall shift the floated boundary potential with the applied bias.'
+  []
+  [jacobian]
+    type = PetscJacobianTester
+    input = 'example_shape_side.i'
+    run_sim = true
+    cli_args = 'Outputs/exodus=false'
+    requirement = 'The system shall provide the nonlocal Jacobian contributions of the shape side boundary condition.'
+  []
+[]
diff --git a/tests/bcs/postprocessor_dirichlet/gold/penalty_out.e b/tests/bcs/postprocessor_dirichlet/gold/penalty_out.e
new file mode 100644
index 0000000..b701a18
Binary files /dev/null and b/tests/bcs/postprocessor_dirichlet/gold/penalty_out.e differ
diff --git a/tests/bcs/postprocessor_dirichlet/gold/postprocessor_dirichlet_out.e b/tests/bcs/postprocessor_dirichlet/gold/postprocessor_dirichlet_out.e
new file mode 100644
index 0000000..4d99528
Binary files /dev/null and b/tests/bcs/postprocessor_dirichlet/gold/postprocessor_dirichlet_out.e differ
diff --git a/tests/bcs/postprocessor_dirichlet/postprocessor_dirichlet.i b/tests/bcs/postprocessor_dirichlet/postprocessor_dirichlet.i
new file mode 100644
index 0000000..6f71221
--- /dev/null
+++ b/tests/bcs/postprocessor_dirichlet/postprocessor_dirichlet.i
@@ -0,0 +1,69 @@
+# Conducting rod whose left end is held at a value derived from a postprocessor
+# (scale*sensor + offset = 2) and whose right end is held at 0, giving a linear
+# profile. The tests select between the two postprocessor-driven Dirichlet BCs.
+
+[Mesh]
+  [rod]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 4
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [T]
+  []
+[]
+
+[Kernels]
+  [conduction]
+    type = Diffusion
+    variable = T
+  []
+[]
+
+[BCs]
+  [left_flexible]
+    type = FlexiblePostprocessorDirichletBC
+    variable = T
+    boundary = left
+    postprocessor = sensor
+    scale = 0.5
+    offset = 1
+  []
+  [left_penalty]
+    type = PostprocessorPenaltyDirichletBC
+    variable = T
+    boundary = left
+    postprocessor = sensor
+    penalty = 1e8
+    enable = false
+  []
+  [right]
+    type = DirichletBC
+    variable = T
+    boundary = right
+    value = 0
+  []
+[]
+
+[Postprocessors]
+  [sensor]
+    type = FunctionValuePostprocessor
+    function = 2
+    execute_on = 'initial'
+    outputs = none
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/bcs/postprocessor_dirichlet/tests b/tests/bcs/postprocessor_dirichlet/tests
new file mode 100644
index 0000000..710fffd
--- /dev/null
+++ b/tests/bcs/postprocessor_dirichlet/tests
@@ -0,0 +1,17 @@
+[Tests]
+  [flexible]
+    type = Exodiff
+    input = 'postprocessor_dirichlet.i'
+    exodiff = 'postprocessor_dirichlet_out.e'
+    requirement = 'The system shall strongly impose a Dirichlet value taken from a scaled and offset postprocessor.'
+  []
+  [penalty]
+    type = Exodiff
+    input = 'postprocessor_dirichlet.i'
+    cli_args = 'BCs/left_flexible/enable=false
+                BCs/left_penalty/enable=true
+                Outputs/file_base=penalty_out'
+    exodiff = 'penalty_out.e'
+    requirement = 'The system shall weakly impose a postprocessor Dirichlet value through a penalty boundary condition.'
+  []
+[]
diff --git a/tests/bcs/robin/gold/faster_flow_out.e b/tests/bcs/robin/gold/faster_flow_out.e
new file mode 100644
index 0000000..e788810
Binary files /dev/null and b/tests/bcs/robin/gold/faster_flow_out.e differ
diff --git a/tests/bcs/robin/gold/robin_out.e b/tests/bcs/robin/gold/robin_out.e
new file mode 100644
index 0000000..8bbdaf4
Binary files /dev/null and b/tests/bcs/robin/gold/robin_out.e differ
diff --git a/tests/bcs/robin/robin.i b/tests/bcs/robin/robin.i
new file mode 100644
index 0000000..c836a23
--- /dev/null
+++ b/tests/bcs/robin/robin.i
@@ -0,0 +1,58 @@
+# Advection-diffusion of u down a channel with a Dirichlet inlet and an open
+# outlet. The RobinBC lets u leave at a rate proportional to its local value,
+# so the profile decays toward the outlet instead of staying flat.
+
+v = 1
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    initial_condition = 1
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = Diffusion
+    variable = u
+  []
+  [advection]
+    type = ConservativeAdvection
+    variable = u
+    velocity_variable = '${v} 0 0'
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = DirichletBC
+    variable = u
+    boundary = left
+    value = 1
+  []
+  [outlet]
+    type = RobinBC
+    variable = u
+    boundary = right
+    velocity = ${v}
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/bcs/robin/tests b/tests/bcs/robin/tests
new file mode 100644
index 0000000..783d381
--- /dev/null
+++ b/tests/bcs/robin/tests
@@ -0,0 +1,16 @@
+[Tests]
+  [robin]
+    type = Exodiff
+    input = 'robin.i'
+    exodiff = 'robin_out.e'
+    requirement = 'The system shall apply a Robin outflow boundary condition proportional to the local solution value.'
+  []
+  [faster_flow]
+    type = Exodiff
+    input = 'robin.i'
+    cli_args = 'v=5
+                Outputs/file_base=faster_flow_out'
+    exodiff = 'faster_flow_out.e'
+    requirement = 'The system shall scale the Robin outflow boundary condition with the specified velocity.'
+  []
+[]
diff --git a/tests/dgkernels/dg_advection/dg_advection.i b/tests/dgkernels/dg_advection/dg_advection.i
new file mode 100644
index 0000000..5ab56a8
--- /dev/null
+++ b/tests/dgkernels/dg_advection/dg_advection.i
@@ -0,0 +1,109 @@
+# Upwind DG advection of u with a uniform source: div(v u) = q, u_in = 1.
+# The cell averages step up by q*h/v along the channel.
+# The ...Temperature... flavours scale every term by rho*cp.
+
+v = 1
+q = 2
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    order = CONSTANT
+    family = MONOMIAL
+    initial_condition = 1
+  []
+[]
+
+[Functions]
+  [vel_x]
+    type = ParsedFunction
+    expression = '${v}'
+  []
+  [zero]
+    type = ParsedFunction
+    expression = '0'
+  []
+[]
+
+[Kernels]
+  [source]
+    type = BodyForce
+    variable = u
+    value = ${q}
+  []
+[]
+
+[DGKernels]
+  [function_advection]
+    type = DGFunctionConvection
+    variable = u
+    vel_x_func = vel_x
+    vel_y_func = zero
+    vel_z_func = zero
+  []
+  [advection]
+    type = DGTemperatureAdvection
+    variable = u
+    velocity = '${v} 0 0'
+    enable = false
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = InflowBC
+    variable = u
+    boundary = left
+    uu = ${v}
+    inlet_conc = 1
+  []
+  [outlet]
+    type = OutflowBC
+    variable = u
+    boundary = right
+    velocity = '${v} 0 0'
+  []
+[]
+
+[Materials]
+  [rho]
+    type = DerivativeParsedMaterial
+    property_name = rho
+    coupled_variables = 'u'
+    expression = '1'
+    derivative_order = 1
+  []
+  [cp]
+    type = DerivativeParsedMaterial
+    property_name = cp
+    coupled_variables = 'u'
+    expression = '1'
+    derivative_order = 1
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/dgkernels/dg_advection/gold/dg_advection_out.e b/tests/dgkernels/dg_advection/gold/dg_advection_out.e
new file mode 100644
index 0000000..ef95e54
Binary files /dev/null and b/tests/dgkernels/dg_advection/gold/dg_advection_out.e differ
diff --git a/tests/dgkernels/dg_advection/gold/function_temperature_out.e b/tests/dgkernels/dg_advection/gold/function_temperature_out.e
new file mode 100644
index 0000000..dd6773e
Binary files /dev/null and b/tests/dgkernels/dg_advection/gold/function_temperature_out.e differ
diff --git a/tests/dgkernels/dg_advection/gold/temperature_out.e b/tests/dgkernels/dg_advection/gold/temperature_out.e
new file mode 100644
index 0000000..32dd622
Binary files /dev/null and b/tests/dgkernels/dg_advection/gold/temperature_out.e differ
diff --git a/tests/dgkernels/dg_advection/tests b/tests/dgkernels/dg_advection/tests
new file mode 100644
index 0000000..a9747d9
--- /dev/null
+++ b/tests/dgkernels/dg_advection/tests
@@ -0,0 +1,33 @@
+[Tests]
+  [function_convection]
+    type = Exodiff
+    input = 'dg_advection.i'
+    exodiff = 'dg_advection_out.e'
+    requirement = 'The system shall advect a scalar with upwind discontinuous Galerkin fluxes using a function velocity.'
+  []
+  [function_temperature_advection]
+    type = Exodiff
+    input = 'dg_advection.i'
+    cli_args = 'DGKernels/function_advection/type=DGFunctionTemperatureAdvection
+                BCs/inlet/type=TemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=2
+                Materials/cp/expression=3
+                Outputs/file_base=function_temperature_out'
+    exodiff = 'function_temperature_out.e'
+    requirement = 'The system shall advect enthalpy with upwind discontinuous Galerkin fluxes using a function velocity scaled by density and specific heat.'
+  []
+  [temperature_advection]
+    type = Exodiff
+    input = 'dg_advection.i'
+    cli_args = 'DGKernels/function_advection/enable=false
+                DGKernels/advection/enable=true
+                BCs/inlet/type=TemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=2
+                Materials/cp/expression=3
+                Outputs/file_base=temperature_out'
+    exodiff = 'temperature_out.e'
+    requirement = 'The system shall advect enthalpy with upwind discontinuous Galerkin fluxes using a constant velocity vector.'
+  []
+[]
diff --git a/tests/dgkernels/dg_coupled_advection/dg_coupled_advection.i b/tests/dgkernels/dg_coupled_advection/dg_coupled_advection.i
new file mode 100644
index 0000000..8f575dd
--- /dev/null
+++ b/tests/dgkernels/dg_coupled_advection/dg_coupled_advection.i
@@ -0,0 +1,79 @@
+# Upwind DG advection of u with a uniform source, where the advecting velocity
+# is supplied as a variable field rather than a constant. The cell averages step
+# up by q*h/v along the channel.
+
+v = 1
+q = 2
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    order = CONSTANT
+    family = MONOMIAL
+    initial_condition = 1
+  []
+[]
+
+[AuxVariables]
+  [vel_x]
+    initial_condition = ${v}
+  []
+[]
+
+[Kernels]
+  [source]
+    type = BodyForce
+    variable = u
+    value = ${q}
+  []
+[]
+
+[DGKernels]
+  [advection]
+    type = DGCoupledAdvection
+    variable = u
+    uvel = vel_x
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = InflowBC
+    variable = u
+    boundary = left
+    uu = ${v}
+    inlet_conc = 1
+  []
+  [outlet]
+    type = OutflowBC
+    variable = u
+    boundary = right
+    velocity = '${v} 0 0'
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/dgkernels/dg_coupled_advection/gold/dg_coupled_advection_out.e b/tests/dgkernels/dg_coupled_advection/gold/dg_coupled_advection_out.e
new file mode 100644
index 0000000..61d0966
Binary files /dev/null and b/tests/dgkernels/dg_coupled_advection/gold/dg_coupled_advection_out.e differ
diff --git a/tests/dgkernels/dg_coupled_advection/gold/faster_flow_out.e b/tests/dgkernels/dg_coupled_advection/gold/faster_flow_out.e
new file mode 100644
index 0000000..e33e6b2
Binary files /dev/null and b/tests/dgkernels/dg_coupled_advection/gold/faster_flow_out.e differ
diff --git a/tests/dgkernels/dg_coupled_advection/tests b/tests/dgkernels/dg_coupled_advection/tests
new file mode 100644
index 0000000..849759f
--- /dev/null
+++ b/tests/dgkernels/dg_coupled_advection/tests
@@ -0,0 +1,16 @@
+[Tests]
+  [dg_coupled_advection]
+    type = Exodiff
+    input = 'dg_coupled_advection.i'
+    exodiff = 'dg_coupled_advection_out.e'
+    requirement = 'The system shall advect a scalar with upwind discontinuous Galerkin fluxes using a velocity supplied as a variable field.'
+  []
+  [faster_flow]
+    type = Exodiff
+    input = 'dg_coupled_advection.i'
+    cli_args = 'v=2
+                Outputs/file_base=faster_flow_out'
+    exodiff = 'faster_flow_out.e'
+    requirement = 'The system shall halve the concentration step across each cell when the coupled advecting velocity doubles.'
+  []
+[]
diff --git a/tests/dgkernels/dg_diffusion_pp_dirichlet/dg_diffusion_pp_dirichlet.i b/tests/dgkernels/dg_diffusion_pp_dirichlet/dg_diffusion_pp_dirichlet.i
new file mode 100644
index 0000000..252e763
--- /dev/null
+++ b/tests/dgkernels/dg_diffusion_pp_dirichlet/dg_diffusion_pp_dirichlet.i
@@ -0,0 +1,90 @@
+# Discontinuous Galerkin conduction with both ends set from a sensor
+# postprocessor: the hot end at scale*sensor + offset = 1.5 and the cold end at
+# 0. The steady profile is the straight line between them.
+
+[Mesh]
+  [rod]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 4
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [T]
+    order = FIRST
+    family = MONOMIAL
+  []
+[]
+
+[Kernels]
+  [conduction]
+    type = MatDiffusion
+    variable = T
+    diffusivity = k
+  []
+[]
+
+[DGKernels]
+  [dg_conduction]
+    type = DGDiffusion
+    variable = T
+    diff = k
+    sigma = 6
+    epsilon = -1
+  []
+[]
+
+[BCs]
+  [hot]
+    type = DGDiffusionPostprocessorDirichletBC
+    variable = T
+    boundary = left
+    postprocessor = sensor
+    scale = 1
+    offset = -0.5
+    D_name = k
+    sigma = 6
+    epsilon = -1
+  []
+  [cold]
+    type = DGDiffusionPostprocessorDirichletBC
+    variable = T
+    boundary = right
+    postprocessor = sensor
+    scale = 0
+    offset = 0
+    D_name = k
+    sigma = 6
+    epsilon = -1
+  []
+[]
+
+[Materials]
+  [conductivity]
+    type = GenericConstantMaterial
+    prop_names = 'k'
+    prop_values = '1'
+  []
+[]
+
+[Postprocessors]
+  [sensor]
+    type = FunctionValuePostprocessor
+    function = 2
+    execute_on = 'initial'
+    outputs = none
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/dgkernels/dg_diffusion_pp_dirichlet/gold/dg_diffusion_pp_dirichlet_out.e b/tests/dgkernels/dg_diffusion_pp_dirichlet/gold/dg_diffusion_pp_dirichlet_out.e
new file mode 100644
index 0000000..5bf4295
Binary files /dev/null and b/tests/dgkernels/dg_diffusion_pp_dirichlet/gold/dg_diffusion_pp_dirichlet_out.e differ
diff --git a/tests/dgkernels/dg_diffusion_pp_dirichlet/gold/hotter_inlet_out.e b/tests/dgkernels/dg_diffusion_pp_dirichlet/gold/hotter_inlet_out.e
new file mode 100644
index 0000000..23e6f2b
Binary files /dev/null and b/tests/dgkernels/dg_diffusion_pp_dirichlet/gold/hotter_inlet_out.e differ
diff --git a/tests/dgkernels/dg_diffusion_pp_dirichlet/tests b/tests/dgkernels/dg_diffusion_pp_dirichlet/tests
new file mode 100644
index 0000000..468514e
--- /dev/null
+++ b/tests/dgkernels/dg_diffusion_pp_dirichlet/tests
@@ -0,0 +1,16 @@
+[Tests]
+  [dg_diffusion_pp_dirichlet]
+    type = Exodiff
+    input = 'dg_diffusion_pp_dirichlet.i'
+    exodiff = 'dg_diffusion_pp_dirichlet_out.e'
+    requirement = 'The system shall weakly impose a Dirichlet value from a postprocessor on a discontinuous Galerkin diffusion problem.'
+  []
+  [hotter_inlet]
+    type = Exodiff
+    input = 'dg_diffusion_pp_dirichlet.i'
+    cli_args = 'BCs/hot/offset=0.5
+                Outputs/file_base=hotter_inlet_out'
+    exodiff = 'hotter_inlet_out.e'
+    requirement = 'The system shall steepen the discontinuous Galerkin diffusion profile when the postprocessor boundary offset increases.'
+  []
+[]
diff --git a/tests/interfacekernels/inter_temperature_advection/gold/heat_removal_out.e b/tests/interfacekernels/inter_temperature_advection/gold/heat_removal_out.e
new file mode 100644
index 0000000..e02cd32
Binary files /dev/null and b/tests/interfacekernels/inter_temperature_advection/gold/heat_removal_out.e differ
diff --git a/tests/interfacekernels/inter_temperature_advection/gold/inter_temperature_advection_out.e b/tests/interfacekernels/inter_temperature_advection/gold/inter_temperature_advection_out.e
new file mode 100644
index 0000000..0fad5ad
Binary files /dev/null and b/tests/interfacekernels/inter_temperature_advection/gold/inter_temperature_advection_out.e differ
diff --git a/tests/interfacekernels/inter_temperature_advection/inter_temperature_advection.i b/tests/interfacekernels/inter_temperature_advection/inter_temperature_advection.i
new file mode 100644
index 0000000..3a6e6a6
--- /dev/null
+++ b/tests/interfacekernels/inter_temperature_advection/inter_temperature_advection.i
@@ -0,0 +1,112 @@
+# Enthalpy advected along a channel split into two blocks by a lumped heat
+# exchanger on the interface. The exchanger deposits -heat_source per unit area
+# into the downstream stream, so heat_source = -1 warms it from 2 K to 3 K.
+
+v = 1
+T_in = 2
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+  [downstream]
+    type = SubdomainBoundingBoxGenerator
+    input = channel
+    block_id = 1
+    bottom_left = '0.5 0 0'
+    top_right = '1 0 0'
+  []
+  [exchanger]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = downstream
+    primary_block = 0
+    paired_block = 1
+    new_boundary = 'exchanger'
+  []
+[]
+
+[Variables]
+  [T_upstream]
+    order = CONSTANT
+    family = MONOMIAL
+    block = 0
+    initial_condition = ${T_in}
+  []
+  [T_downstream]
+    order = CONSTANT
+    family = MONOMIAL
+    block = 1
+    initial_condition = ${T_in}
+  []
+[]
+
+[DGKernels]
+  [upstream_advection]
+    type = DGTemperatureAdvection
+    variable = T_upstream
+    velocity = '${v} 0 0'
+    block = 0
+  []
+  [downstream_advection]
+    type = DGTemperatureAdvection
+    variable = T_downstream
+    velocity = '${v} 0 0'
+    block = 1
+  []
+[]
+
+[InterfaceKernels]
+  [exchanger]
+    type = InterTemperatureAdvection
+    variable = T_upstream
+    neighbor_var = T_downstream
+    boundary = 'exchanger'
+    u_val = ${v}
+    heat_source = -1
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = TemperatureInflowBC
+    variable = T_upstream
+    boundary = left
+    uu = ${v}
+    inlet_conc = ${T_in}
+  []
+  [outlet]
+    type = TemperatureOutflowBC
+    variable = T_downstream
+    boundary = right
+    velocity = '${v} 0 0'
+  []
+[]
+
+[Materials]
+  [properties]
+    type = GenericConstantMaterial
+    prop_names = 'rho cp'
+    prop_values = '1 1'
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/interfacekernels/inter_temperature_advection/tests b/tests/interfacekernels/inter_temperature_advection/tests
new file mode 100644
index 0000000..e8d7308
--- /dev/null
+++ b/tests/interfacekernels/inter_temperature_advection/tests
@@ -0,0 +1,16 @@
+[Tests]
+  [heat_addition]
+    type = Exodiff
+    input = 'inter_temperature_advection.i'
+    exodiff = 'inter_temperature_advection_out.e'
+    requirement = 'The system shall exchange advected enthalpy across a block interface, adding heat to the downstream stream.'
+  []
+  [heat_removal]
+    type = Exodiff
+    input = 'inter_temperature_advection.i'
+    cli_args = 'InterfaceKernels/exchanger/heat_source=1
+                Outputs/file_base=heat_removal_out'
+    exodiff = 'heat_removal_out.e'
+    requirement = 'The system shall remove heat from the downstream stream when the interface heat source sign is reversed.'
+  []
+[]
diff --git a/tests/kernels/conservative_advection/conservative_advection.i b/tests/kernels/conservative_advection/conservative_advection.i
new file mode 100644
index 0000000..ec87471
--- /dev/null
+++ b/tests/kernels/conservative_advection/conservative_advection.i
@@ -0,0 +1,107 @@
+# Steady advection of u along a channel with a uniform source: div(v u) = q.
+# With v = 1, q = 2 and u(0) = 1 the exact solution is the line u = 1 + 2 x.
+# The ...Temperature... objects scale every term by rho*cp.
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    initial_condition = 1
+  []
+[]
+
+[Kernels]
+  [advection]
+    type = ConservativeAdvection
+    variable = u
+    velocity_variable = '1 0 0'
+  []
+  [ctrl_advection]
+    type = CtrlConservativeAdvection
+    variable = u
+    u_val = 1
+    v_val = 0
+    w_val = 0
+    enable = false
+  []
+  [source]
+    type = BodyForce
+    variable = u
+    value = 2
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = InflowBC
+    variable = u
+    boundary = left
+    uu = 1
+    inlet_conc = 1
+  []
+  [pp_inlet]
+    type = PostprocessorInflowBC
+    variable = u
+    boundary = left
+    uu = 1
+    postprocessor = inlet_value
+    enable = false
+  []
+  [outlet]
+    type = OutflowBC
+    variable = u
+    boundary = right
+    velocity = '1 0 0'
+  []
+[]
+
+[Materials]
+  [rho]
+    type = DerivativeParsedMaterial
+    property_name = rho
+    coupled_variables = 'u'
+    expression = '1'
+    derivative_order = 1
+  []
+  [cp]
+    type = DerivativeParsedMaterial
+    property_name = cp
+    coupled_variables = 'u'
+    expression = '1'
+    derivative_order = 1
+  []
+[]
+
+[Postprocessors]
+  [inlet_value]
+    type = FunctionValuePostprocessor
+    function = 1
+    execute_on = 'initial'
+    outputs = none
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/kernels/conservative_advection/gold/conservative_advection_out.e b/tests/kernels/conservative_advection/gold/conservative_advection_out.e
new file mode 100644
index 0000000..7b3e57b
Binary files /dev/null and b/tests/kernels/conservative_advection/gold/conservative_advection_out.e differ
diff --git a/tests/kernels/conservative_advection/gold/ctrl_out.e b/tests/kernels/conservative_advection/gold/ctrl_out.e
new file mode 100644
index 0000000..0e44565
Binary files /dev/null and b/tests/kernels/conservative_advection/gold/ctrl_out.e differ
diff --git a/tests/kernels/conservative_advection/gold/ctrl_temperature_out.e b/tests/kernels/conservative_advection/gold/ctrl_temperature_out.e
new file mode 100644
index 0000000..f4f9b07
Binary files /dev/null and b/tests/kernels/conservative_advection/gold/ctrl_temperature_out.e differ
diff --git a/tests/kernels/conservative_advection/gold/pp_inlet_out.e b/tests/kernels/conservative_advection/gold/pp_inlet_out.e
new file mode 100644
index 0000000..dffc401
Binary files /dev/null and b/tests/kernels/conservative_advection/gold/pp_inlet_out.e differ
diff --git a/tests/kernels/conservative_advection/gold/pp_temperature_inlet_out.e b/tests/kernels/conservative_advection/gold/pp_temperature_inlet_out.e
new file mode 100644
index 0000000..550bf1f
Binary files /dev/null and b/tests/kernels/conservative_advection/gold/pp_temperature_inlet_out.e differ
diff --git a/tests/kernels/conservative_advection/gold/temperature_out.e b/tests/kernels/conservative_advection/gold/temperature_out.e
new file mode 100644
index 0000000..5ab8423
Binary files /dev/null and b/tests/kernels/conservative_advection/gold/temperature_out.e differ
diff --git a/tests/kernels/conservative_advection/tests b/tests/kernels/conservative_advection/tests
new file mode 100644
index 0000000..118e23e
--- /dev/null
+++ b/tests/kernels/conservative_advection/tests
@@ -0,0 +1,90 @@
+[Tests]
+  [scalar]
+    type = Exodiff
+    input = 'conservative_advection.i'
+    exodiff = 'conservative_advection_out.e'
+    requirement = 'The system shall advect a scalar in conservative form with an inflow and an outflow boundary condition.'
+  []
+  [ctrl]
+    type = Exodiff
+    input = 'conservative_advection.i'
+    cli_args = 'Kernels/advection/enable=false
+                Kernels/ctrl_advection/enable=true
+                Outputs/file_base=ctrl_out'
+    exodiff = 'ctrl_out.e'
+    requirement = 'The system shall advect a scalar in conservative form using a controllable constant velocity.'
+  []
+  [postprocessor_inlet]
+    type = Exodiff
+    input = 'conservative_advection.i'
+    cli_args = 'BCs/inlet/enable=false
+                BCs/pp_inlet/enable=true
+                Outputs/file_base=pp_inlet_out'
+    exodiff = 'pp_inlet_out.e'
+    requirement = 'The system shall set the inflow concentration of an advected scalar from a postprocessor value.'
+  []
+  [temperature]
+    type = Exodiff
+    input = 'conservative_advection.i'
+    cli_args = 'Kernels/advection/type=ConservativeTemperatureAdvection
+                BCs/inlet/type=TemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=2
+                Materials/cp/expression=3
+                Outputs/file_base=temperature_out'
+    exodiff = 'temperature_out.e'
+    requirement = 'The system shall advect enthalpy in conservative form, scaling the conserved advection, inflow, and outflow terms by the density and specific heat.'
+  []
+  [ctrl_temperature]
+    type = Exodiff
+    input = 'conservative_advection.i'
+    cli_args = 'Kernels/advection/enable=false
+                Kernels/ctrl_advection/enable=true
+                Kernels/ctrl_advection/type=CtrlConservativeTemperatureAdvection
+                BCs/inlet/type=TemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=2
+                Materials/cp/expression=3
+                Outputs/file_base=ctrl_temperature_out'
+    exodiff = 'ctrl_temperature_out.e'
+    requirement = 'The system shall advect enthalpy in conservative form using a controllable constant velocity.'
+  []
+  [postprocessor_temperature_inlet]
+    type = Exodiff
+    input = 'conservative_advection.i'
+    cli_args = 'Kernels/advection/type=ConservativeTemperatureAdvection
+                BCs/inlet/enable=false
+                BCs/pp_inlet/enable=true
+                BCs/pp_inlet/type=PostprocessorTemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=2
+                Materials/cp/expression=3
+                Outputs/file_base=pp_temperature_inlet_out'
+    exodiff = 'pp_temperature_inlet_out.e'
+    requirement = 'The system shall set the enthalpy inflow of an advected temperature from a postprocessor value.'
+  []
+  [temperature_jacobian]
+    type = PetscJacobianTester
+    input = 'conservative_advection.i'
+    run_sim = true
+    cli_args = 'Kernels/advection/type=ConservativeTemperatureAdvection
+                BCs/inlet/type=TemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=1+0.1*u
+                Materials/cp/expression=1+0.05*u'
+    requirement = 'The system shall provide an exact Jacobian for conservative enthalpy advection with temperature-dependent density and specific heat.'
+  []
+  [ctrl_temperature_jacobian]
+    type = PetscJacobianTester
+    input = 'conservative_advection.i'
+    run_sim = true
+    cli_args = 'Kernels/advection/enable=false
+                Kernels/ctrl_advection/enable=true
+                Kernels/ctrl_advection/type=CtrlConservativeTemperatureAdvection
+                BCs/inlet/type=TemperatureInflowBC
+                BCs/outlet/type=TemperatureOutflowBC
+                Materials/rho/expression=1+0.1*u
+                Materials/cp/expression=1+0.05*u'
+    requirement = 'The system shall provide an exact Jacobian for controllable conservative enthalpy advection with temperature-dependent properties.'
+  []
+[]
diff --git a/tests/kernels/mat_ins_temperature_time_derivative/gold/heated_slab_out.csv b/tests/kernels/mat_ins_temperature_time_derivative/gold/heated_slab_out.csv
new file mode 100644
index 0000000..a1a8630
--- /dev/null
+++ b/tests/kernels/mat_ins_temperature_time_derivative/gold/heated_slab_out.csv
@@ -0,0 +1,4 @@
+time,T_avg
+0,0
+0.5,0.5
+1,1
diff --git a/tests/kernels/mat_ins_temperature_time_derivative/heated_slab.i b/tests/kernels/mat_ins_temperature_time_derivative/heated_slab.i
new file mode 100644
index 0000000..4a4f815
--- /dev/null
+++ b/tests/kernels/mat_ins_temperature_time_derivative/heated_slab.i
@@ -0,0 +1,73 @@
+# Uniformly heated stagnant slab: rho*cp*dT/dt = q.
+# With rho = 2, cp = 3 and q = 6 the slab warms at exactly 1 K/s from T = 0.
+
+[Mesh]
+  [slab]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 1
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [T]
+    initial_condition = 0
+  []
+[]
+
+[Kernels]
+  [time]
+    type = MatINSTemperatureTimeDerivative
+    variable = T
+  []
+  [heating]
+    type = BodyForce
+    variable = T
+    value = 6
+  []
+[]
+
+[Materials]
+  [rho]
+    type = DerivativeParsedMaterial
+    property_name = rho
+    coupled_variables = 'T'
+    expression = '2'
+    derivative_order = 1
+  []
+  [cp]
+    type = DerivativeParsedMaterial
+    property_name = cp
+    coupled_variables = 'T'
+    expression = '3'
+    derivative_order = 1
+  []
+[]
+
+[Postprocessors]
+  [T_avg]
+    type = ElementAverageValue
+    variable = T
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+  dt = 0.5
+  num_steps = 2
+[]
+
+[Outputs]
+  csv = true
+[]
diff --git a/tests/kernels/mat_ins_temperature_time_derivative/tests b/tests/kernels/mat_ins_temperature_time_derivative/tests
new file mode 100644
index 0000000..daef4b2
--- /dev/null
+++ b/tests/kernels/mat_ins_temperature_time_derivative/tests
@@ -0,0 +1,17 @@
+[Tests]
+  [constant_properties]
+    type = CSVDiff
+    input = 'heated_slab.i'
+    csvdiff = 'heated_slab_out.csv'
+    requirement = 'The system shall integrate the transient enthalpy accumulation term with constant density and specific heat.'
+  []
+  [variable_properties_jacobian]
+    type = PetscJacobianTester
+    input = 'heated_slab.i'
+    run_sim = true
+    cli_args = 'Materials/rho/expression=2*(1+0.5*T)
+                Materials/cp/expression=3*(1+0.25*T)
+                Outputs/csv=false'
+    requirement = 'The system shall provide an exact Jacobian for the transient enthalpy accumulation term with temperature-dependent properties.'
+  []
+[]
diff --git a/tests/kernels/nonconservative_advection/gold/nonconservative_advection_out.e b/tests/kernels/nonconservative_advection/gold/nonconservative_advection_out.e
new file mode 100644
index 0000000..3c9e4f8
Binary files /dev/null and b/tests/kernels/nonconservative_advection/gold/nonconservative_advection_out.e differ
diff --git a/tests/kernels/nonconservative_advection/gold/reverse_out.e b/tests/kernels/nonconservative_advection/gold/reverse_out.e
new file mode 100644
index 0000000..f593e83
Binary files /dev/null and b/tests/kernels/nonconservative_advection/gold/reverse_out.e differ
diff --git a/tests/kernels/nonconservative_advection/nonconservative_advection.i b/tests/kernels/nonconservative_advection/nonconservative_advection.i
new file mode 100644
index 0000000..57791a4
--- /dev/null
+++ b/tests/kernels/nonconservative_advection/nonconservative_advection.i
@@ -0,0 +1,52 @@
+# Steady advection of u in non-conservative form with a uniform source: v u' = q.
+# With v = 1, q = 2 and u(0) = 1 the exact solution is the line u = 1 + 2 x.
+
+v = 1
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    initial_condition = 1
+  []
+[]
+
+[Kernels]
+  [advection]
+    type = NonConservativeAdvection
+    variable = u
+    velocity = '${v} 0 0'
+  []
+  [source]
+    type = BodyForce
+    variable = u
+    value = 2
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = DirichletBC
+    variable = u
+    boundary = left
+    value = 1
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/kernels/nonconservative_advection/tests b/tests/kernels/nonconservative_advection/tests
new file mode 100644
index 0000000..a4a76c8
--- /dev/null
+++ b/tests/kernels/nonconservative_advection/tests
@@ -0,0 +1,17 @@
+[Tests]
+  [forward]
+    type = Exodiff
+    input = 'nonconservative_advection.i'
+    exodiff = 'nonconservative_advection_out.e'
+    requirement = 'The system shall advect a scalar in non-conservative form driven from the inlet.'
+  []
+  [reverse]
+    type = Exodiff
+    input = 'nonconservative_advection.i'
+    cli_args = 'v=-1
+                BCs/inlet/boundary=right
+                Outputs/file_base=reverse_out'
+    exodiff = 'reverse_out.e'
+    requirement = 'The system shall advect a scalar in non-conservative form for a reversed velocity driven from the opposite boundary.'
+  []
+[]
diff --git a/tests/kernels/potential_advection/gold/negative_charge_out.e b/tests/kernels/potential_advection/gold/negative_charge_out.e
new file mode 100644
index 0000000..7301b8c
Binary files /dev/null and b/tests/kernels/potential_advection/gold/negative_charge_out.e differ
diff --git a/tests/kernels/potential_advection/gold/potential_advection_out.e b/tests/kernels/potential_advection/gold/potential_advection_out.e
new file mode 100644
index 0000000..dbcca9d
Binary files /dev/null and b/tests/kernels/potential_advection/gold/potential_advection_out.e differ
diff --git a/tests/kernels/potential_advection/potential_advection.i b/tests/kernels/potential_advection/potential_advection.i
new file mode 100644
index 0000000..218a3a6
--- /dev/null
+++ b/tests/kernels/potential_advection/potential_advection.i
@@ -0,0 +1,70 @@
+# Charged scalar u under diffusion and drift in a prescribed potential.
+# The potential falls linearly from 1 to 0, so the field points in +x and
+# positive charges drift toward the right electrode while negatives drift left.
+
+[Mesh]
+  [gap]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    initial_condition = 0.5
+  []
+[]
+
+[AuxVariables]
+  [potential]
+  []
+[]
+
+[AuxKernels]
+  [potential]
+    type = FunctionAux
+    variable = potential
+    function = '1 - x'
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = Diffusion
+    variable = u
+  []
+  [drift]
+    type = PotentialAdvection
+    variable = u
+    potential = potential
+    positive_charge = true
+  []
+[]
+
+[BCs]
+  [left]
+    type = DirichletBC
+    variable = u
+    boundary = left
+    value = 1
+  []
+  [right]
+    type = DirichletBC
+    variable = u
+    boundary = right
+    value = 0
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/kernels/potential_advection/tests b/tests/kernels/potential_advection/tests
new file mode 100644
index 0000000..0bf63a7
--- /dev/null
+++ b/tests/kernels/potential_advection/tests
@@ -0,0 +1,22 @@
+[Tests]
+  [positive_charge]
+    type = Exodiff
+    input = 'potential_advection.i'
+    exodiff = 'potential_advection_out.e'
+    requirement = 'The system shall drift a positively charged scalar down the gradient of a prescribed potential.'
+  []
+  [negative_charge]
+    type = Exodiff
+    input = 'potential_advection.i'
+    cli_args = 'Kernels/drift/positive_charge=false
+                Outputs/file_base=negative_charge_out'
+    exodiff = 'negative_charge_out.e'
+    requirement = 'The system shall drift a negatively charged scalar up the gradient of a prescribed potential.'
+  []
+  [jacobian]
+    type = PetscJacobianTester
+    input = 'potential_advection.i'
+    run_sim = true
+    requirement = 'The system shall provide an exact Jacobian for potential-driven drift of a charged scalar.'
+  []
+[]
diff --git a/tests/kernels/simple_diffusion/gold/simple_diffusion_out.e b/tests/kernels/simple_diffusion/gold/simple_diffusion_out.e
deleted file mode 100644
index 26b05ed..0000000
Binary files a/tests/kernels/simple_diffusion/gold/simple_diffusion_out.e and /dev/null differ
diff --git a/tests/kernels/simple_diffusion/simple_diffusion.i b/tests/kernels/simple_diffusion/simple_diffusion.i
deleted file mode 100644
index 4e46e0a..0000000
--- a/tests/kernels/simple_diffusion/simple_diffusion.i
+++ /dev/null
@@ -1,44 +0,0 @@
-[Mesh]
-  type = GeneratedMesh
-  dim = 2
-  nx = 10
-  ny = 10
-[]
-
-[Variables]
-  [./u]
-  [../]
-[]
-
-[Kernels]
-  [./diff]
-    type = Diffusion
-    variable = u
-  [../]
-[]
-
-[BCs]
-  [./left]
-    type = DirichletBC
-    variable = u
-    boundary = left
-    value = 0
-  [../]
-  [./right]
-    type = DirichletBC
-    variable = u
-    boundary = right
-    value = 1
-  [../]
-[]
-
-[Executioner]
-  type = Steady
-  solve_type = 'PJFNK'
-  petsc_options_iname = '-pc_type -pc_hypre_type'
-  petsc_options_value = 'hypre boomeramg'
-[]
-
-[Outputs]
-  exodus = true
-[]
diff --git a/tests/kernels/simple_diffusion/tests b/tests/kernels/simple_diffusion/tests
deleted file mode 100644
index e4545a4..0000000
--- a/tests/kernels/simple_diffusion/tests
+++ /dev/null
@@ -1,7 +0,0 @@
-[Tests]
-  [./test]
-    type = 'Exodiff'
-    input = 'simple_diffusion.i'
-    exodiff = 'simple_diffusion_out.e'
-  [../]
-[]
diff --git a/tests/kernels/velocity_function_advection/gold/temperature_out.e b/tests/kernels/velocity_function_advection/gold/temperature_out.e
new file mode 100644
index 0000000..ad1282a
Binary files /dev/null and b/tests/kernels/velocity_function_advection/gold/temperature_out.e differ
diff --git a/tests/kernels/velocity_function_advection/gold/velocity_function_advection_out.e b/tests/kernels/velocity_function_advection/gold/velocity_function_advection_out.e
new file mode 100644
index 0000000..358afb9
Binary files /dev/null and b/tests/kernels/velocity_function_advection/gold/velocity_function_advection_out.e differ
diff --git a/tests/kernels/velocity_function_advection/tests b/tests/kernels/velocity_function_advection/tests
new file mode 100644
index 0000000..a8301a2
--- /dev/null
+++ b/tests/kernels/velocity_function_advection/tests
@@ -0,0 +1,31 @@
+[Tests]
+  [scalar]
+    type = Exodiff
+    input = 'velocity_function_advection.i'
+    exodiff = 'velocity_function_advection_out.e'
+    requirement = 'The system shall advect a scalar in conservative form using a spatially varying velocity supplied as functions.'
+  []
+  [temperature]
+    type = Exodiff
+    input = 'velocity_function_advection.i'
+    cli_args = 'Kernels/advection/type=VelocityFunctionTemperatureAdvection
+                BCs/inlet/type=PostprocessorVelocityFunctionTemperatureInflowBC
+                BCs/outlet/type=VelocityFunctionTemperatureOutflowBC
+                Materials/rho/expression=2
+                Materials/cp/expression=3
+                Outputs/file_base=temperature_out'
+    exodiff = 'temperature_out.e'
+    requirement = 'The system shall advect enthalpy using a function velocity, scaling the advection, inflow, and outflow terms by the density and specific heat.'
+  []
+  [temperature_jacobian]
+    type = PetscJacobianTester
+    input = 'velocity_function_advection.i'
+    run_sim = true
+    cli_args = 'Kernels/advection/type=VelocityFunctionTemperatureAdvection
+                BCs/inlet/type=PostprocessorVelocityFunctionTemperatureInflowBC
+                BCs/outlet/type=VelocityFunctionTemperatureOutflowBC
+                Materials/rho/expression=1+0.1*u
+                Materials/cp/expression=1+0.05*u'
+    requirement = 'The system shall provide an exact Jacobian for function-velocity enthalpy advection with temperature-dependent properties.'
+  []
+[]
diff --git a/tests/kernels/velocity_function_advection/velocity_function_advection.i b/tests/kernels/velocity_function_advection/velocity_function_advection.i
new file mode 100644
index 0000000..0477fb3
--- /dev/null
+++ b/tests/kernels/velocity_function_advection/velocity_function_advection.i
@@ -0,0 +1,103 @@
+# Steady advection of u through an accelerating flow v(x) = 1 + x with no source.
+# The flux v u is constant, so u(x) = 1 / (1 + x) and u halves across the domain.
+# The ...Temperature... objects scale every term by rho*cp.
+
+[Mesh]
+  [channel]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 1
+  []
+[]
+
+[Variables]
+  [u]
+    initial_condition = 1
+  []
+[]
+
+[Functions]
+  [vel_x]
+    type = ParsedFunction
+    expression = '1 + x'
+  []
+  [zero]
+    type = ParsedFunction
+    expression = '0'
+  []
+[]
+
+[Kernels]
+  [advection]
+    type = VelocityFunctionConservativeAdvection
+    variable = u
+    vel_x_func = vel_x
+    vel_y_func = zero
+    vel_z_func = zero
+  []
+[]
+
+[BCs]
+  [inlet]
+    type = PostprocessorVelocityFunctionInflowBC
+    variable = u
+    boundary = left
+    vel_x_func = vel_x
+    vel_y_func = zero
+    vel_z_func = zero
+    postprocessor = inlet_value
+  []
+  [outlet]
+    type = VelocityFunctionOutflowBC
+    variable = u
+    boundary = right
+    vel_x_func = vel_x
+    vel_y_func = zero
+    vel_z_func = zero
+  []
+[]
+
+[Materials]
+  [rho]
+    type = DerivativeParsedMaterial
+    property_name = rho
+    coupled_variables = 'u'
+    expression = '1'
+    derivative_order = 1
+  []
+  [cp]
+    type = DerivativeParsedMaterial
+    property_name = cp
+    coupled_variables = 'u'
+    expression = '1'
+    derivative_order = 1
+  []
+[]
+
+[Postprocessors]
+  [inlet_value]
+    type = FunctionValuePostprocessor
+    function = 1
+    execute_on = 'initial'
+    outputs = none
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+[]
+
+[Outputs]
+  exodus = true
+[]