rework the flowcell logic and added a docs/example/flowcell_sampling.py

Author: MartinPdeS

FlowCyPy/cpp/fluidics/flow_cell/flow_cell.cpp

@@ -9,6 +9,7 @@ FlowCell::FlowCell(
     int N_terms,
     int n_int,
     const std::string& event_scheme,
+    const std::string& transverse_sampling_scheme,
     bool perfectly_aligned
 )
     : width(width),
@@ -20,34 +21,72 @@ FlowCell::FlowCell(
       N_terms(N_terms),
       n_int(n_int),
       event_scheme(event_scheme),
+      transverse_sampling_scheme(transverse_sampling_scheme),
       perfectly_aligned(perfectly_aligned)
 {
+    this->validate_physical_parameters();
     this->validate_event_scheme();
+    this->validate_transverse_sampling_scheme();
     this->initialize();
 }
 
 std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>
 FlowCell::sample_transverse_profile(int n_samples) const {
+    if (n_samples < 0) {
+        throw std::runtime_error("n_samples must be non negative.");
+    }
+
     std::vector<double> y_samples;
     std::vector<double> z_samples;
     std::vector<double> velocity_samples;
 
-    y_samples.reserve(n_samples);
-    z_samples.reserve(n_samples);
-    velocity_samples.reserve(n_samples);
+    y_samples.reserve(static_cast<std::size_t>(n_samples));
+    z_samples.reserve(static_cast<std::size_t>(n_samples));
+    velocity_samples.reserve(static_cast<std::size_t>(n_samples));
 
     std::random_device random_device;
     std::mt19937 generator(random_device());
+
     std::uniform_real_distribution<double> y_distribution(-sample.width / 2.0, sample.width / 2.0);
     std::uniform_real_distribution<double> z_distribution(-sample.height / 2.0, sample.height / 2.0);
+    std::uniform_real_distribution<double> acceptance_distribution(0.0, 1.0);
 
     for (int sample_index = 0; sample_index < n_samples; ++sample_index) {
-        const double y = perfectly_aligned ? 0.0 : y_distribution(generator);
-        const double z = perfectly_aligned ? 0.0 : z_distribution(generator);
+        double y = 0.0;
+        double z = 0.0;
+        double velocity = u_center;
+
+        if (!perfectly_aligned) {
+            bool accepted_sample = false;
+
+            while (!accepted_sample) {
+                y = y_distribution(generator);
+                z = z_distribution(generator);
+                velocity = this->get_velocity(y, z, dpdx);
+
+                if (transverse_sampling_scheme == "uniform-random") {
+                    accepted_sample = true;
+                } else if (transverse_sampling_scheme == "velocity-weighted") {
+                    const double acceptance_probability = velocity / u_center;
+
+                    if (
+                        acceptance_probability >= 1.0 ||
+                        acceptance_distribution(generator) <= acceptance_probability
+                    ) {
+                        accepted_sample = true;
+                    }
+                } else {
+                    throw std::runtime_error(
+                        "Unsupported transverse_sampling_scheme '" + transverse_sampling_scheme +
+                        "'. Expected one of: 'velocity-weighted', 'uniform-random'."
+                    );
+                }
+            }
+        }
 
         y_samples.push_back(y);
         z_samples.push_back(z);
-        velocity_samples.push_back(this->get_velocity(y, z, dpdx));
+        velocity_samples.push_back(velocity);
     }
 
     return std::make_tuple(y_samples, z_samples, velocity_samples);
@@ -70,20 +109,33 @@ void FlowCell::initialize() {
         throw std::runtime_error("Computed center velocity is non positive.");
     }
 
-    const double sample_area = sample_volume_flow / u_center;
-    const double sample_height = std::sqrt(sample_area * height / width);
-    const double sample_width = (width / height) * sample_height;
-    const double average_flow_speed_sample = sample_volume_flow / sample_area;
+    const double sample_core_scale = this->compute_sample_core_scale();
+    const double sample_width = sample_core_scale * width;
+    const double sample_height = sample_core_scale * height;
+    const double sample_area = sample_width * sample_height;
+    const double sample_average_flow_speed = sample_volume_flow / sample_area;
+    const double sample_max_flow_speed = this->get_velocity(0.0, 0.0, dpdx);
 
     this->sample = FluidRegion(
         sample_height,
         sample_width,
         sample_volume_flow,
-        u_center,
-        average_flow_speed_sample
+        sample_max_flow_speed,
+        sample_average_flow_speed
     );
 
-    this->sheath = FluidRegion(height, width, sheath_volume_flow);
+    const double sheath_area = area - sample_area;
+    const double sheath_average_flow_speed = sheath_area > 0.0
+        ? sheath_volume_flow / sheath_area
+        : 0.0;
+
+    this->sheath = FluidRegion(
+        height,
+        width,
+        sheath_volume_flow,
+        u_center,
+        sheath_average_flow_speed
+    );
 }
 
 void FlowCell::validate_event_scheme() const {
@@ -99,6 +151,99 @@ void FlowCell::validate_event_scheme() const {
     }
 }
 
+void FlowCell::validate_transverse_sampling_scheme() const {
+    if (
+        transverse_sampling_scheme != "velocity-weighted" &&
+        transverse_sampling_scheme != "uniform-random"
+    ) {
+        throw std::runtime_error(
+            "Unsupported transverse_sampling_scheme '" + transverse_sampling_scheme +
+            "'. Expected one of: 'velocity-weighted', 'uniform-random'."
+        );
+    }
+}
+
+void FlowCell::validate_physical_parameters() const {
+    if (width <= 0.0) {
+        throw std::runtime_error("width must be strictly positive.");
+    }
+
+    if (height <= 0.0) {
+        throw std::runtime_error("height must be strictly positive.");
+    }
+
+    if (sample_volume_flow < 0.0) {
+        throw std::runtime_error("sample_volume_flow must be non negative.");
+    }
+
+    if (sheath_volume_flow < 0.0) {
+        throw std::runtime_error("sheath_volume_flow must be non negative.");
+    }
+
+    if (sample_volume_flow + sheath_volume_flow <= 0.0) {
+        throw std::runtime_error("The total volume flow must be strictly positive.");
+    }
+
+    if (viscosity <= 0.0) {
+        throw std::runtime_error("viscosity must be strictly positive.");
+    }
+
+    if (N_terms <= 0) {
+        throw std::runtime_error("N_terms must be strictly positive.");
+    }
+
+    if (n_int < 2) {
+        throw std::runtime_error("n_int must be greater than or equal to 2.");
+    }
+}
+
+double FlowCell::compute_sample_core_scale() const {
+    if (sample_volume_flow == 0.0) {
+        return 0.0;
+    }
+
+    if (sheath_volume_flow == 0.0) {
+        return 1.0;
+    }
+
+    const double full_channel_flow = this->compute_region_flow(width, height, dpdx);
+
+    if (full_channel_flow <= 0.0) {
+        throw std::runtime_error("Computed full channel flow is non positive.");
+    }
+
+    if (sample_volume_flow > full_channel_flow) {
+        throw std::runtime_error("sample_volume_flow exceeds the computed full channel flow.");
+    }
+
+    double lower_scale = 0.0;
+    double upper_scale = 1.0;
+
+    constexpr int maximum_iterations = 80;
+    constexpr double relative_tolerance = 1e-10;
+
+    for (int iteration = 0; iteration < maximum_iterations; ++iteration) {
+        const double middle_scale = 0.5 * (lower_scale + upper_scale);
+        const double region_width = middle_scale * width;
+        const double region_height = middle_scale * height;
+        const double region_flow = this->compute_region_flow(region_width, region_height, dpdx);
+
+        const double relative_error = std::abs(region_flow - sample_volume_flow) / sample_volume_flow;
+
+        if (relative_error < relative_tolerance) {
+            return middle_scale;
+        }
+
+        if (region_flow < sample_volume_flow) {
+            lower_scale = middle_scale;
+        } else {
+            upper_scale = middle_scale;
+        }
+    }
+
+    return 0.5 * (lower_scale + upper_scale);
+}
+
 double FlowCell::get_velocity(double y, double z, double dpdx_local) const {
     double velocity = 0.0;
     const double prefactor =
@@ -117,11 +262,36 @@ double FlowCell::get_velocity(double y, double z, double dpdx_local) const {
     return prefactor * velocity;
 }
 
-double FlowCell::compute_channel_flow(double dpdx_input) const {
-    const double y_min = -width / 2.0;
-    const double y_max = width / 2.0;
-    const double z_min = -height / 2.0;
-    const double z_max = height / 2.0;
+double FlowCell::compute_region_flow(
+    const double region_width,
+    const double region_height,
+    const double dpdx_input
+) const
+{
+    if (region_width < 0.0) {
+        throw std::runtime_error("region_width must be non negative.");
+    }
+
+    if (region_height < 0.0) {
+        throw std::runtime_error("region_height must be non negative.");
+    }
+
+    if (region_width > width) {
+        throw std::runtime_error("region_width cannot exceed the channel width.");
+    }
+
+    if (region_height > height) {
+        throw std::runtime_error("region_height cannot exceed the channel height.");
+    }
+
+    if (region_width == 0.0 || region_height == 0.0) {
+        return 0.0;
+    }
+
+    const double y_min = -region_width / 2.0;
+    const double y_max = region_width / 2.0;
+    const double z_min = -region_height / 2.0;
+    const double z_max = region_height / 2.0;
 
     const double dy = (y_max - y_min) / static_cast<double>(n_int - 1);
     const double dz = (z_max - z_min) / static_cast<double>(n_int - 1);
@@ -140,6 +310,10 @@ double FlowCell::compute_channel_flow(double dpdx_input) const {
     return sum * dy * dz;
 }
 
+double FlowCell::compute_channel_flow(double dpdx_input) const {
+    return this->compute_region_flow(width, height, dpdx_input);
+}
+
 std::vector<double>
 FlowCell::sample_arrival_times_uniform_random(
     const std::size_t n_events,
@@ -250,4 +424,4 @@ FlowCell::sample_arrival_times(
         "Unsupported event_scheme '" + event_scheme +
         "'. Expected one of: 'uniform-random', 'linear', 'poisson'."
     );
-}
+}