diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C
index 33ab95644c68..e51aa4a9bbdd 100644
--- a/framework/src/systems/SystemBase.C
+++ b/framework/src/systems/SystemBase.C
@@ -1356,7 +1356,9 @@ SystemBase::solutionState(const unsigned int state,
                " was requested in ",
                name(),
                " but only up to state ",
-               _solution_states[static_cast<unsigned short>(iteration_type)].size() - 1,
+               (_solution_states[static_cast<unsigned short>(iteration_type)].size() == 0)
+                   ? 0
+                   : _solution_states[static_cast<unsigned short>(iteration_type)].size() - 1,
                " is available.");
 
   const auto & solution_states = _solution_states[static_cast<unsigned short>(iteration_type)];
@@ -1385,6 +1387,9 @@ SystemBase::needSolutionState(const unsigned int state,
                               const Moose::SolutionIterationType iteration_type)
 {
   libmesh_parallel_only(this->comm());
+  mooseAssert(!Threads::in_threads,
+              "This routine is not thread-safe. Request the solution state before using it in "
+              "a threaded region.");
 
   if (hasSolutionState(state, iteration_type))
     return;
diff --git a/framework/src/variables/MooseVariableFV.C b/framework/src/variables/MooseVariableFV.C
index 214683214519..8e36c026dc79 100644
--- a/framework/src/variables/MooseVariableFV.C
+++ b/framework/src/variables/MooseVariableFV.C
@@ -461,9 +461,10 @@ MooseVariableFV<OutputType>::getElemValue(const Elem * const elem, const StateAr
   // which is wrong during things like finite difference Jacobian evaluation, e.g. when PETSc
   // perturbs the solution vector we feed these perturbations into the current_local_solution
   // while the libMesh solution is frozen in the non-perturbed state
-  const auto & global_soln = (state.state == 0)
-                                 ? *this->_sys.currentSolution()
-                                 : this->_sys.solutionState(state.state, state.iteration_type);
+  const auto & global_soln =
+      (state.state == 0)
+          ? *this->_sys.currentSolution()
+          : std::as_const(this->_sys).solutionState(state.state, state.iteration_type);
 
   ADReal value = global_soln(index);
 
diff --git a/modules/navier_stokes/src/base/NavierStokesMethods.C b/modules/navier_stokes/src/base/NavierStokesMethods.C
index 4af46984c142..74d525576f04 100644
--- a/modules/navier_stokes/src/base/NavierStokesMethods.C
+++ b/modules/navier_stokes/src/base/NavierStokesMethods.C
@@ -156,7 +156,7 @@ computeSpeed(const ADRealVectorValue & velocity)
 
 /// Bounded element maps for wall treatment
 void
-getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_name,
+getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_names,
                        const FEProblemBase & fe_problem,
                        const SubProblem & subproblem,
                        const std::set<SubdomainID> & block_ids,
@@ -164,6 +164,7 @@ getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_name,
 {
 
   wall_bounded_map.clear();
+  const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_names);
 
   for (const auto & elem : fe_problem.mesh().getMesh().active_element_ptr_range())
   {
@@ -171,9 +172,8 @@ getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_name,
       for (const auto i_side : elem->side_index_range())
       {
         const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
-        for (const auto & name : wall_boundary_name)
+        for (const auto & wall_id : wall_boundary_ids)
         {
-          const auto wall_id = subproblem.mesh().getBoundaryID(name);
           for (const auto side_id : side_bnds)
             if (side_id == wall_id)
               wall_bounded_map[elem] = true;