Add Alex's idea on setting the systems. (#25599)

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

@@ -188,8 +188,11 @@ FVInterfaceKernel::computeResidual(const FaceInfo & fi)
 
   const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
 
-  addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true);
-  if (_var1.sys().number() == _var2.sys().number())
+  // If the two variables belong to two different nonlinear systems, we only contribute to the one
+  // which is being assembled right now
+  if (_var1.sys().number() == _subproblem.currentNlSysNum())
+    addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true);
+  if (_var2.sys().number() == _subproblem.currentNlSysNum())
     addResidual(_elem_is_one ? -r : r, _var2.number(), _elem_is_one ? true : false);
 }
 
@@ -206,11 +209,14 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi)
 
   const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
 
-  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())
+  // If the two variables belong to two different nonlinear systems, we only contribute to the one
+  // which is being assembled right now
+  if (_var1.sys().number() == _subproblem.currentNlSysNum())
+    addResidualsAndJacobian(_assembly,
+                            std::array<ADReal, 1>{{_elem_is_one ? r : -r}},
+                            _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(),
+                            _var1.scalingFactor());
+  if (_var2.sys().number() == _subproblem.currentNlSysNum())
     addResidualsAndJacobian(_assembly,
                             std::array<ADReal, 1>{{_elem_is_one ? -r : r}},
                             _elem_is_one ? _var2.dofIndicesNeighbor() : _var2.dofIndices(),

framework/src/problems/FEProblemBase.C

@@ -2929,18 +2929,10 @@ FEProblemBase::addFVInterfaceKernel(const std::string & fv_ik_name,
                                     const std::string & name,
                                     InputParameters & parameters)
 {
-  /// Checking the nonlinear system number of variable1
-  /// If variable1 and variable2 live on different systems, the interface kernel
-  /// is added to the system of variable1 alone. Another interface kernel needs to be
-  /// created with flipped variables for the system of variable2
-  const auto nl_sys_num =
-      determineNonlinearSystem(parameters.varName("variable1", name), true).first
-          ? determineNonlinearSystem(parameters.varName("variable1", name), true).second
-          : (unsigned int)0;
-
-  /// This is set because the Attribute system filters IKs based on this
-  parameters.set<SystemBase *>("_sys") = _nl[nl_sys_num].get();
-  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 ////
@@ -3567,21 +3559,21 @@ 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 bool system_already_set = parameters.get<SystemBase *>("_sys");
   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") &&
       parameters.get<bool>("use_displaced_mesh"))
   {
     parameters.set<SubProblem *>("_subproblem") = _displaced_problem.get();
-    if (!system_already_set)
-      parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys(nl_sys_num);
+    parameters.set<SystemBase *>("_sys") = &_displaced_problem->nlSys(nl_sys_num);
   }
   else
   {
@@ -3593,8 +3585,7 @@ FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::st
       parameters.set<bool>("use_displaced_mesh") = false;
 
     parameters.set<SubProblem *>("_subproblem") = this;
-    if (!system_already_set)
-      parameters.set<SystemBase *>("_sys") = _nl[nl_sys_num].get();
+    parameters.set<SystemBase *>("_sys") = _nl[nl_sys_num].get();
   }
 }
 

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,7 +37,8 @@ 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(
@@ -51,9 +52,8 @@ FVConvectionCorrelationInterface::computeQpResidual()
   // 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
-  bool var1_is_fluid = ("wraps_" + var1().name() == _temp_fluid.functorName());
   const Elem * elem_on_fluid_side =
-      (elemIsOne() && var1_is_fluid) || (!elemIsOne() && !var1_is_fluid)
+      (elemIsOne() && _var1_is_fluid) || (!elemIsOne() && !_var1_is_fluid)
           ? &_face_info->elem()
           : _face_info->neighborPtr();
 
@@ -77,11 +77,11 @@ FVConvectionCorrelationInterface::computeQpResidual()
   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_is_fluid ? var1() : var2()).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 fluid_side = singleSidedFaceArg(var1_is_fluid ? var1() : var2(), _face_info);
-  const auto solid_side = singleSidedFaceArg(var1_is_fluid ? var2() : var1(), _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);