Repository: OpenMDAO/OpenMDAO

Issue:
AutoIVC and src_indices errors upon run_model
### Description

I've ran into two separate (but maybe related) bugs lately, which I was able to reproduce with the provided example script. Both occur during run_model, not setup.

The first is that, when some group's component input is unconnected, I get the error
```
Traceback (most recent call last):
  File "/lustre3/hpnobackup2/asthelen/Codes/Testing/OM_bugs/ivc/main.py", line 22, in <module>
    prob.run_model()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/problem.py", line 680, in run_model
    self.final_setup()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/problem.py", line 1088, in final_setup
    model._setup_part2()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/group.py", line 1103, in _setup_part2
    self._resolve_group_input_defaults()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/group.py", line 1909, in _resolve_group_input_defaults
    top_origin = paths[0][1]
IndexError: list index out of range
```
It looks like in the n2, the two components' ```in1``` inputs get turned into AutoIVCs as they're supposed to... the model just can't be ran.

If you uncomment the two connect lines in the example, then I get this other error:
```
Traceback (most recent call last):
  File "/lustre3/hpnobackup2/asthelen/Codes/Testing/OM_bugs/ivc/main.py", line 22, in <module>
    prob.run_model()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/problem.py", line 680, in run_model
    self.final_setup()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/problem.py", line 1089, in final_setup
    self._check_collected_errors()
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/problem.py", line 638, in _check_collected_errors
    raise exc_type('\n'.join(final_msg)).with_traceback(tback)
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/core/group.py", line 1641, in _resolve_src_indices
    sinds = convert_src_inds(current_pinfo.src_indices, current_pinfo.src_shape,
  File "/hpnobackup2/asthelen/Codes/Repositories/OM_test/OpenMDAO_2025jan24/openmdao/utils/general_utils.py", line 1184, in convert_src_inds
    return parent_src_inds.shaped_array(flat=True)[my_src_inds.flat()]
IndexError: 
Collected errors for problem 'main':
   <model> <class Top>: When connecting 'ivc.in1' to 'comp2.in1': input 'comp2.in1' src_indices are [1] and indexing into those failed using src_indices [1] from input 'comp2.in1'. Error was: index 1 is out of bounds for axis 0 with size 1
```
It's only an issue with src_indices=[1] and not [0], so if you set them both to [0] then the error doesn't occur.

### Example

```
import numpy as np
import openmdao.api as om

class Comp(om.ExplicitComponent):
    def setup(self):
        self.add_input("in1", val=1.)
        self.add_output("out1", val=1.)

class Top(om.Group):
    def setup(self):
        self.add_subsystem("ivc", om.IndepVarComp(), promotes=["*"])
        self.ivc.add_output("in1", [1.,1.])
        self.add_subsystem("comp1", Comp())
        self.add_subsystem("comp2", Comp())

        #self.connect("in1", "comp1.in1", src_indices=[0])
        #self.connect("in1", "comp2.in1", src_indices=[1])

prob = om.Problem()
prob.model = Top()
prob.setup(mode='rev')
om.n2(prob, show_browser=False, outfile='n2.html')
prob.run_model()
```

### OpenMDAO Version

3.38.1-dev

### Relevant environment information

_No response_

Hints from issue comments:

Solution patch:
diff --git a/openmdao/components/explicit_func_comp.py b/openmdao/components/explicit_func_comp.py
index d38d1bad4..4298012c1 100644
--- a/openmdao/components/explicit_func_comp.py
+++ b/openmdao/components/explicit_func_comp.py
@@ -201,7 +201,7 @@ class ExplicitFuncComp(ExplicitComponent):
             else:
                 j = [np.asarray(a).reshape((a.shape[0], shape_to_len(a.shape[1:])))
                      for a in jac_reverse(func, argnums, tangents)(*invals)]
-                j = coloring.expand_jac(np.hstack(j), 'rev')
+                j = coloring._expand_jac(np.hstack(j), 'rev').toarray()
         else:
             tangents = self._get_tangents(invals, 'fwd', coloring, argnums)
             if coloring is None:
