Merge pull request #25605 from grmnptr/multisys-fv-interface-25599

Enable FV interface kernels for multiple nonlinear systems.

framework/doc/content/source/fviks/FVDiffusionInterface.md

@@ -5,6 +5,10 @@
 The diffusive flux is obtained from a two point gradient, and the diffusivity is
 interpolated to the interface.
 
+!alert note
+This kernel supports interfaces between variables which belong to different nonlinear systems.
+For instructions on how to set these cases up, visit the [FVInterfaceKernels syntax page](syntax/FVInterfaceKernels/index.md).
+
 ## Example input file syntax
 
 In this example, two diffusion problems with a source terms are solved on each side

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

@@ -28,6 +28,7 @@ side. So making use of the `subdomain` parameters, we provide a protected method
 called `elemIsOne()` that returns a boolean indicating whether the
 `FaceInfo::elem` side of the interface corresponds to the `subdomain1` side of
 the interface. This allows the developer to write code like the following:
+
 ```
 FVFooInterface::FVFooInterface(const InputParameters & params)
   : FVInterfaceKernel(params),
@@ -47,5 +48,16 @@ FVFooInterface::computeQpResidual()
   /// Code that uses coef_elem and coef_neighbor
 }
 ```
+
 and have confidence that they have good data in `coef_elem` and `coef_neighbor`
 and have clarity about what is happening in their code.
+
+!alert! note
+When using an FVInterfaceKernel which connects variables that belong to different nonlinear systems,
+create two kernels with flipped variable and material property parameters. The reason behind this
+is that the interface kernel will only contribute to the system which `variable1` belongs to.
+For an example, see:
+
+!listing /test/tests/fviks/diffusion/multisystem.i
+
+!alert-end!

framework/include/fviks/FVInterfaceKernel.h

@@ -99,12 +99,6 @@ class FVInterfaceKernel : public MooseObject,
    */
   const std::set<SubdomainID> & sub2() const { return _subdomain2; }
 
-  /**
-   * @return The system associated with this object. Either an undisplaced or displaced nonlinear
-   * system
-   */
-  const SystemBase & sys() const { return _sys; }
-
   /**
    * @return Whether the \p FaceInfo element is on the 1st side of the interface
    */
@@ -192,16 +186,13 @@ class FVInterfaceKernel : public MooseObject,
   /// The SubProblem
   SubProblem & _subproblem;
 
-  /// the system object
-  SystemBase & _sys;
+  MooseVariableFV<Real> & _var1;
+  MooseVariableFV<Real> & _var2;
 
   /// The Assembly object
   Assembly & _assembly;
 
 private:
-  MooseVariableFV<Real> & _var1;
-  MooseVariableFV<Real> & _var2;
-
   std::set<SubdomainID> _subdomain1;
   std::set<SubdomainID> _subdomain2;
 

framework/include/problems/FEProblemBase.h

@@ -851,13 +851,15 @@ class FEProblemBase : public SubProblem, public Restartable
    * @param name Name for the object to be created
    * @param parameters InputParameters for the object
    * @param threaded Whether or not to create n_threads copies of the object
