Repository: NREL/floris

Issue:
Feature request: Enable Cumulative Curl model in `solve_for_points()`
Over in #1038, @reos-aburns pointed out that we have not yet enabled the Cumulative Curl model to be used with the `solve_for_points()` procedure. I believe that this should not be too difficult to enable, and I'm going to look into it to see what I can do.

Discussion:

Solution patch:
diff --git a/docs/advanced_concepts.ipynb b/docs/advanced_concepts.ipynb
index 45fb9c017..023ca3466 100644
--- a/docs/advanced_concepts.ipynb
+++ b/docs/advanced_concepts.ipynb
@@ -95,16 +95,148 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## FLORIS as a library\n",
+    "\n",
+    "FLORIS is commonly used as a library in other software packages.\n",
+    "In cases where the calling-code will create inputs for FLORIS rather than require them from the\n",
+    "user, it can be helpful to initialize the FLORIS model with default inputs and then\n",
+    "change them in code.\n",
+    "In this case, the following workflow is recommended."
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "import floris\n",
+    "\n",
+    "# Initialize FLORIS with defaults\n",
+    "fmodel = floris.FlorisModel(\"defaults\")\n",
+    "\n",
+    "# Within the calling-code's setup step, update FLORIS as needed\n",
+    "fmodel.set(\n",
+    "    wind_directions=[i for i in range(10)],\n",
+    "    wind_speeds=[5 + i for i in range(10)],\n",
+    "    turbulence_intensities=[i for i in range(10)],\n",
+    "    # turbine_library_path=\"path/to/turbine_library\",   # Shown here for reference\n",
+    "    # turbine_type=[\"my_turbine\"]\n",
+    ")\n",
+    "\n",
+    "# Within the calling code's computation, run FLORIS\n",
+    "fmodel.run()"
+   ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
+   "source": [
+    "Alternatively, the calling-code can import the FLORIS default inputs as a Python dictionary\n",
+    "and modify them directly before initializing the FLORIS model.\n",
+    "This is especially helpful when the calling-code will modify a parameter that isn't\n",
+    "supported by the `FlorisModel.set(...)` command.\n",
+    "In particular, the wake model parameters are not directly accessible, so these can be updated\n",
+    "externally, as shown below.\n",
+    "Note that the `FlorisModel.get_defaults()` function returns a deep copy of the default inputs,\n",
+    "so these can be modified directly without side effects."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "solver\n",
+      "    type\n",
+      "        turbine_grid\n",
+      "    turbine_grid_points\n",
+      "        3\n",
+      "wake\n",
+      "    model_strings\n",
+      "        combination_model\n",
+      "            sosfs\n",
+      "        deflection_model\n",
+      "            gauss\n",
+      "        turbulence_model\n",
+      "            crespo_hernandez\n",
+      "        velocity_model\n",
+      "            jensen\n",
+      "farm\n",
+      "    layout_x\n",
+      "        [0.0]\n",
+      "    layout_y\n",
+      "        [0.0]\n",
+      "    turbine_type\n",
+      "        ['nrel_5MW']\n",
+      "    turbine_library_path\n",
+      "        /Users/rmudafor/Development/floris/floris/turbine_library\n",
+      "flow_field\n",
+      "    wind_speeds\n",
+      "        [5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0]\n",
+      "    wind_directions\n",
+      "        [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]\n",
+      "    wind_veer\n",
+      "        0.0\n",
+      "    wind_shear\n",
+      "        0.12\n",
+      "    air_density\n",
+      "        1.225\n",
+      "    turbulence_intensities\n",
+      "        [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]\n",
+      "    reference_wind_height\n",
+      "        90.0\n",
+      "name\n",
+      "    GCH\n",
+      "description\n",
+      "    Default initialization: Gauss-Curl hybrid model (GCH)\n",
+      "floris_version\n",
+      "    v4\n"
+     ]
+    }
+   ],
+   "source": [
+    "import floris\n",
+    "\n",
+    "# Retrieve the default parameters\n",
+    "fdefaults = floris.FlorisModel.get_defaults()\n",
+    "\n",
+    "# Update wake model parameters\n",
+    "fdefaults[\"wake\"][\"model_strings\"][\"velocity_model\"] = \"jensen\"\n",
+    "fdefaults[\"wake\"][\"wake_velocity_parameters\"][\"jensen\"][\"we\"] = 0.05\n",
+    "\n",
+    "# Initialize FLORIS with modified parameters\n",
+    "fmodel = floris.FlorisModel(configuration=fdefaults)\n",
+    "\n",
+    "# Within the calling-code's setup step, update FLORIS as needed\n",
+    "fmodel.set(\n",
+    "    wind_directions=[i for i in range(10)],\n",
+    "    wind_speeds=[5 + i for i in range(10)],\n",
+    "    turbulence_intensities=[i for i in range(10)],\n",
+    "    # turbine_library_path=\"path/to/turbine_library\",   # Shown here for reference\n",
+    "    # turbine_type=[\"my_turbine\"]\n",
+    ")\n",
+    "\n",
+    "# Verify settings are correct\n",
+    "fmodel.show_config()  # Shows truncated set of inputs; show all with fmodel.show_config(full=True)\n",
+    "\n",
+    "# Within the calling code's computation, run FLORIS\n",
+    "fmodel.run()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": []
   }
  ],
