diff --git a/HEU_MET_FAST_003/HEU_MET_FAST_003.py b/HEU_MET_FAST_003/HEU_MET_FAST_003.py
index 79efe11..974884c 100644
--- a/HEU_MET_FAST_003/HEU_MET_FAST_003.py
+++ b/HEU_MET_FAST_003/HEU_MET_FAST_003.py
@@ -104,7 +104,6 @@
             {"block_ids": [1], "xs": xs_oralloy},
             {"block_ids": [2], "xs": xs_tuballoy},
         ],
-        scattering_order=3,
         options={
             "use_precursors": False,
             "verbose_inner_iterations": True,
@@ -129,7 +128,7 @@
         print(f"Computed k-eigenvalue: {k}")
 
     # export
-    fflist = phys.GetScalarFieldFunctionList(only_scalar_flux=False)
+    fflist = phys.GetScalarFluxFieldFunction(only_scalar_flux=False)
     vtk_basename = "HEU_MET_FAST_003"
     # export only the flux of group g (first []), moment 0 (second [])
     FieldFunctionGridBased.ExportMultipleToPVTU(
diff --git a/OpenSn_Logo_CAD/opensn.py b/OpenSn_Logo_CAD/opensn.py
index 99e79da..8ee1118 100644
--- a/OpenSn_Logo_CAD/opensn.py
+++ b/OpenSn_Logo_CAD/opensn.py
@@ -59,7 +59,6 @@
             {"block_ids": [1], "xs": xs_source},
             {"block_ids": [2], "xs": xs_block},
         ],
-        scattering_order=0,
         volumetric_sources=[src],
     )
     ss_solver = SteadyStateSourceSolver(problem=phys)
@@ -67,5 +66,5 @@
     ss_solver.Execute()
 
     # Export results to VTK
-    fflist = phys.GetScalarFieldFunctionList()
+    fflist = phys.GetScalarFluxFieldFunction()
     FieldFunctionGridBased.ExportMultipleToPVTU(fflist, "flux")
diff --git a/Six_1g_spherical_benchmarks/Problem_1.py b/Six_1g_spherical_benchmarks/Problem_1.py
index 6b3ff26..2c54a67 100644
--- a/Six_1g_spherical_benchmarks/Problem_1.py
+++ b/Six_1g_spherical_benchmarks/Problem_1.py
@@ -102,7 +102,6 @@
         xs_map=[
             {"block_ids": [1], "xs": xs_mat},
         ],
-        scattering_order=0,
         boundary_conditions=[
             {"name": "sphere_surface", "type": "isotropic", "group_strength": bsrc},
         ],
@@ -113,7 +112,7 @@
     ss_solver.Execute()
 
     # Export
-    fflist = phys.GetFieldFunctions()
+    fflist = phys.GetScalarFluxFieldFunction()
 
     # Compute the average flux over a logical volume
     def average_flx_logvol(logvol):
diff --git a/Six_1g_spherical_benchmarks/Problem_2.py b/Six_1g_spherical_benchmarks/Problem_2.py
index da96c6e..a8b0688 100644
--- a/Six_1g_spherical_benchmarks/Problem_2.py
+++ b/Six_1g_spherical_benchmarks/Problem_2.py
@@ -66,7 +66,6 @@
         xs_map=[
             {"block_ids": [1, 2], "xs": xs_void},
         ],
-        scattering_order=0,
         volumetric_sources=[mg_src],
     )
 
@@ -75,7 +74,7 @@
     ss_solver.Execute()
 
     # export