@@ -222,7 +222,7 @@ class ExplicitFuncComp(ExplicitComponent):
             else:
                 j = [np.asarray(a).reshape((shape_to_len(a.shape[:-1]), a.shape[-1]))
                      for a in jac_forward(func, argnums, tangents)(*invals)]
-                j = coloring.expand_jac(np.vstack(j), 'fwd')
+                j = coloring._expand_jac(np.vstack(j), 'fwd').toarray()
 
         self._jacobian.set_dense_jac(self, j)
 
diff --git a/openmdao/components/implicit_func_comp.py b/openmdao/components/implicit_func_comp.py
index 37fb11509..0fc66a585 100644
--- a/openmdao/components/implicit_func_comp.py
+++ b/openmdao/components/implicit_func_comp.py
@@ -251,7 +251,7 @@ class ImplicitFuncComp(ImplicitComponent):
                 j = [np.asarray(a).reshape((a.shape[0], shape_to_len(a.shape[1:])))
                      for a in jac_reverse(self._apply_nonlinear_func_jax, argnums,
                                           tangents)(*invals)]
-                j = coloring.expand_jac(np.hstack(self._reorder_col_chunks(j)), 'rev')
+                j = coloring._expand_jac(np.hstack(self._reorder_col_chunks(j)), 'rev').toarray()
             else:
                 j = []
                 for a in jac_reverse(self._apply_nonlinear_func_jax, argnums, tangents)(*invals):
@@ -270,7 +270,7 @@ class ImplicitFuncComp(ImplicitComponent):
                 j = [np.asarray(a).reshape((shape_to_len(a.shape[:-1]), a.shape[-1]))
                      for a in jac_forward(self._apply_nonlinear_func_jax, argnums,
                                           tangents)(*invals)]
-                j = coloring.expand_jac(np.vstack(j), 'fwd')
+                j = coloring._expand_jac(np.vstack(j), 'fwd').toarray()
             else:
                 tangents = self._get_tangents(invals, 'fwd', coloring, argnums)
                 j = []
diff --git a/openmdao/components/jax_explicit_comp.py b/openmdao/components/jax_explicit_comp.py
index a2a37915b..0fd750a7a 100644
--- a/openmdao/components/jax_explicit_comp.py
+++ b/openmdao/components/jax_explicit_comp.py
@@ -11,7 +11,7 @@ from openmdao.core.explicitcomponent import ExplicitComponent
 from openmdao.utils.om_warnings import issue_warning
 from openmdao.utils.jax_utils import jax, jit, \
     _jax_register_pytree_class, _compute_sparsity, get_vmap_tangents, \
-    _update_subjac_sparsity, _jax_derivs2partials, _uncompress_jac, _jax2np, \
+    _update_subjac_sparsity, _jax_derivs2partials, _jax2np, \
     _ensure_returns_tuple, _compute_output_shapes, _update_add_input_kwargs, \
     _update_add_output_kwargs, _get_differentiable_compute_primal, _re_init
 from openmdao.utils.code_utils import get_return_names, get_function_deps
@@ -385,7 +385,12 @@ class JaxExplicitComponent(ExplicitComponent):
             The partials to compute.
         """
         J = self._jac_func_(self._tangents['fwd'], tuple(inputs.values()))
-        partials.set_dense_jac(self, _uncompress_jac(self, _jax2np(J), 'fwd'))
+        J = _jax2np(J)
+        if self._coloring_info.coloring is not None:
+            J = self._coloring_info.coloring._expand_jac(J, 'fwd')
+            partials.set_csc_jac(self, J)
+        else:
+            partials.set_dense_jac(self, J)
 
     def _jacrev_colored(self, inputs, partials):
         """