diff --git a/examples/examples_get_flow/002_extract_wind_speed_at_points.py b/examples/examples_get_flow/002_extract_wind_speed_at_points.py
index aaf086f4b..663124140 100644
--- a/examples/examples_get_flow/002_extract_wind_speed_at_points.py
+++ b/examples/examples_get_flow/002_extract_wind_speed_at_points.py
@@ -20,8 +20,8 @@
 
 
 # User options
-# FLORIS model to use (limited to Gauss/GCH, Jensen, and empirical Gauss)
-floris_model = "gch"  # Try "gch", "jensen", "emgauss"
+# FLORIS model to use (legacy Turbopark not available)
+floris_model = "gch"  # Try "gch", "cc", "jensen", "emgauss", "turboparkgauss"
 # Option to try different met mast locations
 met_mast_option = 0  # Try 0, 1, 2, 3
 
diff --git a/floris/core/core.py b/floris/core/core.py
index d12005ba5..857b90fa4 100644
--- a/floris/core/core.py
+++ b/floris/core/core.py
@@ -252,13 +252,15 @@ def solve_for_points(self, x, y, z):
 
         vel_model = self.wake.model_strings["velocity_model"]
 
-        if vel_model == "cc" or vel_model == "turbopark":
+        if vel_model == "turbopark":
             raise NotImplementedError(
-                "solve_for_points is currently only available with the "+\
-                "gauss, jensen, and empirical_gauss models."
+                "solve_for_points is not available for the legacy \'turbopark\' model. "
+                "However, it is available for \'turboparkgauss\'."
             )
         elif vel_model == "empirical_gauss":
             full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake)
+        elif vel_model == "cc":
+            full_flow_cc_solver(self.farm, self.flow_field, field_grid, self.wake)
         else:
             full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake)
 
