Merge pull request #28978 from GiudGiud/PR_aux_ik

Support auxiliary neighbor variables for FE IKs

framework/doc/content/source/timeintegrators/Ralston.md

@@ -35,7 +35,7 @@ All kernels except time-(derivative)-kernels should have the parameter `implicit
 time integrator.
 
 !alert warning
-ExplicitRK2-derived TimeIntegrators [ExplicitMidpoint.md], [Heun.md], [Ralston.md]) and other multistage
+ExplicitRK2-derived TimeIntegrators ([ExplicitMidpoint.md], [Heun.md], [Ralston.md]) and other multistage
 TimeIntegrators are known not to work with Materials/AuxKernels that accumulate 'state' and
 should be used with caution.
 

framework/doc/content/syntax/InterfaceKernels/index.md

@@ -23,6 +23,13 @@ Interface kernels can be used to provide any general flux condition at an interf
 
 For additional information about the interface kernel system, don't hesitate to contact the [MOOSE Discussion forum](https://github.com/idaholab/moose/discussions).
 
+## Multiple system support
+
+Using multiple nonlinear system with interface kernels is currently not supported.
+The only feature supported, which can at times suffice, notably when using
+[MultiApps](syntax/MultiApps/index.md) to couple equations, is to use an auxiliary
+variable as the `neighbor_var`.
+
 !syntax list /InterfaceKernels objects=True actions=False subsystems=False
 
 !syntax list /InterfaceKernels objects=False actions=False subsystems=True

framework/include/interfacekernels/InterfaceKernel.h

@@ -52,8 +52,7 @@ class InterfaceKernelTempl : public InterfaceKernelBase, public NeighborMooseVar
 
   /**
    * Using the passed DGJacobian type, selects the correct test function and trial function spaces
-   * and
-   * jacobian block, and then calls computeQpJacobian
+   * and jacobian block, and then calls computeQpJacobian
    */
   virtual void computeElemNeighJacobian(Moose::DGJacobianType type);
 
@@ -144,6 +143,10 @@ class InterfaceKernelTempl : public InterfaceKernelBase, public NeighborMooseVar
   /// Gradient of side neighbor shape function
   const TemplateVariableTestGradient & _grad_test_neighbor;
 
+  /// Whether the variable and the neighbor variables are part of the same system
+  /// (whether from two different nonlinear systems or a nonlinear and an auxiliary system)
+  const bool _same_system;
+
   /// Holds residual entries as they are accumulated by this InterfaceKernel
   /// This variable is temporarily reserved for RattleSnake
   DenseMatrix<Number> _local_kxx;

framework/src/interfacekernels/InterfaceKernel.C

@@ -13,6 +13,7 @@
 #include "Assembly.h"
 #include "MooseVariableFE.h"
 #include "SystemBase.h"
+#include "AuxiliarySystem.h"
 
 #include "libmesh/quadrature.h"
 
@@ -53,9 +54,11 @@ InterfaceKernelTempl<T>::InterfaceKernelTempl(const InputParameters & parameters
     _phi_neighbor(_assembly.phiFaceNeighbor(_neighbor_var)),
     _grad_phi_neighbor(_assembly.gradPhiFaceNeighbor(_neighbor_var)),
     _test_neighbor(_neighbor_var.phiFaceNeighbor()),
-    _grad_test_neighbor(_neighbor_var.gradPhiFaceNeighbor())
-
+    _grad_test_neighbor(_neighbor_var.gradPhiFaceNeighbor()),
+    _same_system(_var.sys().number() == _neighbor_var.sys().number())
 {
+  // Neighbor variable dependency is added by
+  // NeighborCoupleableMooseVariableDependencyIntermediateInterface
   addMooseVariableDependency(this->mooseVariable());
 
   if (!parameters.isParamValid("boundary"))
@@ -177,6 +180,7 @@ InterfaceKernelTempl<T>::computeElemNeighResidual(Moose::DGResidualType type)
 
   accumulateTaggedLocalResidual();
 
+  // To save the diagonal of the Jacobian
   if (_has_primary_residuals_saved_in && is_elem)
   {
     Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
@@ -216,7 +220,9 @@ InterfaceKernelTempl<T>::computeResidual()
   computeElemNeighResidual(Moose::Element);
 
   // Compute the residual for the neighbor
-  computeElemNeighResidual(Moose::Neighbor);
+  // This also prevents computing a residual if the neighbor variable is auxiliary
+  if (_same_system)
+    computeElemNeighResidual(Moose::Neighbor);
 }
 
 template <typename T>
@@ -263,6 +269,7 @@ InterfaceKernelTempl<T>::computeElemNeighJacobian(Moose::DGJacobianType type)
 
   accumulateTaggedLocalMatrix();
 
+  // To save the diagonal of the Jacobian
   if (_has_primary_jacobians_saved_in && type == Moose::ElementElement)
   {
     auto rows = _local_ke.m();
@@ -307,7 +314,8 @@ InterfaceKernelTempl<T>::computeJacobian()
   precalculateJacobian();
 
   computeElemNeighJacobian(Moose::ElementElement);
-  computeElemNeighJacobian(Moose::NeighborNeighbor);
+  if (_same_system)
+    computeElemNeighJacobian(Moose::NeighborNeighbor);
 }
 
 template <typename T>
@@ -369,7 +377,7 @@ InterfaceKernelTempl<T>::computeElementOffDiagJacobian(unsigned int jvar)
     computeElemNeighJacobian(Moose::ElementElement);
     is_jvar_not_interface_var = false;
   }
-  if (jvar == _neighbor_var.number())
+  if (jvar == _neighbor_var.number() && _same_system)
   {
     precalculateJacobian();
     computeElemNeighJacobian(Moose::ElementNeighbor);
@@ -401,6 +409,11 @@ InterfaceKernelTempl<T>::computeNeighborOffDiagJacobian(unsigned int jvar)
       !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
     return;
 
+  // We don't care about any contribution to the neighbor Jacobian rows if it's not in the system
+  // we are currently working with (the variable's system)
+  if (!_same_system)
+    return;
+
   bool is_jvar_not_interface_var = true;
   if (jvar == _var.number())
   {
@@ -446,8 +459,9 @@ InterfaceKernelTempl<T>::computeResidualAndJacobian()
     if (_var.number() == ivar)
       computeElementOffDiagJacobian(jvar);
 
-    if (_neighbor_var.number() == ivar)
-      computeNeighborOffDiagJacobian(jvar);
+    if (_same_system)
+      if (_neighbor_var.number() == ivar)
+        computeNeighborOffDiagJacobian(jvar);
   }
 }
 

framework/src/interfacekernels/InterfaceKernelBase.C

@@ -62,7 +62,7 @@ InterfaceKernelBase::validParams()
       "length as diag_save_in. This vector specifies whether the corresponding aux_var should "
       "save-in jacobian contributions from the primary ('p') or secondary side ('s').");
   params.addParamNamesToGroup("diag_save_in save_in save_in_var_side diag_save_in_var_side",
-                              "Advanced");
+                              "Residual and Jacobian debug output");
 
   // InterfaceKernels always need one layer of ghosting.
   params.addRelationshipManager("ElementSideNeighborLayers",
@@ -98,7 +98,6 @@ InterfaceKernelBase::InterfaceKernelBase(const InputParameters & parameters)
     _save_in_strings(parameters.get<std::vector<AuxVariableName>>("save_in")),
     _diag_save_in_var_side(parameters.get<MultiMooseEnum>("diag_save_in_var_side")),
     _diag_save_in_strings(parameters.get<std::vector<AuxVariableName>>("diag_save_in"))
-
 {
 }
 

framework/src/outputs/AdvancedOutput.C

@@ -484,7 +484,7 @@ AdvancedOutput::initAvailableLists()
 void
 AdvancedOutput::initExecutionTypes(const std::string & name, ExecFlagEnum & input)
 {
-  // Build the input paramemter name
+  // Build the input parameter name
   std::string param_name = "execute_";
   param_name += name + "_on";
 

framework/src/userobjects/PropertyReadFile.C

@@ -303,7 +303,7 @@ PropertyReadFile::getElementData(const Elem * elem, unsigned int prop_num) const
                jelem,
                " greater than than total number of element in mesh: ",
                _mesh.nElem(),
-               ". Elements should be numbered consecutively.");
+               ". Elements should be numbered consecutively, and ids should start from 0.");
   return _reader.getData(jelem)[prop_num];
 }
 
@@ -316,7 +316,7 @@ PropertyReadFile::getNodeData(const Node * const node, const unsigned int prop_n
                jnode,
                " greater than than total number of nodes in mesh: ",
                _mesh.nNodes(),
-               ". Nodes should be numbered consecutively.");
+               ". Nodes should be numbered consecutively, with ids starting from 0.");
   return _reader.getData(jnode)[prop_num];
 }
 
@@ -334,7 +334,7 @@ PropertyReadFile::getBlockData(const Elem * elem, unsigned int prop_num) const
                elem_subdomain_id,
                " greater than than total number of blocks in mesh: ",
                _nblock,
-               ". Blocks should be numbered consecutively.");
+               ". Blocks should be numbered consecutively, starting from 0.");
   return _reader.getData(elem_subdomain_id - offset)[prop_num];
 }
 

modules/navier_stokes/doc/content/source/fvbcs/NSFVHeatFluxBC.md

@@ -70,7 +70,7 @@ must be provided:
 - `splitting = effective_thermal_conductivity`, `locality = global`: `average_kappa` and `average_kappa_solid` postprocessors,
    as well as `average_eps` postprocessor due to internal representation of $\kappa_f$ as $\epsilon\tilde{\kappa}_f$.
 
-! alert note
+!alert note
 To protect against cases where at the first time step
 the thermal conductivity or effective thermal conductivity might not have yet
 been initialized, or cases where the coupled postprocessors have not yet been

modules/solid_mechanics/doc/content/source/materials/ADComputeFiniteStrain.md

@@ -42,7 +42,7 @@ where $\boldsymbol{x}_{n+1}$ is the position vector of materials points in
 $\kappa_{n+1}$, and $\boldsymbol{x}_{n}$ is the position vector of materials
 points in $\kappa_{n}$.
 
-alert note title=Incremental vs Total Deformation Gradient
+!alert note title=Incremental vs Total Deformation Gradient
 Note that $\hat{\boldsymbol{F}}$ is NOT the deformation gradient, but rather the
 incremental deformation gradient of $\kappa_{n+1}$ with respect to $\kappa_n$.
 Thus $\hat{\boldsymbol{F}} = \boldsymbol{F}_{n+1} \boldsymbol{F}_n^{-1}$, where

test/tests/fviks/auxiliary_variables/fv_reaction_1D.i

@@ -0,0 +1,88 @@
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 2
+  []
+  [subdomain1]
+    input = gen
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '1.0 0 0'
+    block_id = 1
+    top_right = '2.0 1.0 0'
+  []
+  [interface]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = 'subdomain1'
+    primary_block = '0'
+    paired_block = '1'
+    new_boundary = 'primary0_interface'
+  []
+[]
+
+[Variables]
+  [u]
+    type = MooseVariableFVReal
+    block = '0'
+  []
+[]
+
+[AuxVariables]
+  [v]
+    type = MooseVariableFVReal
+    block = '1'
+    initial_condition = 4
+  []
+[]
+
+[FVKernels]
+  [diff_u]
+    type = FVDiffusion
+    variable = u
+    block = '0'
+    coeff = 1
+  []
+[]
+
+[FVInterfaceKernels]
+  [interface]
+    type = FVDiffusionInterface
+    variable1 = u
+    variable2 = 'v'
+    boundary = 'primary0_interface'
+    coeff1 = 1
+    coeff2= 2
+    subdomain1 = 0
+    subdomain2 = 1
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = PJFNK
+  nl_rel_tol = 1e-10
+  nl_forced_its = 2
+[]
+
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Outputs]
+  csv = true
+[]
+
+[Postprocessors]
+  [min]
+    type = ElementExtremeValue
+    variable = 'u'
+    value_type = 'min'
+    block = '0'
+  []
+  [max]
+    type = ElementExtremeValue
+    variable = 'u'
+    block = '0'
+  []
+[]

test/tests/fviks/auxiliary_variables/gold/fv_reaction_1D_out.csv

@@ -0,0 +1,3 @@
+time,max,min
+0,0,0
+1,4,4

test/tests/fviks/auxiliary_variables/tests

@@ -0,0 +1,20 @@
+[Tests]
+  issues = '#28977'
+  design = 'FVInterfaceKernels/index.md'
+  [aux]
+    requirement = 'The system shall be able to solve problems with auxiliary variables coupled through interface kernels'
+    [fv]
+      type = CSVDiff
+      input = 'fv_reaction_1D.i'
+      csvdiff= 'fv_reaction_1D_out.csv'
+      detail = 'to finite volume variables,'
+    []
+    [fv_res_jac]
+      type = CSVDiff
+      input = 'fv_reaction_1D.i'
+      csvdiff= 'fv_reaction_1D_out.csv'
+      cli_args = 'Executioner/residual_and_jacobian_together=true'
+      detail = 'to finite volume variables when computing the residual and Jacobian together,'
+    []
+  []
+[]

test/tests/interfacekernels/auxiliary_variables/gold/reaction_1D_steady_out.csv

@@ -0,0 +1,3 @@
+time,max,min
+0,0,0
+1,4,4