@@ -399,7 +404,12 @@ class JaxExplicitComponent(ExplicitComponent):
             The partials to compute.
         """
         J = self._jac_func_(self._tangents['rev'], tuple(inputs.values()))
-        partials.set_dense_jac(self, _uncompress_jac(self, _jax2np(J).T, 'rev'))
+        J = _jax2np(J).T
+        if self._coloring_info.coloring is not None:
+            J = self._coloring_info.coloring._expand_jac(J, 'rev')
+            partials.set_csc_jac(self, J)
+        else:
+            partials.set_dense_jac(self, J)
 
     def compute_sparsity(self, direction=None, num_iters=1, perturb_size=1e-9):
         """
diff --git a/openmdao/components/jax_implicit_comp.py b/openmdao/components/jax_implicit_comp.py
index 904373773..2b344386f 100644
--- a/openmdao/components/jax_implicit_comp.py
+++ b/openmdao/components/jax_implicit_comp.py
@@ -13,7 +13,7 @@ from openmdao.utils.om_warnings import issue_warning
 from openmdao.utils.jax_utils import jax, jit, _jax_register_pytree_class, \
     _compute_sparsity, get_vmap_tangents, _update_subjac_sparsity, \
     _jax_derivs2partials, _ensure_returns_tuple, _update_add_input_kwargs, \
-    _update_add_output_kwargs, _re_init, _get_differentiable_compute_primal, _uncompress_jac, \
+    _update_add_output_kwargs, _re_init, _get_differentiable_compute_primal, \
     _jax2np, _compute_output_shapes
 from openmdao.utils.code_utils import get_return_names, get_function_deps
 import openmdao.utils.coloring as coloring_mod
@@ -375,7 +375,12 @@ class JaxImplicitComponent(ImplicitComponent):
             The partials to compute.
         """
         J = self._jac_func_(self._tangents['fwd'], tuple(chain(inputs.values(), outputs.values())))
-        partials.set_dense_jac(self, _uncompress_jac(self, _jax2np(J), 'fwd'))
+        J = _jax2np(J)
+        if self._coloring_info.coloring is not None:
+            J = self._coloring_info.coloring._expand_jac(J, 'fwd')
+            partials.set_csc_jac(self, J)
+        else:
+            partials.set_dense_jac(self, J)
 
     def _jacrev_colored(self, inputs, outputs, partials):
         """
@@ -391,7 +396,12 @@ class JaxImplicitComponent(ImplicitComponent):
             The partials to compute.
         """
         J = self._jac_func_(self._tangents['rev'], tuple(chain(inputs.values(), outputs.values())))
-        partials.set_dense_jac(self, _uncompress_jac(self, _jax2np(J).T, 'rev'))
+        J = _jax2np(J).T
+        if self._coloring_info.coloring is not None:
+            J = self._coloring_info.coloring._expand_jac(J, 'rev')
+            partials.set_csc_jac(self, J)
+        else:
+            partials.set_dense_jac(self, J)
 
     def compute_sparsity(self, direction=None, num_iters=1, perturb_size=1e-9):
         """
diff --git a/openmdao/core/component.py b/openmdao/core/component.py
index a56162fb6..b54548f94 100644
--- a/openmdao/core/component.py
+++ b/openmdao/core/component.py
@@ -4,7 +4,7 @@ import sys
 import inspect
 from collections import defaultdict
 from collections.abc import Iterable
-from itertools import product
+from itertools import product, chain
 from io import StringIO
 
 from numbers import Integral
@@ -78,7 +78,7 @@ class Component(System):
     _var_rel2meta : dict
         Dictionary mapping relative names to metadata.
         This is only needed while adding inputs and outputs. During setup, these are used to
-        build the dictionaries of metadata.
+        build the dictionaries of metadata.  This contains both continuous and discrete variables.
     _static_var_rel2meta : dict
         Static version of above - stores data for variables added outside of setup.
     _var_rel_names : {'input': [str, ...], 'output': [str, ...]}
@@ -1638,6 +1638,17 @@ class Component(System):
 
             self._subjacs_info[abs_key] = meta
 
