diff --git a/test/conftest.py b/test/conftest.py
index 60b9dfa9b4..ec03bfc6c6 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -17,8 +17,10 @@
 import importlib
 import os
 import pathlib
+import random
 from collections import defaultdict
 
+import numpy as np
 import pytest
 import torch
 
@@ -171,3 +173,26 @@ def requires_module(names):
 def device(request):
     """Device fixture that automatically skips CUDA tests when not available."""
     return request.param
+
+
+@pytest.fixture(autouse=True, scope="function")
+def seed_random_state():
+    """Reset all random number generators to a fixed seed before each test.
+
+    This ensures test reproducibility and isolation - each test starts with
+    identical RNG state regardless of test execution order or subset.
+
+    Tests that need a specific seed can still call torch.manual_seed() etc.
+    explicitly, which will override this fixture's seeding.
+    """
+    SEED = 95051
+
+    random.seed(SEED)
+    np.random.seed(SEED)
+    torch.manual_seed(SEED)
+
+    # CUDA seeding (no-op if CUDA unavailable)
+    if torch.cuda.is_available():
+        torch.cuda.manual_seed_all(SEED)
+
+    yield
diff --git a/test/mesh/io/io_pyvista/test_data_array_shapes.py b/test/mesh/io/io_pyvista/test_data_array_shapes.py
index 504f8b9a52..937b3a2336 100644
--- a/test/mesh/io/io_pyvista/test_data_array_shapes.py
+++ b/test/mesh/io/io_pyvista/test_data_array_shapes.py
@@ -30,7 +30,6 @@ class TestDataArrayShapes:
 
     def test_scalar_data(self):
         """Test scalar data (1D array per point/cell)."""
-        np.random.seed(0)
         pv_mesh = pv.Sphere(radius=1.0, theta_resolution=10, phi_resolution=10)
 
         # Add scalar data
@@ -53,7 +52,6 @@ def test_scalar_data(self):
 
     def test_vector_data(self):
         """Test vector data (Nx3 arrays)."""
-        np.random.seed(0)
         pv_mesh = pv.Sphere(radius=1.0, theta_resolution=10, phi_resolution=10)
 
         # Add vector data
@@ -81,7 +79,6 @@ def test_matrix_data(self):
         For higher-dimensional data like 3x3 stress tensors, you must
         flatten them to (n, 9) before adding to PyVista.
         """
-        np.random.seed(0)
         pv_mesh = pv.Sphere(radius=1.0, theta_resolution=10, phi_resolution=10)
 
         # For tensor data, must be pre-flattened to 2D
@@ -113,7 +110,6 @@ def test_large_2d_array_data(self):
         NOTE: PyVista only accepts arrays with dimensionality ≤ 2.
         Higher-order tensors must be pre-flattened before adding to PyVista.
         """
-        np.random.seed(0)
         pv_mesh = pv.Sphere(radius=1.0, theta_resolution=10, phi_resolution=10)
 
         # For higher-dimensional data, flatten to 2D before adding to PyVista
@@ -137,7 +133,6 @@ def test_large_2d_array_data(self):
 
     def test_mixed_data_types(self):
         """Test mesh with multiple data arrays of different shapes and types."""
-        np.random.seed(0)
         pv_mesh = pv.Sphere(radius=1.0, theta_resolution=10, phi_resolution=10)
 
         # Clear default data to have a clean slate
diff --git a/test/mesh/io/io_pyvista/test_data_preservation.py b/test/mesh/io/io_pyvista/test_data_preservation.py
index 8c0d1f8b6c..996c38fb83 100644
--- a/test/mesh/io/io_pyvista/test_data_preservation.py
+++ b/test/mesh/io/io_pyvista/test_data_preservation.py
@@ -49,7 +49,6 @@ class TestDataPreservation:
 
     def test_point_data_preserved(self):
         """Test that point_data is preserved during conversion."""
-        np.random.seed(0)
         pv_mesh = pv.Sphere()
 
         # Explicitly create point data