diff --git a/floris/core/solver.py b/floris/core/solver.py
index a7c3d8796..679a97c4a 100644
--- a/floris/core/solver.py
+++ b/floris/core/solver.py
@@ -694,7 +694,7 @@ def cc_solver(
     flow_field.turbulence_intensity_field_sorted_avg = np.mean(
         turbine_turbulence_intensity,
         axis=(2,3)
-    )
+    )[:, :, None, None]
 
 
 def full_flow_cc_solver(
@@ -774,7 +774,7 @@ def full_flow_cc_solver(
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
             power_setpoints=turbine_grid_farm.power_setpoints_sorted,
-            awc_modes=turbine_grid_farm.awc_modes,
+            awc_modes=turbine_grid_farm.awc_modes_sorted,
             awc_amplitudes=turbine_grid_farm.awc_amplitudes_sorted,
             thrust_coefficient_functions=turbine_grid_farm.turbine_thrust_coefficient_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
@@ -794,7 +794,7 @@ def full_flow_cc_solver(
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
             power_setpoints=turbine_grid_farm.power_setpoints_sorted,
-            awc_modes=turbine_grid_farm.awc_modes,
+            awc_modes=turbine_grid_farm.awc_modes_sorted,
             awc_amplitudes=turbine_grid_farm.awc_amplitudes_sorted,
             axial_induction_functions=turbine_grid_farm.turbine_axial_induction_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
@@ -937,7 +937,7 @@ def turbopark_solver(
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
             power_setpoints=farm.power_setpoints_sorted,
-            awc_modes=farm.awc_modes,
+            awc_modes=farm.awc_modes_sorted,
             awc_amplitudes=farm.awc_amplitudes_sorted,
             thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
@@ -956,7 +956,7 @@ def turbopark_solver(
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
             power_setpoints=farm.power_setpoints_sorted,
-            awc_modes=farm.awc_modes,
+            awc_modes=farm.awc_modes_sorted,
             awc_amplitudes=farm.awc_amplitudes_sorted,
             thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
@@ -978,7 +978,7 @@ def turbopark_solver(
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
             power_setpoints=farm.power_setpoints_sorted,
-            awc_modes=farm.awc_modes,
+            awc_modes=farm.awc_modes_sorted,
             awc_amplitudes=farm.awc_amplitudes_sorted,
             axial_induction_functions=farm.turbine_axial_induction_functions,
             tilt_interps=farm.turbine_tilt_interps,
@@ -1027,7 +1027,7 @@ def turbopark_solver(
                     yaw_angles=farm.yaw_angles_sorted,
                     tilt_angles=farm.tilt_angles_sorted,
                     power_setpoints=farm.power_setpoints_sorted,
-                    awc_modes=farm.awc_modes,
+                    awc_modes=farm.awc_modes_sorted,
                     awc_amplitudes=farm.awc_amplitudes_sorted,
                     thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
                     tilt_interps=farm.turbine_tilt_interps,
@@ -1122,7 +1122,7 @@ def turbopark_solver(
     flow_field.turbulence_intensity_field_sorted_avg = np.mean(
         turbine_turbulence_intensity,
         axis=(2, 3)
-    )
+    )[:, :, None, None]
 
 
 def full_flow_turbopark_solver(
diff --git a/floris/default_inputs.yaml b/floris/default_inputs.yaml
new file mode 100644
index 000000000..67e199559
--- /dev/null
+++ b/floris/default_inputs.yaml
@@ -0,0 +1,106 @@
+
+name: GCH
+description: "Default initialization: Gauss-Curl hybrid model (GCH)"
+floris_version: v4
+
+logging:
+  console:
+    enable: true
+    level: WARNING
+  file:
+    enable: false
+    level: WARNING
+
+solver:
+  type: turbine_grid
+  turbine_grid_points: 3
+
+farm:
+  layout_x:
+  - 0.0
+  layout_y:
+  - 0.0
+  turbine_type:
+  - nrel_5MW
+
+flow_field:
+  air_density: 1.225
+  reference_wind_height: -1
+  turbulence_intensities: []
+  wind_directions: []
+  wind_shear: 0.12
+  wind_speeds: []
+  wind_veer: 0.0
+
+wake:
+  model_strings:
+    combination_model: sosfs
+    deflection_model: gauss
+    turbulence_model: crespo_hernandez
+    velocity_model: gauss
+
+  enable_secondary_steering: true
+  enable_yaw_added_recovery: true
+  enable_transverse_velocities: true
+  enable_active_wake_mixing: false
+
+  wake_deflection_parameters:
+    gauss:
+      ad: 0.0
+      alpha: 0.58
+      bd: 0.0
+      beta: 0.077
+      dm: 1.0
+      ka: 0.38
+      kb: 0.004
+    jimenez:
+      ad: 0.0
+      bd: 0.0
+      kd: 0.05
+    empirical_gauss:
+      horizontal_deflection_gain_D: 3.0
+      vertical_deflection_gain_D: -1
+      deflection_rate: 22
+      mixing_gain_deflection: 0.0
+      yaw_added_mixing_gain: 0.0
+
+  wake_velocity_parameters:
+    gauss:
+      alpha: 0.58
+      beta: 0.077
+      ka: 0.38
+      kb: 0.004
+    jensen:
+      we: 0.05
+    cc:
+      a_s: 0.179367259
+      b_s: 0.0118889215
+      c_s1: 0.0563691592
+      c_s2: 0.13290157
+      a_f: 3.11
+      b_f: -0.68
+      c_f: 2.41
+      alpha_mod: 1.0
+    turbopark:
+      A: 0.04
+      sigma_max_rel: 4.0
+    turboparkgauss:
+      A: 0.04
+      include_mirror_wake: True
+    empirical_gauss:
+      wake_expansion_rates: [0.023, 0.008]
+      breakpoints_D: [10]
+      sigma_0_D: 0.28
+      smoothing_length_D: 2.0
+      mixing_gain_velocity: 2.0
+      awc_wake_exp: 1.2
+      awc_wake_denominator: 400
+
+  wake_turbulence_parameters:
+    crespo_hernandez:
+      initial: 0.1
+      constant: 0.5
+      ai: 0.8
+      downstream: -0.32
+    wake_induced_mixing:
+      atmospheric_ti_gain: 0.0
diff --git a/floris/floris_model.py b/floris/floris_model.py
index f521b5faf..edf987081 100644
--- a/floris/floris_model.py
+++ b/floris/floris_model.py
@@ -33,6 +33,7 @@
     NDArrayStr,
 )
 from floris.utilities import (
+    load_yaml,
     nested_get,
     nested_set,
     print_nested_dict,
@@ -63,7 +64,15 @@ class FlorisModel(LoggingManager):
                 - **logging**: See `floris.simulation.core.Core` for more details.
     """
 
+    @staticmethod
+    def get_defaults() -> dict:
+        return copy.deepcopy(load_yaml(Path(__file__).parent / "default_inputs.yaml"))
+
     def __init__(self, configuration: dict | str | Path):
+
+        if configuration == "defaults":
+            configuration = FlorisModel.get_defaults()
+
         self.configuration = configuration
 
         if isinstance(self.configuration, (str, Path)):
@@ -1627,11 +1636,29 @@ def get_turbine_layout(self, z=False):
         else:
             return xcoords, ycoords
 
-    def print_dict(self) -> None:
+    def show_config(self, full=False) -> None:
         """Print the FlorisModel dictionary.
         """
-        print_nested_dict(self.core.as_dict())
+        config_dict = self.core.as_dict()
+        if not full:
+            del config_dict["logging"]
+            del config_dict["wake"]["enable_secondary_steering"]
+            del config_dict["wake"]["enable_yaw_added_recovery"]
+            del config_dict["wake"]["enable_transverse_velocities"]
+            del config_dict["wake"]["enable_active_wake_mixing"]
+            del config_dict["wake"]["wake_deflection_parameters"]
+            del config_dict["wake"]["wake_velocity_parameters"]
+            del config_dict["wake"]["wake_turbulence_parameters"]
+        print_nested_dict(config_dict)
 
+    def print_dict(self) -> None:
+        """Print the FlorisModel dictionary.
+        """
+        self.logger.warning(
+            "The print_dict() method has been deprecated."
+            " Please use the show_config() method instead."
+        )
+        self.show_config(full=True)
 
     ### Properties
 
diff --git a/floris/par_floris_model.py b/floris/par_floris_model.py
index bfb60b300..8c3e5dce8 100644
--- a/floris/par_floris_model.py
+++ b/floris/par_floris_model.py
@@ -8,6 +8,9 @@
 
 from floris.core import State
 from floris.floris_model import FlorisModel
+from floris.type_dec import (
+    NDArrayFloat,
+)
 
 
 class ParFlorisModel(FlorisModel):
@@ -123,81 +126,111 @@ def run(self) -> None:
             t0 = timerpc()
             super().run()
             t1 = timerpc()
-        elif self.interface == "multiprocessing":
+            self._print_timings(t0, t1, None, None)
+        else:
             t0 = timerpc()
             self.core.initialize_domain()
             parallel_run_inputs = self._preprocessing()
             t1 = timerpc()
-            if self.return_turbine_powers_only:
-                with self._PoolExecutor(self.max_workers) as p:
-                    self._turbine_powers_split = p.starmap(
-                        _parallel_run_powers_only,
+            if self.interface == "multiprocessing":
+                if self.return_turbine_powers_only:
+                    with self._PoolExecutor(self.max_workers) as p:
+                        self._turbine_powers_split = p.starmap(
+                            _parallel_run_powers_only,
+                            parallel_run_inputs
+                        )
+                else:
+                    with self._PoolExecutor(self.max_workers) as p:
+                        self._fmodels_split = p.starmap(_parallel_run, parallel_run_inputs)
+            elif self.interface == "pathos":
+                if self.return_turbine_powers_only:
+                    self._turbine_powers_split = self.pathos_pool.map(
+                        _parallel_run_powers_only_map,
                         parallel_run_inputs
                     )
-            else:
-                with self._PoolExecutor(self.max_workers) as p:
-                    self._fmodels_split = p.starmap(_parallel_run, parallel_run_inputs)
+                else:
+                    self._fmodels_split = self.pathos_pool.map(
+                        _parallel_run_map,
+                        parallel_run_inputs
+                    )
+            elif self.interface == "concurrent":
+                if self.return_turbine_powers_only:
+                    with self._PoolExecutor(self.max_workers) as p:
+                        self._turbine_powers_split = p.map(
+                            _parallel_run_powers_only_map,
+                            parallel_run_inputs
+                        )
+                        self._turbine_powers_split = list(self._turbine_powers_split)
+                else:
+                    with self._PoolExecutor(self.max_workers) as p:
+                        self._fmodels_split = p.map(
+                            _parallel_run_map,
+                            parallel_run_inputs
+                        )
+                        self._fmodels_split = list(self._fmodels_split)
             t2 = timerpc()
             self._postprocessing()
             self.core.farm.finalize(self.core.grid.unsorted_indices)
             self.core.state = State.USED
             t3 = timerpc()
-        elif self.interface == "pathos":
+            self._print_timings(t0, t1, t2, t3)
+
+    def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloat):
+        """
+        Sample the flow field at specified points.
+
+        Args:
+            x: The x-coordinates of the points.
+            y: The y-coordinates of the points.
+            z: The z-coordinates of the points.
+
+        Returns:
+            NDArrayFloat: The wind speeds at the specified points.
+        """
+        if self.return_turbine_powers_only:
+            raise NotImplementedError(
+                "Sampling flow at points is not supported when "
+                "return_turbine_powers_only is set to True on ParFlorisModel."
+            )
+
+        if self.interface is None:
             t0 = timerpc()
-            self.core.initialize_domain()
-            parallel_run_inputs = self._preprocessing()
+            sampled_wind_speeds = super().sample_flow_at_points(x, y, z)
             t1 = timerpc()
-            if self.return_turbine_powers_only:
-                self._turbine_powers_split = self.pathos_pool.map(
-                    _parallel_run_powers_only_map,
-                    parallel_run_inputs
-                )
-            else:
-                self._fmodels_split = self.pathos_pool.map(
-                    _parallel_run_map,
-                    parallel_run_inputs
-                )
-            t2 = timerpc()
-            self._postprocessing()
-            self.core.farm.finalize(self.core.grid.unsorted_indices)
-            self.core.state = State.USED
-            t3 = timerpc()
-        elif self.interface == "concurrent":
+            self._print_timings(t0, t1, None, None)
+        else:
             t0 = timerpc()
             self.core.initialize_domain()
             parallel_run_inputs = self._preprocessing()
+            parallel_sample_flow_at_points_inputs = [
+                (fmodel_dict, control_setpoints, x, y, z)
+                for fmodel_dict, control_setpoints in parallel_run_inputs
+            ]
             t1 = timerpc()
-            if self.return_turbine_powers_only:
+            if self.interface == "multiprocessing":
                 with self._PoolExecutor(self.max_workers) as p:
-                    self._turbine_powers_split = p.map(
-                        _parallel_run_powers_only_map,
-                        parallel_run_inputs
+                    sampled_wind_speeds_p = p.starmap(
+                        _parallel_sample_flow_at_points,
+                        parallel_sample_flow_at_points_inputs
                     )
-                    self._turbine_powers_split = list(self._turbine_powers_split)
-            else:
+            elif self.interface == "pathos":
+                sampled_wind_speeds_p = self.pathos_pool.map(
+                    _parallel_sample_flow_at_points_map,
+                    parallel_sample_flow_at_points_inputs
+                )
+            elif self.interface == "concurrent":
                 with self._PoolExecutor(self.max_workers) as p:
-                    self._fmodels_split = p.map(
-                        _parallel_run_map,
-                        parallel_run_inputs
+                    sampled_wind_speeds_p = p.map(
+                        _parallel_sample_flow_at_points_map,
+                        parallel_sample_flow_at_points_inputs
                     )
-                    self._fmodels_split = list(self._fmodels_split)
+                    sampled_wind_speeds_p = list(sampled_wind_speeds_p)
             t2 = timerpc()
-            self._postprocessing()
-            self.core.farm.finalize(self.core.grid.unsorted_indices)
-            self.core.state = State.USED
+            sampled_wind_speeds = np.concatenate(sampled_wind_speeds_p, axis=0)
             t3 = timerpc()
-        if self.print_timings:
-            print("===============================================================================")
-            if self.interface is None:
-                print(f"Total time spent for serial calculation (interface=None): {t1 - t0:.3f} s")
-            else:
-                print(
-                    "Total time spent for parallel calculation "
-                    f"({self.max_workers} workers): {t3-t0:.3f} s"
-                )
-                print(f"  Time spent in parallel preprocessing: {t1-t0:.3f} s")
-                print(f"  Time spent in parallel loop execution: {t2-t1:.3f} s.")
-                print(f"  Time spent in parallel postprocessing: {t3-t2:.3f} s")
+            self._print_timings(t0, t1, t2, t3)
+
+        return sampled_wind_speeds
 
     def _preprocessing(self):
         """
@@ -278,6 +311,23 @@ def _postprocessing(self):
                     axis=0
                 )
 
+    def _print_timings(self, t0, t1, t2, t3):
+        """
+        Print the timings for the parallel execution.
+        """
+        if self.print_timings:
+            print("===============================================================================")
+            if self.interface is None:
+                print(f"Total time spent for serial calculation (interface=None): {t1 - t0:.3f} s")
+            else:
+                print(
+                    "Total time spent for parallel calculation "
+                    f"({self.max_workers} workers): {t3-t0:.3f} s"
+                )
+                print(f"  Time spent in parallel preprocessing: {t1-t0:.3f} s")
+                print(f"  Time spent in parallel loop execution: {t2-t1:.3f} s.")
+                print(f"  Time spent in parallel postprocessing: {t3-t2:.3f} s")
+
     def _get_turbine_powers(self):
         """
         Calculates the power at each turbine in the wind farm.
@@ -364,3 +414,14 @@ def _parallel_run_powers_only_map(x):
     Wrapper for unpacking inputs to _parallel_run_powers_only() for use with map().
     """
     return _parallel_run_powers_only(*x)
+
+def _parallel_sample_flow_at_points(fmodel_dict, set_kwargs, x, y, z):
+    fmodel = FlorisModel(fmodel_dict)
+    fmodel.set(**set_kwargs)
+    return fmodel.sample_flow_at_points(x, y, z)
+
+def _parallel_sample_flow_at_points_map(x):
+    """
+    Wrapper for unpacking inputs to _parallel_sample_flow_at_points() for use with map().
+    """
+    return _parallel_sample_flow_at_points(*x)
diff --git a/pyproject.toml b/pyproject.toml
index e64d5cc37..23fe7aead 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -9,8 +9,8 @@ description = "A controls-oriented engineering wake model."
 readme = "README.md"
 requires-python = ">=3.9"
 authors = [
-    { name = "Rafael Mudafort", email = "rafael.mudafort@nrel.gov" },
-    { name = "Paul Fleming", email = "paul.fleming@nrel.gov" },
+    { name = "Rafael Mudafort", email = "Rafael.Mudafort@nrel.gov" },
+    { name = "Paul Fleming", email = "Paul.Fleming@nrel.gov" },
     { name = "Michael (Misha) Sinner", email = "Michael.Sinner@nrel.gov" },
     { name = "Eric Simley", email = "Eric.Simley@nrel.gov" },
     { name = "Christopher Bay", email = "Christopher.Bay@nrel.gov" },
@@ -18,16 +18,16 @@ authors = [
 license = { file = "LICENSE.txt" }
 keywords = ["floris"]
 classifiers = [
-        "License :: OSI Approved :: BSD License",
-        "Programming Language :: Python",
-        "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.9",
-        "Programming Language :: Python :: 3.10",
-        "Programming Language :: Python :: 3.11",
-        "Programming Language :: Python :: 3.12",
-        "Programming Language :: Python :: 3.13",
-        "Programming Language :: Python :: Implementation :: CPython",
-        "Programming Language :: Python :: Implementation :: PyPy"
+    "License :: OSI Approved :: BSD License",
+    "Programming Language :: Python",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: Implementation :: CPython",
+    "Programming Language :: Python :: Implementation :: PyPy"
 ]
 dependencies = [
     "attrs",
@@ -63,7 +63,8 @@ include = ["floris*"]
 [tool.setuptools.package-data]
 floris = [
     "turbine_library/*.yaml",
-    "core/wake_velocity/turbopark_lookup_table.mat"
+    "core/wake_velocity/turbopark_lookup_table.mat",
+    "default_inputs.yaml"
 ]
 
 [project.urls]
@@ -76,7 +77,6 @@ Documentation = "https://nrel.github.io/floris/"
 branch = true
 source = "floris/*"
 omit = [
-    "setup.py",
     "tests/*"
 ]
 


Test patch:
diff --git a/tests/floris_model_integration_test.py b/tests/floris_model_integration_test.py
index eabe764db..3ba94b955 100644
--- a/tests/floris_model_integration_test.py
+++ b/tests/floris_model_integration_test.py
@@ -17,6 +17,24 @@
 YAML_INPUT = TEST_DATA / "input_full.yaml"
 
 
+def test_default_init():
+    # Test getting the default dict
+    defaults = FlorisModel.get_defaults()
+    fmodel1 = FlorisModel(defaults)
+    assert isinstance(fmodel1, FlorisModel)
+
+    # Test that there are no side effects from changing the default dict
+    defaults2 = FlorisModel.get_defaults()
+    defaults2["farm"]["layout_x"] = [0, 1000]
+    defaults2["farm"]["layout_y"] = [0, 0]
+    fmodel2 = FlorisModel(defaults2)
+    assert fmodel2.core.as_dict() != FlorisModel(defaults).core.as_dict()
+
+    # Test using the "default" string
+    # This checks that the init works and that the default dictionary hasn't changed
+    fmodel3 = FlorisModel("defaults")
+    assert fmodel1.core.as_dict() == fmodel3.core.as_dict()
+
 def test_read_yaml():
     fmodel = FlorisModel(configuration=YAML_INPUT)
     assert isinstance(fmodel, FlorisModel)
@@ -488,7 +506,6 @@ def test_expected_farm_value_regression():
     expected_farm_value = fmodel.get_expected_farm_value()
     assert np.allclose(expected_farm_value,75108001.05154414 , atol=1e-1)
 
-
 def test_get_farm_avp(caplog):
     fmodel = FlorisModel(configuration=YAML_INPUT)
 
diff --git a/tests/par_floris_model_unit_test.py b/tests/par_floris_model_unit_test.py
index 9e56e4d8c..9f5669a3c 100644
--- a/tests/par_floris_model_unit_test.py
+++ b/tests/par_floris_model_unit_test.py
@@ -286,3 +286,29 @@ def test_control_setpoints(sample_inputs_fixture):
 
     assert powers_fmodel.shape == powers_pfmodel.shape
     assert np.allclose(powers_fmodel, powers_pfmodel)
+
+def test_sample_flow_at_points(sample_inputs_fixture):
+
+    sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL
+    sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL
+
+    fmodel = FlorisModel(sample_inputs_fixture.core)
+
+    wind_speeds = np.array([8.0, 8.0, 9.0])
+    wind_directions = np.array([270.0, 270.0, 270.0])
+    fmodel.set(
+        wind_directions=wind_speeds.flatten(),
+        wind_speeds=wind_directions.flatten(),
+        turbulence_intensities=0.06 * np.ones_like(wind_speeds.flatten()),
+    )
+
+    x_test = np.array([500.0, 750.0, 1000.0, 1250.0, 1500.0])
+    y_test = np.array([0.0, 0.0, 0.0, 0.0, 0.0])
+    z_test = np.array([90.0, 90.0, 90.0, 90.0, 90.0])
+
+    ws_base = fmodel.sample_flow_at_points(x_test, y_test, z_test)
+
+    for interface in ["multiprocessing", "pathos", "concurrent"]:
+        pfmodel = ParFlorisModel(fmodel, max_workers=2, interface=interface)
+        ws_test = pfmodel.sample_flow_at_points(x_test, y_test, z_test)
+        assert np.allclose(ws_base, ws_test)