+    def _column_iotypes(self):
+        """
+        Return a tuple of the iotypes that make up columns of the jacobian.
+
+        Returns
+        -------
+        tuple of the form ('input',)
+            The iotypes that make up columns of the jacobian.
+        """
+        return ('input',)
+
     def _get_partials_wrts(self):
         """
         Get list of 'wrt' variables that form the partial jacobian.
@@ -1647,11 +1658,8 @@ class Component(System):
         list
             List of 'wrt' relative variable names.
         """
-        # filter out any discrete inputs or outputs
-        if self._discrete_inputs:
-            return [n for n in self._var_rel_names['input'] if n not in self._discrete_inputs]
-
-        return list(self._var_rel_names['input'])
+        namelists = [self._var_rel_names[io] for io in self._column_iotypes()]
+        return [n for n in chain(*namelists)]
 
     def _get_partials_ofs(self, use_resname=False):
         """
@@ -1667,10 +1675,6 @@ class Component(System):
         list
             List of 'of' relative variable names.
         """
-        # filter out any discrete inputs or outputs
-        if self._discrete_outputs:
-            return [n for n in self._var_rel_names['output'] if n not in self._discrete_outputs]
-
         return list(self._var_rel_names['output'])
 
     def _matching_key_iter(self, of_patterns, wrt_patterns, use_resname=False):
@@ -1740,8 +1744,8 @@ class Component(System):
             List of tuples of the form (abs_name, meta) where abs_name is the absolute name of the
             matching variable and meta is the metadata for that variable.
         """
-        wrt_list = [pattern] if isinstance(pattern, str) else pattern
-        return [(pattern, find_matches(pattern, self._get_partials_wrts())) for pattern in wrt_list]
+        patterns = [pattern] if isinstance(pattern, str) else pattern
+        return [(pattern, find_matches(pattern, self._get_partials_wrts())) for pattern in patterns]
 
     def _check_partials_meta(self, abs_key, val, shape):
         """
diff --git a/openmdao/core/group.py b/openmdao/core/group.py
index e4986ffe2..6b43935e0 100644
--- a/openmdao/core/group.py
+++ b/openmdao/core/group.py
@@ -1906,8 +1906,9 @@ class Group(System):
                 continue
             try:
                 paths = [(i, m['path']) for i, m in enumerate(metalist) if not m['auto']]
-                top_origin = paths[0][1]
-                top_prom = metalist[paths[0][0]]['prom']
+                if paths:
+                    top_origin = paths[0][1]
+                    top_prom = metalist[paths[0][0]]['prom']
             except KeyError:
                 issue_warning("No auto IVCs found", prefix=self.msginfo, category=PromotionWarning)
             allmeta = set()
diff --git a/openmdao/core/implicitcomponent.py b/openmdao/core/implicitcomponent.py
index 1ff9ea6a2..06cd155d6 100644
--- a/openmdao/core/implicitcomponent.py
+++ b/openmdao/core/implicitcomponent.py
@@ -16,7 +16,6 @@ from openmdao.utils.units import simplify_unit
 from openmdao.utils.rangemapper import RangeMapper
 from openmdao.utils.om_warnings import issue_warning
 from openmdao.utils.coloring import _ColSparsityJac
-from openmdao.utils.iter_utils import meta2range_iter
 
 _tuplist = (tuple, list)
 
@@ -736,16 +735,16 @@ class ImplicitComponent(Component):
             raise ValueError(f"{self.msginfo}: residual units '{units}' for residual '{resid}' != "
                              f"output res_units '{res_units}' for output '{output}'.")
 
-    def _get_partials_wrts(self):
+    def _column_iotypes(self):
         """
-        Get the list of wrt variables that form the partial jacobian.
+        Return a tuple of the iotypes that make up columns of the jacobian.
 
         Returns
         -------
-        list
-            List of wrt variable names (relative names).
+        tuple of the form ('output', 'input')
+            The iotypes that make up columns of the jacobian.
         """