@@ -78,7 +77,6 @@ def test_point_data_preserved(self):
 
     def test_cell_data_preserved(self):
         """Test that cell_data is preserved as cell_data."""
-        np.random.seed(0)
         pv_mesh = pv.Sphere()
 
         # Explicitly create cell data
@@ -161,7 +159,6 @@ class TestDataPreservationParametrized:
 
     def test_data_preservation_with_device_transfer(self, device):
         """Test that data is preserved when transferring to different device."""
-        np.random.seed(42)
         pv_mesh = pv.Sphere(theta_resolution=5, phi_resolution=5)
         pv_mesh.point_data["temp"] = np.random.rand(pv_mesh.n_points).astype(np.float32)
         pv_mesh.cell_data["pressure"] = np.random.rand(pv_mesh.n_cells).astype(
diff --git a/test/mesh/io/io_pyvista/test_from_pyvista_0d.py b/test/mesh/io/io_pyvista/test_from_pyvista_0d.py
index f49d5fa575..ed37d7de7f 100644
--- a/test/mesh/io/io_pyvista/test_from_pyvista_0d.py
+++ b/test/mesh/io/io_pyvista/test_from_pyvista_0d.py
@@ -30,7 +30,6 @@ class TestFromPyvista0D:
 
     def test_pointset_auto_detection(self):
         """Test automatic detection of 0D manifold from PointSet."""
-        np.random.seed(0)
         points = np.random.rand(100, 3).astype(np.float32)
         pv_mesh = pv.PointSet(points)
 
@@ -50,7 +49,6 @@ def test_pointset_auto_detection(self):
 
     def test_pointset_explicit_dim(self):
         """Test explicit manifold_dim specification for point cloud."""
-        np.random.seed(0)
         points = np.random.rand(50, 3).astype(np.float32)
         pv_mesh = pv.PointSet(points)
 
@@ -65,7 +63,6 @@ def test_polydata_points_only(self):
 
         PolyData can represent point clouds using vertex cells.
         """
-        np.random.seed(0)
         points = np.random.rand(25, 3).astype(np.float32)
         pv_mesh = pv.PolyData(points)
 
diff --git a/test/mesh/io/io_pyvista/test_from_pyvista_1d.py b/test/mesh/io/io_pyvista/test_from_pyvista_1d.py
index 94fd5881d0..bdba5f5999 100644
--- a/test/mesh/io/io_pyvista/test_from_pyvista_1d.py
+++ b/test/mesh/io/io_pyvista/test_from_pyvista_1d.py
@@ -105,7 +105,6 @@ def test_spline_constructed(self):
 
         Create a spline through specific points and verify it converts correctly.
         """
-        np.random.seed(0)
         # Create control points for the spline
         control_points = np.array(
             [
diff --git a/test/mesh/io/io_pyvista/test_round_trip.py b/test/mesh/io/io_pyvista/test_round_trip.py
index c8654d810c..4ec5d7bce9 100644
--- a/test/mesh/io/io_pyvista/test_round_trip.py
+++ b/test/mesh/io/io_pyvista/test_round_trip.py
@@ -94,7 +94,6 @@ def test_round_trip_1d_spline(self):
 
     def test_round_trip_0d_pointset(self):
         """Test round-trip conversion for 0D mesh."""
-        np.random.seed(0)
         points_orig = np.random.rand(25, 3).astype(np.float32)
         pv_original = pv.PointSet(points_orig)
 
@@ -108,7 +107,6 @@ def test_round_trip_0d_pointset(self):
 
     def test_round_trip_with_data(self):
         """Test round-trip conversion preserves data arrays."""
-        np.random.seed(0)
         pv_original = pv.Sphere(theta_resolution=10, phi_resolution=10)
         pv_original.clear_data()
 
@@ -225,7 +223,6 @@ def test_round_trip_spline_device_parametrized(self, device):
 
     def test_device_transfer_preserves_data(self, device):
         """Test that device transfer preserves all data."""
-        np.random.seed(42)
         # Create mesh with data
         pv_mesh = pv.Sphere(theta_resolution=5, phi_resolution=5)
         pv_mesh.point_data["temp"] = np.random.rand(pv_mesh.n_points).astype(np.float32)
diff --git a/test/mesh/io/io_pyvista/test_to_pyvista.py b/test/mesh/io/io_pyvista/test_to_pyvista.py
index a4f4784fa5..b6b66f6fb2 100644
--- a/test/mesh/io/io_pyvista/test_to_pyvista.py
+++ b/test/mesh/io/io_pyvista/test_to_pyvista.py
@@ -102,7 +102,6 @@ def test_1d_mesh_to_polydata(self):
 
     def test_0d_mesh_to_pointset(self):
         """Test converting 0D mesh to PointSet."""
-        np.random.seed(0)
         points = torch.from_numpy(np.random.rand(50, 3).astype(np.float32))
         cells = torch.empty((0, 1), dtype=torch.long)
 
@@ -116,9 +115,6 @@ def test_0d_mesh_to_pointset(self):
 
     def test_data_preservation_to_pyvista(self):
         """Test that point_data, cell_data, and global_data are preserved."""
-        np.random.seed(0)
-        torch.manual_seed(42)
-
         # Create a mesh with data
         points = torch.rand(10, 3)
         cells = torch.tensor([[0, 1, 2], [2, 3, 4]], dtype=torch.long)
diff --git a/test/mesh/visualization/test_visualization.py b/test/mesh/visualization/test_visualization.py
index eddc3cd2b4..c38ab0e504 100644
--- a/test/mesh/visualization/test_visualization.py
+++ b/test/mesh/visualization/test_visualization.py
@@ -215,7 +215,6 @@ def test_pyvista_points_padded_to_3d():
 
 def test_unsupported_spatial_dims():
     """Test that meshes with >3 spatial dimensions raise error."""
-    torch.manual_seed(42)
     # Create a 4D mesh
     points = torch.randn(10, 4)
     cells = torch.randint(0, 10, (5, 2))
@@ -241,7 +240,6 @@ def test_no_scalars():
 
 def test_point_scalars_tensor():
     """Test point scalars with direct tensor."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
     point_scalars = torch.rand(mesh.n_points)
     ax = mesh.draw(show=False, backend="matplotlib", point_scalars=point_scalars)
@@ -251,7 +249,6 @@ def test_point_scalars_tensor():
 
 def test_cell_scalars_tensor():
     """Test cell scalars with direct tensor."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
     cell_scalars = torch.rand(mesh.n_cells)
     ax = mesh.draw(show=False, backend="matplotlib", cell_scalars=cell_scalars)
@@ -261,7 +258,6 @@ def test_cell_scalars_tensor():
 
 def test_point_scalars_key():
     """Test point scalars with key lookup."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
     mesh.point_data["temperature"] = torch.rand(mesh.n_points)
     ax = mesh.draw(show=False, backend="matplotlib", point_scalars="temperature")
@@ -271,7 +267,6 @@ def test_point_scalars_key():
 
 def test_cell_scalars_key():
     """Test cell scalars with key lookup."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
     mesh.cell_data["pressure"] = torch.rand(mesh.n_cells)
     ax = mesh.draw(show=False, backend="matplotlib", cell_scalars="pressure")
@@ -283,7 +278,6 @@ def test_nested_tensordict_key():
     """Test scalar lookup with nested TensorDict key."""
     from tensordict import TensorDict
 
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
 
     # Create nested structure
@@ -300,7 +294,6 @@ def test_nested_tensordict_key():
 
 def test_multidimensional_scalars_norm():
     """Test that multidimensional scalars are L2-normed."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
 
     # Create 3D vector field
@@ -313,7 +306,6 @@ def test_multidimensional_scalars_norm():
 
 def test_mutual_exclusivity():
     """Test that point_scalars and cell_scalars are mutually exclusive."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
 
     with pytest.raises(ValueError, match="mutually exclusive"):
@@ -326,7 +318,6 @@ def test_mutual_exclusivity():
 
 def test_scalar_wrong_shape():
     """Test that scalars with wrong shape raise error."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
 
     with pytest.raises(ValueError, match="wrong first dimension"):
@@ -350,7 +341,6 @@ def test_scalar_key_not_found():
 
 def test_colormap():
     """Test custom colormap."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
     mesh.cell_data["data"] = torch.rand(mesh.n_cells)
 
@@ -361,7 +351,6 @@ def test_colormap():
 
 def test_vmin_vmax():
     """Test colormap range specification."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
     mesh.cell_data["data"] = torch.rand(mesh.n_cells)
 
@@ -447,7 +436,6 @@ def test_draw_1d_in_2d():
 
 def test_draw_empty_mesh():
     """Test drawing mesh with no cells."""
-    torch.manual_seed(42)
     points = torch.randn(10, 2)
     cells = torch.empty((0, 3), dtype=torch.long)
     mesh = Mesh(points=points, cells=cells)
@@ -459,7 +447,6 @@ def test_draw_empty_mesh():
 
 def test_pyvista_with_scalars():
     """Test PyVista backend with scalar coloring."""
-    torch.manual_seed(42)
     mesh = create_3d_surface_mesh()
     mesh.cell_data["pressure"] = torch.rand(mesh.n_cells)
 
@@ -472,7 +459,6 @@ def test_pyvista_with_scalars():
 
 def test_pyvista_with_point_scalars():
     """Test PyVista backend with point scalar coloring."""
-    torch.manual_seed(42)
     mesh = create_3d_surface_mesh()
     mesh.point_data["temperature"] = torch.rand(mesh.n_points)
 
@@ -488,7 +474,6 @@ def test_pyvista_with_point_scalars():
 
 def test_full_workflow_matplotlib():
     """Test complete workflow with matplotlib backend."""
-    torch.manual_seed(42)
     mesh = create_2d_triangle_mesh()
 
     # Add some data
@@ -512,7 +497,6 @@ def test_full_workflow_matplotlib():
 
 def test_full_workflow_pyvista():
     """Test complete workflow with PyVista backend."""
-    torch.manual_seed(42)
     mesh = create_3d_surface_mesh()
 
     # Add some data
@@ -608,7 +592,6 @@ def test_basic_visualization_parametrized(
     @pytest.mark.parametrize("backend", ["matplotlib", "pyvista"])
     def test_visualization_with_scalars_parametrized(self, backend):
         """Test visualization with scalar data across backends."""
-        torch.manual_seed(42)
         if backend == "pyvista":
             # Use 3D mesh for PyVista
             points = torch.tensor([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
diff --git a/test/models/domino/test_domino_encodings.py b/test/models/domino/test_domino_encodings.py
index f3a1485993..089d2401a1 100644
--- a/test/models/domino/test_domino_encodings.py
+++ b/test/models/domino/test_domino_encodings.py
@@ -28,8 +28,6 @@ def test_fourier_mlp(device, fourier_features, num_modes):
     """Test FourierMLP with various configurations"""
     from physicsnemo.nn import FourierMLP
 
-    torch.manual_seed(0)
-
     model = FourierMLP(
         input_features=3,
         base_layer=64,
@@ -48,8 +46,6 @@ def test_fourier_encode_vectorized(device):
     """Test fourier encoding function"""
     from physicsnemo.nn import fourier_encode
 
-    torch.manual_seed(0)
-
     coords = torch.randn(4, 20, 3).to(device)
     freqs = torch.exp(torch.linspace(0, math.pi, 5)).to(device)
 
@@ -66,8 +62,6 @@ def test_local_geometry_encoding(device):
 
     BATCH_SIZE = 1
 
-    torch.manual_seed(0)
-
     N_ENCODING_CHANNELS = 3
     N_NEIGHBORS = 32
     N_MESH_POINTS = 50
@@ -99,8 +93,6 @@ def test_multi_geometry_encoding(device, geo_encoding_type):
     from physicsnemo.models.domino.encodings import MultiGeometryEncoding
     from physicsnemo.models.domino.model import get_activation
 
-    torch.manual_seed(0)
-
     BATCH_SIZE = 1
     N_MESH_POINTS = 50
     GRID_RESOLUTION = (32, 32, 32)
diff --git a/test/models/domino/test_domino_geometry_rep.py b/test/models/domino/test_domino_geometry_rep.py
index c1e2208059..c8502e9284 100644
--- a/test/models/domino/test_domino_geometry_rep.py
+++ b/test/models/domino/test_domino_geometry_rep.py
@@ -29,8 +29,6 @@ def test_geo_conv_out(device, act, fourier_features):
     """Test GeoConvOut layer"""
     from physicsnemo.models.domino.geometry_rep import GeoConvOut
 
-    torch.manual_seed(0)
-
     @dataclass
     class TestParams:
         base_neurons: int = 32
@@ -71,8 +69,6 @@ def test_geo_processor(device, act):
     """Test GeoProcessor CNN"""
     from physicsnemo.models.domino.geometry_rep import GeoProcessor
 
-    torch.manual_seed(0)
-
     @dataclass
     class TestParams:
         base_filters: int = 8
@@ -98,8 +94,6 @@ def test_geometry_rep(
     """Test GeometryRep module with different configurations"""
     from physicsnemo.models.domino.geometry_rep import GeometryRep
 
-    torch.manual_seed(0)
-
     # Modify params for this test
     params = base_model_params()
     params.geometry_encoding_type = geometry_encoding_type
diff --git a/test/models/domino/test_domino_mlps.py b/test/models/domino/test_domino_mlps.py
index 0d22d35753..a7cb2f76bc 100644
--- a/test/models/domino/test_domino_mlps.py
+++ b/test/models/domino/test_domino_mlps.py
@@ -26,8 +26,6 @@ def test_aggregation_model(device, activation):
     from physicsnemo.models.domino.mlps import AggregationModel
     from physicsnemo.models.domino.model import get_activation
 
-    torch.manual_seed(0)
-
     model = AggregationModel(
         input_features=100,
         output_features=1,
@@ -47,8 +45,6 @@ def test_local_point_conv(device, activation):
     from physicsnemo.models.domino.mlps import LocalPointConv
     from physicsnemo.models.domino.model import get_activation
 
-    torch.manual_seed(0)
-
     model = LocalPointConv(
         input_features=50,
         base_layer=128,
diff --git a/test/models/domino/test_domino_solutions.py b/test/models/domino/test_domino_solutions.py
index a6eba70a1b..2ac04281c6 100644
--- a/test/models/domino/test_domino_solutions.py
+++ b/test/models/domino/test_domino_solutions.py
@@ -32,8 +32,6 @@ def test_solution_calculator_volume(
     from physicsnemo.models.domino.solutions import SolutionCalculatorVolume
     from physicsnemo.nn import FourierMLP, get_activation
 
-    torch.manual_seed(0)
-
     activation = get_activation("relu")
 
     # Create parameter model if needed
@@ -116,8 +114,6 @@ def test_solution_calculator_surface(
     from physicsnemo.models.domino.solutions import SolutionCalculatorSurface
     from physicsnemo.nn import FourierMLP, get_activation
 
-    torch.manual_seed(0)
-
     activation = get_activation("relu")
 
     # Determine input features based on surface configuration
@@ -199,8 +195,6 @@ def test_sample_sphere(device, r, num_points):
     """Test sphere sampling function"""
     from physicsnemo.models.domino.solutions import sample_sphere
 
-    torch.manual_seed(0)
-
     center = torch.randn(2, 30, 3).to(device)
     output = sample_sphere(center, r, num_points)
 
@@ -215,8 +209,6 @@ def test_sample_sphere_shell(device):
     """Test spherical shell sampling function"""
     from physicsnemo.models.domino.solutions import sample_sphere_shell
 
-    torch.manual_seed(0)
-
     center = torch.randn(2, 30, 3).to(device)
     r_inner, r_outer = 0.5, 1.5
     num_points = 50