+   * @param var_param_name The name of the parameter on the object which holds the primary variable.
    * @return A vector of shared_ptrs to the added objects
    */
   template <typename T>
   std::vector<std::shared_ptr<T>> addObject(const std::string & type,
                                             const std::string & name,
                                             InputParameters & parameters,
-                                            const bool threaded = true);
+                                            const bool threaded = true,
+                                            const std::string & var_param_name = "variable");
 
   // Postprocessors /////
   virtual void addPostprocessor(const std::string & pp_name,
@@ -2295,7 +2297,9 @@ class FEProblemBase : public SubProblem, public Restartable
    *
    * This is needed due to header includes/forward declaration issues
    */
-  void addObjectParamsHelper(InputParameters & params, const std::string & object_name);
+  void addObjectParamsHelper(InputParameters & params,
+                             const std::string & object_name,
+                             const std::string & var_param_name = "variable");
 
 #ifdef LIBMESH_ENABLE_AMR
   Adaptivity _adaptivity;
@@ -2585,12 +2589,13 @@ std::vector<std::shared_ptr<T>>
 FEProblemBase::addObject(const std::string & type,
                          const std::string & name,
                          InputParameters & parameters,
-                         const bool threaded)
+                         const bool threaded,
+                         const std::string & var_param_name)
 {
   parallel_object_only();
 
   // Add the _subproblem and _sys parameters depending on use_displaced_mesh
-  addObjectParamsHelper(parameters, name);
+  addObjectParamsHelper(parameters, name, var_param_name);
 
   const auto n_threads = threaded ? libMesh::n_threads() : 1;
   std::vector<std::shared_ptr<T>> objects(n_threads);

framework/src/fviks/FVInterfaceKernel.C

@@ -98,13 +98,19 @@ FVInterfaceKernel::FVInterfaceKernel(const InputParameters & parameters)
     ADFunctorInterface(this),
     _tid(getParam<THREAD_ID>("_tid")),
     _subproblem(*getCheckedPointerParam<SubProblem *>("_subproblem")),
-    _sys(*getCheckedPointerParam<SystemBase *>("_sys")),
-    _assembly(_subproblem.assembly(_tid, _sys.number())),
-    _var1(_sys.getFVVariable<Real>(_tid, getParam<NonlinearVariableName>("variable1"))),
-    _var2(_sys.getFVVariable<Real>(_tid,
+    _var1(_subproblem.getVariable(_tid, getParam<NonlinearVariableName>("variable1"))
+              .sys()
+              .getFVVariable<Real>(_tid, getParam<NonlinearVariableName>("variable1"))),
+    _var2(_subproblem
+              .getVariable(_tid,
+                           isParamValid("variable2") ? getParam<NonlinearVariableName>("variable2")
+                                                     : getParam<NonlinearVariableName>("variable1"))
+              .sys()
+              .getFVVariable<Real>(_tid,
                                    isParamValid("variable2")
                                        ? getParam<NonlinearVariableName>("variable2")
                                        : getParam<NonlinearVariableName>("variable1"))),
+    _assembly(_subproblem.assembly(_tid, _var1.sys().number())),
     _mesh(_subproblem.mesh())
 {
   if (getParam<bool>("use_displaced_mesh"))
@@ -180,13 +186,15 @@ FVInterfaceKernel::computeResidual(const FaceInfo & fi)
 {
   setupData(fi);
 
-  const auto var_elem_num = _elem_is_one ? _var1.number() : _var2.number();
-  const auto var_neigh_num = _elem_is_one ? _var2.number() : _var1.number();
-
   const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
 
-  addResidual(r, var_elem_num, false);
-  addResidual(-r, var_neigh_num, true);
+  // If the two variables belong to two different nonlinear systems, we only contribute to the one
+  // which is being assembled right now
+  mooseAssert(_var1.sys().number() == _subproblem.currentNlSysNum(),
+              "The interface kernel should contribute to the system which variable1 belongs to!");
+  addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true);
+  if (_var1.sys().number() == _var2.sys().number())
+    addResidual(_elem_is_one ? -r : r, _var2.number(), _elem_is_one ? true : false);
 }
 
 void
@@ -200,20 +208,21 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi)
 {
   setupData(fi);
 
-  const auto & elem_dof_indices = _elem_is_one ? _var1.dofIndices() : _var2.dofIndices();
-  const auto & neigh_dof_indices =
-      _elem_is_one ? _var2.dofIndicesNeighbor() : _var1.dofIndicesNeighbor();
-  mooseAssert((elem_dof_indices.size() == 1) && (neigh_dof_indices.size() == 1),
-              "We're currently built to use CONSTANT MONOMIALS");
-  const auto elem_scaling_factor = _elem_is_one ? _var1.scalingFactor() : _var2.scalingFactor();
-  const auto neigh_scaling_factor = _elem_is_one ? _var2.scalingFactor() : _var1.scalingFactor();
-
   const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
 
-  addResidualsAndJacobian(
-      _assembly, std::array<ADReal, 1>{{r}}, elem_dof_indices, elem_scaling_factor);
-  addResidualsAndJacobian(
-      _assembly, std::array<ADReal, 1>{{-r}}, neigh_dof_indices, neigh_scaling_factor);
+  // If the two variables belong to two different nonlinear systems, we only contribute to the one
+  // which is being assembled right now
+  mooseAssert(_var1.sys().number() == _subproblem.currentNlSysNum(),
+              "The interface kernel should contribute to the system which variable1 belongs to!");
+  addResidualsAndJacobian(_assembly,
+                          std::array<ADReal, 1>{{_elem_is_one ? r : -r}},
+                          _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(),
+                          _var1.scalingFactor());
+  if (_var1.sys().number() == _var2.sys().number())
+    addResidualsAndJacobian(_assembly,
+                            std::array<ADReal, 1>{{_elem_is_one ? -r : r}},
+                            _elem_is_one ? _var2.dofIndicesNeighbor() : _var2.dofIndices(),
+                            _var2.scalingFactor());
 }
 
 Moose::ElemArg

framework/src/fviks/FVScalarLagrangeMultiplierInterface.C

@@ -28,6 +28,8 @@ FVScalarLagrangeMultiplierInterface::FVScalarLagrangeMultiplierInterface(
     _lambda_var(*getScalarVar("lambda", 0)),
     _lambda(adCoupledScalarValue("lambda"))
 {
+  if (var1().sys().number() != var2().sys().number())
+    mooseError(this->type(), " does not support multiple nonlinear systems!");
 }
 
 void

framework/src/problems/FEProblemBase.C

@@ -2929,7 +2929,10 @@ FEProblemBase::addFVInterfaceKernel(const std::string & fv_ik_name,
                                     const std::string & name,
                                     InputParameters & parameters)
 {
-  addObject<FVInterfaceKernel>(fv_ik_name, name, parameters);
+  /// We assume that variable1 and variable2 can live on different systems, in this case
+  /// the user needs to create two interface kernels with flipped variables and parameters
+  addObject<FVInterfaceKernel>(
+      fv_ik_name, name, parameters, /*threaded=*/true, /*variable_param_name=*/"variable1");
 }
 
 // InterfaceKernels ////
@@ -3556,12 +3559,14 @@ FEProblemBase::swapBackMaterialsNeighbor(THREAD_ID tid)
 }
 
 void
-FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::string & object_name)
+FEProblemBase::addObjectParamsHelper(InputParameters & parameters,
+                                     const std::string & object_name,
+                                     const std::string & var_param_name)
 {
   const auto nl_sys_num =
-      parameters.isParamValid("variable") &&
-              determineNonlinearSystem(parameters.varName("variable", object_name)).first
-          ? determineNonlinearSystem(parameters.varName("variable", object_name)).second
+      parameters.isParamValid(var_param_name) &&
+              determineNonlinearSystem(parameters.varName(var_param_name, object_name)).first
+          ? determineNonlinearSystem(parameters.varName(var_param_name, object_name)).second
           : (unsigned int)0;
 
   if (_displaced_problem && parameters.have_parameter<bool>("use_displaced_mesh") &&

modules/navier_stokes/doc/content/source/fviks/FVConvectionCorrelationInterface.md

@@ -10,6 +10,10 @@ with $q_s$ the surface convective heat flux, $h_{correlation}$ the heat transfer
 defined by the correlation as a material property and $T_{solid/fluid}$ the temperature of the
 adjacent solid and fluid.
 
+!alert note
+This kernel supports interfaces between variables which belong to different nonlinear systems.
+For instructions on how to set these cases up, visit the [FVInterfaceKernels syntax page](syntax/FVInterfaceKernels/index.md).
+
 ## Example input file syntax
 
 In this example, a cold fluid is flowing next to a centrally-heated solid region. The heat diffuses

modules/navier_stokes/include/fviks/FVConvectionCorrelationInterface.h

@@ -35,4 +35,7 @@ class FVConvectionCorrelationInterface : public FVInterfaceKernel
 
   /// libmesh object to find points in the mesh
   std::unique_ptr<PointLocatorBase> _pl;
+
+  /// Boolean to see if variable1 is the fluid
+  const bool _var1_is_fluid;
 };

modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C

@@ -37,52 +37,58 @@ FVConvectionCorrelationInterface::FVConvectionCorrelationInterface(const InputPa
     _htc(getFunctor<ADReal>("h")),
     _bulk_distance(getParam<Real>("bulk_distance")),
     _use_wall_cell(getParam<bool>("wall_cell_is_bulk")),
-    _pl(mesh().getPointLocator())
+    _pl(mesh().getPointLocator()),
+    _var1_is_fluid("wraps_" + var1().name() == _temp_fluid.functorName())
 {
   if (!_use_wall_cell && (_bulk_distance < 0))
     mooseError(
         "The bulk distance should be specified or 'wall_cell_is_bulk' should be set to true for "
         "the FVTwoVarConvectionCorrelationInterface");
-  if ("wraps_" + var1().name() != _temp_fluid.functorName())
-    paramError(NS::T_fluid, "variable1 must be equal to T_fluid parameter.");
-  if ("wraps_" + var2().name() != _temp_solid.functorName())
-    paramError(NS::T_solid, "variable2 must be equal to T_solid parameter.");
 }
 
 ADReal
 FVConvectionCorrelationInterface::computeQpResidual()
 {
-  const Elem * current_elem = elemIsOne() ? &_face_info->elem() : _face_info->neighborPtr();
+  // If variable1 is fluid and variable 1 is on elem or
+  // if variable2 is fluid and variable 2 is on elem
+  // the fluid element will be elem otherwise it is the neighbor
+  const Elem * elem_on_fluid_side =
+      (elemIsOne() && _var1_is_fluid) || (!elemIsOne() && !_var1_is_fluid)
+          ? &_face_info->elem()
+          : _face_info->neighborPtr();
+
   const Elem * bulk_elem;
   const auto state = determineState();
   if (!_use_wall_cell)
   {
     Point p = _face_info->faceCentroid();
     Point du = Point(MetaPhysicL::raw_value(_normal));
     du *= _bulk_distance;
-    if (elemIsOne())
+    // The normal always points outwards from the elem (towards the neighbor)
+    if (elem_on_fluid_side == &_face_info->elem())
       p -= du;
     else
       p += du;
     bulk_elem = (*_pl)(p);
   }
   else
-    bulk_elem = current_elem;
+    bulk_elem = elem_on_fluid_side;
 
   mooseAssert(bulk_elem,
               "The element at bulk_distance from the wall was not found in the mesh. "
               "Increase the number of ghost layers with the 'ghost_layers' parameter.");
-  mooseAssert(var1().hasBlocks(bulk_elem->subdomain_id()),
+  mooseAssert((_var1_is_fluid ? var1() : var2()).hasBlocks(bulk_elem->subdomain_id()),
               "The fluid temperature is not defined at bulk_distance from the wall.");
 
-  const auto face_arg_side1 = singleSidedFaceArg(var1(), _face_info);
-  const auto face_arg_side2 = singleSidedFaceArg(var2(), _face_info);
+  const auto fluid_side = singleSidedFaceArg(_var1_is_fluid ? var1() : var2(), _face_info);
+  const auto solid_side = singleSidedFaceArg(_var1_is_fluid ? var2() : var1(), _face_info);
+
   const auto bulk_elem_arg = makeElemArg(bulk_elem);
 
   /// We make sure that the gradient*normal part is addressed
   auto multipler =
       _normal * (_face_info->faceCentroid() - bulk_elem->vertex_average()) > 0 ? 1 : -1;
 
-  return multipler * _htc(face_arg_side1, state) *
-         (_temp_fluid(bulk_elem_arg, state) - _temp_solid(face_arg_side2, state));
+  return multipler * _htc(fluid_side, state) *
+         (_temp_fluid(bulk_elem_arg, state) - _temp_solid(solid_side, state));
 }

test/include/executioners/SteadySolve2.h

@@ -42,4 +42,5 @@ class SteadySolve2 : public Executioner
   const unsigned int _first_nl_sys;
   const unsigned int _second_nl_sys;
   bool _last_solve_converged;
+  unsigned int _number_of_iterations;
 };

test/src/executioners/SteadySolve2.C

@@ -27,6 +27,10 @@ SteadySolve2::validParams()
                                                "The first nonlinear system to solve");
   params.addRequiredParam<NonlinearSystemName>("second_nl_sys_to_solve",
                                                "The second nonlinear system to solve");
+  params.addRangeCheckedParam<unsigned int>("number_of_iterations",
+                                            1,
+                                            "number_of_iterations>0",
+                                            "The number of iterations between the two systems.");
   return params;
 }
 
@@ -37,7 +41,8 @@ SteadySolve2::SteadySolve2(const InputParameters & parameters)
     _time_step(_problem.timeStep()),
     _time(_problem.time()),
     _first_nl_sys(_problem.nlSysNum(getParam<NonlinearSystemName>("first_nl_sys_to_solve"))),
-    _second_nl_sys(_problem.nlSysNum(getParam<NonlinearSystemName>("second_nl_sys_to_solve")))
+    _second_nl_sys(_problem.nlSysNum(getParam<NonlinearSystemName>("second_nl_sys_to_solve"))),
+    _number_of_iterations(getParam<unsigned int>("number_of_iterations"))
 {
   _fixed_point_solve->setInnerSolve(_feproblem_solve);
 
@@ -107,12 +112,16 @@ SteadySolve2::execute()
       return _problem.nlConverged(sys_num);
     };
 
-    const bool converged = solve(_first_nl_sys) && solve(_second_nl_sys);
+    bool converged = true;
 
-    if (!converged)
+    for (unsigned int i = 0; i < _number_of_iterations; i++)
     {
-      _console << "Aborting as solve did not converge" << std::endl;
-      break;
+      converged = converged && solve(_first_nl_sys) && solve(_second_nl_sys);
+      if (!converged)
+      {
+        _console << "Aborting as solve did not converge" << std::endl;
+        break;
+      }
     }
 
     _problem.computeIndicators();

test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C

@@ -27,6 +27,8 @@ FVOnlyAddDiffusionToOneSideOfInterface::FVOnlyAddDiffusionToOneSideOfInterface(
     const InputParameters & params)
   : FVInterfaceKernel(params), _coeff2(getFunctor<ADReal>("coeff2"))
 {
+  if (var1().sys().number() != var2().sys().number())
+    mooseError(this->type(), " does not support multiple nonlinear systems!");
 }
 
 void