-    fflist = phys.GetFieldFunctions()
+    fflist = phys.GetScalarFluxFieldFunction()
     """
         vtk_basename = "Problem_2_"
         # export only the scalar flux of group g
diff --git a/Six_1g_spherical_benchmarks/Problem_3.py b/Six_1g_spherical_benchmarks/Problem_3.py
index 313927d..42fa49f 100644
--- a/Six_1g_spherical_benchmarks/Problem_3.py
+++ b/Six_1g_spherical_benchmarks/Problem_3.py
@@ -63,7 +63,6 @@
         xs_map=[
             {"block_ids": [1, 2, 3], "xs": xs_void},
         ],
-        scattering_order=0,
         volumetric_sources=[mg_src],
     )
 
@@ -72,7 +71,7 @@
     ss_solver.Execute()
 
     # Export
-    fflist = phys.GetFieldFunctions()
+    fflist = phys.GetScalarFluxFieldFunction()
 
     # Compute the average flux over a logical volume
     def average_flx_logvol(logvol):
diff --git a/Six_1g_spherical_benchmarks/Problem_4.py b/Six_1g_spherical_benchmarks/Problem_4.py
index 8b1afb6..383fe18 100644
--- a/Six_1g_spherical_benchmarks/Problem_4.py
+++ b/Six_1g_spherical_benchmarks/Problem_4.py
@@ -66,7 +66,6 @@
             {"block_ids": [1], "xs": xs_mat},
             {"block_ids": [2], "xs": xs_void},
         ],
-        scattering_order=0,
         boundary_conditions=[
             {"name": "sphere_surface", "type": "isotropic", "group_strength": bsrc},
         ],
@@ -77,7 +76,7 @@
     ss_solver.Execute()
 
     # Export
-    fflist = phys.GetFieldFunctions()
+    fflist = phys.GetScalarFluxFieldFunction()
 
     # Compute the average flux over a logical volume
     def average_flx_logvol(logvol):
diff --git a/Six_1g_spherical_benchmarks/Problem_5.py b/Six_1g_spherical_benchmarks/Problem_5.py
index 1d5c18c..eb51d20 100644
--- a/Six_1g_spherical_benchmarks/Problem_5.py
+++ b/Six_1g_spherical_benchmarks/Problem_5.py
@@ -66,7 +66,6 @@
             {"block_ids": [2], "xs": xs_mat},
             {"block_ids": [1, 3], "xs": xs_void},
         ],
-        scattering_order=0,
         boundary_conditions=[
             {"name": "sphere_surface", "type": "isotropic", "group_strength": bsrc},
         ],
@@ -77,7 +76,7 @@
     ss_solver.Execute()
 
     # Export
-    fflist = phys.GetFieldFunctions()
+    fflist = phys.GetScalarFluxfieldFunction()
 
     # Compute the average flux over a logical volume
     def average_flx_logvol(logvol):
diff --git a/Six_1g_spherical_benchmarks/Problem_6.py b/Six_1g_spherical_benchmarks/Problem_6.py
index 3ef913d..dd16ca6 100644
--- a/Six_1g_spherical_benchmarks/Problem_6.py
+++ b/Six_1g_spherical_benchmarks/Problem_6.py
@@ -66,7 +66,6 @@
             {"block_ids": [1], "xs": xs_mat},
             {"block_ids": [2], "xs": xs_void},
         ],
-        scattering_order=0,
         volumetric_sources=[mg_src],
     )
 
@@ -75,7 +74,7 @@
     ss_solver.Execute()
 
     # Export
-    fflist = phys.GetFieldFunctions()
+    fflist = phys.GetScalarFluxFieldFunction()
 
     # Compute the average flux over a logical volume
     def average_flx_logvol(logvol):
diff --git a/Urban_Source/urban_source.py b/Urban_Source/urban_source.py
index 2e2afe7..ae0abd0 100644
--- a/Urban_Source/urban_source.py
+++ b/Urban_Source/urban_source.py
@@ -62,7 +62,6 @@
             {"block_ids": [2], "xs": xs_source},
             {"block_ids": [3], "xs": xs_air},
         ],
-        scattering_order=0,
         volumetric_sources=[src],
     )
     ss_solver = SteadyStateSourceSolver(problem=phys)
@@ -70,5 +69,5 @@
     ss_solver.Execute()
 
     # Export results to VTK
-    fflist = phys.GetScalarFieldFunctionList()
+    fflist = phys.GetScalarFluxFieldFunction()
     FieldFunctionGridBased.ExportMultipleToPVTU(fflist, "flux")