-        return list(self._var_rel_names['output']) + list(self._var_rel_names['input'])
+        return ('output', 'input')
 
     def _get_partials_ofs(self, use_resname=False):
         """
@@ -1010,31 +1009,6 @@ class ImplicitComponent(Component):
         return jac.get_sparsity()
 
 
-def _overlap_range_iter(meta_dict1, meta_dict2, names1=None, names2=None):
-    """
-    Yield names and ranges of overlapping variables from two metadata dictionaries.
-
-    The metadata dicts are assumed to contain a 'shape' entry, and the total size of the
-    variables in meta_dict1 must equal the total size of the variables in meta_dict2.
-    """
-    iter2 = meta2range_iter(meta_dict2, names=names2)
-    start2 = end2 = -1
-
-    for name1, start1, end1 in meta2range_iter(meta_dict1, names=names1):
-        try:
-            while not (start2 <= start1 < end2 or start2 <= end1 < end2):
-                name2, start2, end2 = next(iter2)
-
-            if end1 < end2:
-                yield name1, start1, end1, name2, start2, end2
-            else:
-                while end1 >= end2:
-                    yield name1, start1, end1, name2, start2, end2
-                    name2, start2, end2 = next(iter2)
-        except StopIteration:
-            return
-
-
 class _ResidsWrapper(object):
     def __init__(self, vec, name2slice_shape):
         self.__dict__['_vec'] = vec
diff --git a/openmdao/core/problem.py b/openmdao/core/problem.py
index 6cceb746e..74ad284f0 100644
--- a/openmdao/core/problem.py
+++ b/openmdao/core/problem.py
@@ -1085,8 +1085,9 @@ class Problem(object, metaclass=ProblemMetaclass):
             first = True
             self._metadata['static_mode'] = False
             try:
-                model._setup_part2()
-                self._check_collected_errors()
+                if self._metadata['setup_status'] < _SetupStatus.POST_SETUP2:
+                    model._setup_part2()
+                    self._check_collected_errors()
 
                 responses = model.get_responses(recurse=True, use_prom_ivc=True)
                 designvars = model.get_design_vars(recurse=True, use_prom_ivc=True)
@@ -1180,24 +1181,23 @@ class Problem(object, metaclass=ProblemMetaclass):
         **setup_kwargs : dict
             Keyword arguments to pass to the setup method if it hasn't already been called.
         """
-        current_status = self._metadata['setup_status']
-        if current_status >= status:
+        if self._metadata['setup_status'] >= status:
             return
 
         if status >= _SetupStatus.POST_SETUP:
-            if current_status < _SetupStatus.POST_SETUP:
+            if self._metadata['setup_status'] < _SetupStatus.POST_SETUP:
                 self.setup(**setup_kwargs)
-                current_status = _SetupStatus.POST_SETUP
+                self._metadata['setup_status'] = _SetupStatus.POST_SETUP
 
         if status >= _SetupStatus.POST_SETUP2:
-            if current_status < _SetupStatus.POST_SETUP2:
+            if self._metadata['setup_status'] < _SetupStatus.POST_SETUP2:
                 self.model._setup_part2()
-                current_status = _SetupStatus.POST_SETUP2
+                self._metadata['setup_status'] = _SetupStatus.POST_SETUP2
 
         if status >= _SetupStatus.POST_FINAL_SETUP:
-            if current_status < _SetupStatus.POST_FINAL_SETUP:
+            if self._metadata['setup_status'] < _SetupStatus.POST_FINAL_SETUP:
                 self.final_setup()
