diff --git a/.github/workflows/unit_tests.yml b/.github/workflows/unit_tests.yml
index 9809d1b54..af5d2b843 100644
--- a/.github/workflows/unit_tests.yml
+++ b/.github/workflows/unit_tests.yml
@@ -128,8 +128,9 @@ jobs:
         run: |
           # Install petsc for openmdao tests (there's a problem with 3.22.3 on mac runner)
           mamba install -c conda-forge "petsc>=3.19,!=3.22.3" petsc4py "numpy>=2" -q -y;
+          pip list
           # Install pygeo for mach tests
-          mamba install -c smdogroup pygeo -q -y;
+          mamba install -c smdogroup pygeo "numpy>=2" -q -y;
           # Install ESP/CAPS prebuilt for shape derivatives, only in Real mode
            if [[ ${{ matrix.NAME }} == 'Real-Ubuntu' ]]; then
             cd $TACS_DIR/extern/
@@ -159,6 +160,7 @@ jobs:
           if [[ ${{ matrix.NAME }} == 'Real-Ubuntu' ]] && [[ -f ${TACS_DIR}/extern/ESP122/ESPenv.sh ]]; then
             source ${TACS_DIR}/extern/ESP122/ESPenv.sh
           fi
+          pip list
           testflo -v --timeout 420 .;
       - name: Build docs
         run: |
diff --git a/tacs/mach/struct_problem.py b/tacs/mach/struct_problem.py
index 7244d17d1..cc5cf4b48 100644
--- a/tacs/mach/struct_problem.py
+++ b/tacs/mach/struct_problem.py
@@ -1524,6 +1524,11 @@ def readExternalForceFile(self, fileName):
         # Step 2: Overwrite bdfInfo loads with forceInfo loads
         bdfInfo.loads = copy.deepcopy(forceInfo.loads)
         bdfInfo.load_combinations = copy.deepcopy(forceInfo.load_combinations)
+        # The newly assigned loads are not cross-referenced (node_ref etc. are None),
+        # so force a re-cross-reference on the next addLoadFromBDF call. bdfInfo is
+        # shared across all problems from the same meshLoader, so it may already be
+        # marked as cross-referenced from a previously constructed problem.
+        bdfInfo.is_xrefed = False
 
         # Create a copy of the internal loads already added to model
         F = self.staticProblem.F
@@ -1540,3 +1545,6 @@ def readExternalForceFile(self, fileName):
         # Step 3: Restore original loads back into bdfInfo
         bdfInfo.loads = originalLoads
         bdfInfo.load_combinations = originalLoadCombinations
+        # The restored loads are an un-cross-referenced deepcopy, so force a
+        # re-cross-reference on the next addLoadFromBDF call.
+        bdfInfo.is_xrefed = False
diff --git a/tests/integration_tests/test_mach_beam_load_file.py b/tests/integration_tests/test_mach_beam_load_file.py
index 95fa355f0..84195896d 100644
--- a/tests/integration_tests/test_mach_beam_load_file.py
+++ b/tests/integration_tests/test_mach_beam_load_file.py
@@ -31,10 +31,16 @@ class TestMACHBeamExample(MACHStructProblemTestCase.MACHStructProblemTest):
 
     N_PROCS = 2
 
-    # Reference values for regression testing
+    # Reference values for regression testing.
+    # The second problem is identical to the first (same load file, same assembler),
+    # so it should produce the same function values. It exists to guard against a
+    # regression where loading forces from a file only worked for the first
+    # StructProblem sharing a given assembler/meshLoader (see readExternalForceFile).
     FUNC_REFS = {
         "tip_shear_ks_vmfailure": 2.492124644887184,
         "tip_shear_mass": 2.7799999999999963,
+        "tip_shear_2_ks_vmfailure": 2.492124644887184,
+        "tip_shear_2_mass": 2.7799999999999963,
     }
 
     def setup_struct_problems(self, comm):
@@ -81,4 +87,18 @@ def element_callback(
         # Create MACH StructProblem
         structProb = StructProblem(staticProb, FEAAssembler, loadFile=load_file)
 
-        return [structProb]
+        # Create a second, identical static problem from the same assembler and
+        # wrap it in another StructProblem that also loads forces from a file.
+        # Because the bdfInfo is shared across problems from the same meshLoader,
+        # this exercises the code path where the load file must be applied to more
+        # than one StructProblem.
+        staticProb2 = FEAAssembler.createStaticProblem("tip_shear_2")
+        staticProb2.addFunction("mass", functions.StructuralMass)
+        staticProb2.addFunction(
+            "ks_vmfailure", functions.KSFailure, safetyFactor=1.0, ksWeight=ksweight
+        )
+        staticProb2.setOption("L2Convergence", 1e-20)
+        staticProb2.setOption("L2ConvergenceRel", 1e-20)
+        structProb2 = StructProblem(staticProb2, FEAAssembler, loadFile=load_file)
+
+        return [structProb, structProb2]