-                current_status = _SetupStatus.POST_FINAL_SETUP
+                self._metadata['setup_status'] = _SetupStatus.POST_FINAL_SETUP
 
     def check_partials(self, out_stream=_DEFAULT_OUT_STREAM, includes=None, excludes=None,
                        compact_print=False, abs_err_tol=0.0, rel_err_tol=1e-6,
diff --git a/openmdao/jacobians/jacobian.py b/openmdao/jacobians/jacobian.py
index 9ae4df6ee..75178c287 100644
--- a/openmdao/jacobians/jacobian.py
+++ b/openmdao/jacobians/jacobian.py
@@ -385,6 +385,30 @@ class Jacobian(object):
                     if match_inds.size > 0:
                         subjac['val'][match_inds] = column[start:end][subjac['rows'][match_inds]]
 
+    def set_csc_jac(self, system, jac):
+        """
+        Assign a CSC jacobian to this jacobian.
+
+        Parameters
+        ----------
+        system : System
+            The system that owns this jacobian.
+        jac : csc_matrix
+            CSC jacobian.
+        """
+        ofiter = list(system._jac_of_iter())
+        for wrt, wstart, wend, _, _, _ in system._jac_wrt_iter():
+            wjac = jac[:, wstart:wend]
+            for of, start, end, _, _ in ofiter:
+                key = (of, wrt)
+                if key in self._subjacs_info:
+                    subjac = self.get_metadata(key)
+                    if subjac['cols'] is None:  # dense
+                        subjac['val'][:, :] = wjac[start:end, :]
+                    else:  # our COO format
+                        subj = wjac[start:end, :]
+                        subjac['val'][:] = subj[subjac['rows'], subjac['cols']]
+
     def set_dense_jac(self, system, jac):
         """
         Assign a dense jacobian to this jacobian.
diff --git a/openmdao/utils/coloring.py b/openmdao/utils/coloring.py
index 828afb09a..f66cb08c6 100644
--- a/openmdao/utils/coloring.py
+++ b/openmdao/utils/coloring.py
@@ -699,6 +699,8 @@ class Coloring(object):
         Dictionary mapping absolute names to promoted names.
     _subtractions : list
         List of subtraction tuples. Only used for the substitution method of bidirectional coloring.
+    _color_arrays : dict
+        Dictionary of color arrays for each direction.
     """
 
     def __init__(self, sparsity, row_vars=None, row_var_sizes=None, col_vars=None,
@@ -735,6 +737,7 @@ class Coloring(object):
 
         self._abs2prom = None
         self._subtractions = None
+        self._color_arrays = {}
 
     def get_renamed_copy(self, row_translate, col_translate):
         """
@@ -787,6 +790,40 @@ class Coloring(object):
 
         yield from colors
 
+    def _get_color_array(self, direction):
+        """
+        Return the color array for the given direction.
+
+        The size of the array is the number of nonzero values in the sparsity matrix.
+
+        Parameters
+        ----------
+        direction : str
+            Derivative direction ('fwd' or 'rev').
+
+        Returns
+        -------
+        ndarray
+            Color array for the given direction.
+        """
+        if direction not in self._color_arrays:
+            if direction == 'fwd':
+                color_array = np.zeros(self._shape[1], dtype=int)
+                for i, color_list in enumerate(self._fwd[0]):
+                    color_array[color_list] = i
+                colors_nz = color_array[self._nzcols]
+            elif direction == 'rev':
+                color_array = np.zeros(self._shape[0], dtype=int)
+                for i, color_list in enumerate(self._rev[0]):
+                    color_array[color_list] = i
+                colors_nz = color_array[self._nzrows]
+            else:
+                raise RuntimeError(f"Invalid direction '{direction}' in _get_color_array.")
+
+            self._color_arrays[direction] = colors_nz
+
+        return self._color_arrays[direction]
+
     def color_nonzero_iter(self, direction):
         """
         Given a direction, yield an iterator over (columns, nz_rows) or (rows, nz_columns).
@@ -841,9 +878,9 @@ class Coloring(object):
                 for jac_irow, nzpart in zip(nzs, nzparts):
                     yield compressed_j[i, nzpart], nzpart, jac_irow
 
-    def expand_jac(self, compressed_j, direction):
+    def _expand_jac(self, compressed_j, direction):
         """
-        Expand the given compressed jacobian into a full jacobian.
+        Expand the given compressed dense jacobian into a full CSC jacobian.
 
         Parameters
         ----------
@@ -854,11 +891,16 @@ class Coloring(object):
 
         Returns
         -------
-        ndarray
+        csc_matrix
             The full jacobian.
         """
-        sparsity = self.get_sparsity_of_type(bool).toarray()
-        return colored2full_jac(compressed_j, sparsity, self.color_iter(direction), direction)
+        colors_coo = self._get_color_array(direction)
+        if direction == 'fwd':
+            data = compressed_j[self._nzrows, colors_coo]
+        elif direction == 'rev':
+            data = compressed_j[colors_coo, self._nzcols]
+
+        return csc_matrix((data, (self._nzrows, self._nzcols)), shape=self._shape)
 
     def get_row_col_map(self, direction):
         """
@@ -2209,105 +2251,6 @@ class Coloring(object):
         return _sort_subtractions(subtractions)
 
 
-def _get_colormap(color_array):
-    """
-    Get a map of colors to column indices.
-
-    Parameters
-    ----------
-    color_array : ndarray
-        Array of colors.
-
-    Returns
-    -------
-    dict
-        Map of colors to column indices.
-    """
-    colmap = {}
-    for i, c in enumerate(color_array):
-        if c not in colmap:
-            colmap[c] = []
-        colmap[c].append(i)
-    return colmap
-
-
-def _get_expander_matrix(sparsity_shape, color_iter, direction):
-    """
-    Get a matrix to expand a colored jacobian to an expanded colored jacobian.
-
-    The expanded colored jacobian is the shape of the full jacobian but columns sharing
-    the same color are summed together.
-
-    Parameters
-    ----------
-    sparsity_shape : tuple
-        Shape of the sparsity pattern of the jacobian.
-    color_iter : iterator
-        Iterator over columns for each color.
-    direction : str
-        Direction of the jacobian to expand.
-
-    Returns
-    -------
-    csc_matrix
-        Matrix to expand a colored jacobian to an expanded colored jacobian.
-    """
-    # Ic looks like an identity matrix (ncols x nncols for fwd) with all columns corresponding
-    # to the same color being summed together and combined into a single column.
-    # Ic = csc_matrix((np.array([]), (np.array([], dtype=np.int), np.array([], dtype=np.int))),
-    #                 shape=(I.shape[1], np.unique(color_array).size))
-    Ic_rows = []
-    Ic_cols = []
-    ncolors = 0
-    for c, cols in enumerate(color_iter):
-        ncolors += 1
-        for col in cols:
-            Ic_cols.append(c)
-            Ic_rows.append(col)
-
-    if direction == 'rev':
-        Ic_cols, Ic_rows = Ic_rows, Ic_cols
-        nrows = ncolors
-        ncolors = sparsity_shape[0]
-    else:
-        nrows = sparsity_shape[1]
-
-    Ic = coo_matrix((np.ones(len(Ic_rows)), (Ic_rows, Ic_cols)),
-                    shape=(nrows, ncolors)).T
-
-    if direction == 'fwd':
-        return Ic.tocsc()
-    else:
-        return Ic.tocsr()
-
-
-def colored2full_jac(Jc, sparsity, color_iter, direction):
-    """
-    Expand a colored jacobian to a full jacobian.
-
-    Parameters
-    ----------
-    Jc : csc_matrix
-        Colored jacobian.
-    sparsity : csc_matrix
-        Sparsity pattern of the jacobian.
-    color_iter : iterator
-        Iterator over columns for each color.
-    direction : str
-        Direction of the jacobian to expand.
-
-    Returns
-    -------
-    csc_matrix
-        Full jacobian.
-    """
-    X = _get_expander_matrix(sparsity.shape, color_iter, direction).toarray()
-    if direction == 'fwd':
-        return np.where(sparsity, Jc @ X, 0.)
-    else:
-        return np.where(sparsity, X @ Jc, 0.)
-
-
 def _order_by_ID(col_adj_matrix):
     """
     Return columns in order of incidence degree (ID).
@@ -2335,7 +2278,7 @@ def _order_by_ID(col_adj_matrix):
     colored_degrees = np.zeros(ncols, dtype=INT_DTYPE)
     colored_degrees[col_adj_matrix.indices] = 1  # make sure zero cols aren't considered
 
-    for i in range(np.nonzero(colored_degrees)[0].size):
+    for _ in range(np.nonzero(colored_degrees)[0].size):
         col = colored_degrees.argmax()
         colnzrows = col_adj_matrix.getcol(col).indices
         colored_degrees[colnzrows] += 1
@@ -2599,7 +2542,8 @@ def _color_partition(J, Jprows, Jpcols, direct=True, overlap=None):
     Jpcols : ndarray
         Nonzero columns of a partition of the matrix being colored.
     direct : bool
-        If True, use the direct method to compute the column adjacency matrix.
+        If True, use the direct method to compute the column adjacency matrix. Otherwise, use the
+        substitution method.
     overlap : set or None
         Nonzero entries indicate overlapping columns.
 
diff --git a/openmdao/utils/jax_utils.py b/openmdao/utils/jax_utils.py
index 08a1e5ece..48245b72d 100644
--- a/openmdao/utils/jax_utils.py
+++ b/openmdao/utils/jax_utils.py
@@ -1169,29 +1169,6 @@ def _jax2np(J):
         return np.asarray(J).reshape(-1, J.shape[-1])
 
 
-def _uncompress_jac(self, J, direction):
-    """
-    Uncompress the Jacobian using the coloring information.
-
-    Parameters
-    ----------
-    self : Component
-        The component to uncompress the Jacobian for.
-    J : ndarray
-        The Jacobian to uncompress.
-    direction : str
-        The direction to uncompress the Jacobian in.
-
-    Returns
-    -------
-    ndarray
-        The uncompressed Jacobian.
-    """
-    if self._coloring_info.coloring is not None:
-        return self._coloring_info.coloring.expand_jac(J, direction)
-    return J
-
-
 def _jax_derivs2partials(self, deriv_vals, partials, ofnames, wrtnames):
     """
     Copy JAX derivatives into partials.


Test patch:
diff --git a/openmdao/core/tests/test_group.py b/openmdao/core/tests/test_group.py
index c5726cecc..5afea3431 100644
--- a/openmdao/core/tests/test_group.py
+++ b/openmdao/core/tests/test_group.py
@@ -2034,6 +2034,35 @@ class TestGroupPromotes(unittest.TestCase):
         self.assertEqual(str(cm.exception),
                            "\"<model> <class BranchGroup>: Could not find 'Branch1.G1.comp1.a'. Perhaps you meant one of the following variables: ['Branch1.G1.comp1.b', 'Branch1.G1.a']\"")
 
+    def test_promotes_ivc_bug(self):
+        # This used to fail with an index error
+
+        prob = om.Problem()
+        model = prob.model
+        model.add_subsystem("ivc", om.IndepVarComp(), promotes=["*"])
+        model.ivc.add_output("in1", [1.,1.])
+        model.add_subsystem("comp1", om.ExecComp("out1 = in1 * 2.0"))
+        model.add_subsystem("comp2", om.ExecComp("out1 = in1 * 3.0"))
+
+        model.connect("in1", "comp1.in1", src_indices=[0])
+        model.connect("in1", "comp2.in1", src_indices=[1])
+
+        prob.setup(mode='rev')
+        om.n2(prob, show_browser=False, outfile='n2.html')
+        prob.run_model()
+
+    def test_promotes_auto_ivc_bug(self):
+        # This used to fail with an index error
+
+        prob = om.Problem()
+        model = prob.model
+        model.add_subsystem("comp1", om.ExecComp("out1 = in1 * 2.0"))
+        model.add_subsystem("comp2", om.ExecComp("out1 = in1 * 3.0"))
+
+        prob.setup(mode='rev')
+        om.n2(prob, show_browser=False, outfile='n2.html')
+        prob.run_model()
+
     def test_multiple_promotes_collision(self):
 
         class SubGroup(om.Group):