Butterworth convolution filter for synthetic melt pool generation

examples/synthetic_api_example/rve.py

@@ -0,0 +1,106 @@
+# =============================================================================
+# Copyright (c) 2025 Oak Ridge National Laboratory
+#
+# All rights reserved.
+#
+# This file is part of Raptor.
+#
+# For details, see the top-level LICENSE file at:
+# https://github.com/ORNL-MDF/Raptor/LICENSE
+# =============================================================================
+import numpy as np
+from pathlib import Path
+from raptor.io import read_data
+from raptor.api import (
+    create_path_vectors,
+    create_melt_pool,
+    create_grid,
+    compute_porosity,
+    write_vtk,
+    compute_morphology,
+    write_morphology,
+    visualize,
+)
+from raptor.utilities import ScanPathBuilder, MeltPoolFilter
+
+# 1. Create voxel grid for the representative volume element (RVE)
+min_point = np.array([0.0, 0.0, 0.0])
+max_point = np.array([5.0e-4, 5.0e-4, 5.0e-4])
+bound_box = np.array([min_point, max_point])
+voxel_resolution = 5.0e-6
+
+grid = create_grid(voxel_resolution, bound_box=bound_box)
+
+# 2. Create path vectors through the representative volume element (RVE)
+power = 370
+velocity = 1.7
+hatch_spacing = 140e-6
+layer_height = 30e-6
+rotation = 67
+scan_extension = max(max_point - min_point)
+extra_layers = 0
+
+scan_path_builder = ScanPathBuilder(
+    bound_box,
+    power,
+    velocity,
+    hatch_spacing,
+    layer_height,
+    rotation,
+    scan_extension,
+    extra_layers,
+)
+
+scan_path_builder.generate_layers()
+path_vectors = scan_path_builder.process_vectors()
+
+# 3. Create melt pools from convolution filter
+mean = 148.0e-6
+std_dev = 18.0e-6
+frequency = 250000
+duration = 0.08
+
+# Instantiate object
+mp_filter = MeltPoolFilter(mean, std_dev, velocity, [frequency, duration])
+
+# Define physical scales
+mp_filter.add_effect("melt_pool", [800e-6, None, 1])
+
+# Generate stochastic melt pool
+mp_filter.initialize()
+width_data = mp_filter.generate_fluctuations(1)
+
+n_modes = 50
+
+# scale melt pool data by constant factor
+width_scale = 1.0
+depth_scale = 0.8
+height_scale = 0.4
+
+# assign shape to melt pool and cap (1 = parabola, 2 = ellipse)
+width_shape = 2  # placeholder
+height_shape = 1
+depth_shape = 1
+
+melt_pool_dict = {
+    "width": (width_data, n_modes, width_scale, width_shape),
+    "depth": (width_data, n_modes, depth_scale, depth_shape),
+    "height": (width_data, n_modes, height_scale, height_shape),
+}
+
+melt_pool = create_melt_pool(melt_pool_dict, enable_random_phases=True)
+
+# 4. Compute porosity using conic section / superellipse curves for melt pool mask
+porosity = compute_porosity(grid, path_vectors, melt_pool, jit_warmup=True)
+
+# 5. Write porosity field to .VTI
+write_vtk(grid.origin, grid.resolution, porosity, "rve.vti")
+
+# 6. Compute morphology
+morphology = compute_morphology(
+    porosity, voxel_resolution, ["area", "equivalent_diameter_area"]
+)
+write_morphology(morphology, "rve_morphology.csv")
+
+# 7. Visualize using PyVista
+visualize("./rve.vti")

src/raptor/utilities.py

@@ -11,6 +11,7 @@
 import numpy as np
 from typing import List, Tuple
 from .structures import PathVector
+from scipy.signal import lfilter, butter
 
 
 class ScanPathBuilder:
@@ -216,3 +217,107 @@ def write_layers(self, output_name, mode="layers"):
                 comments="",
             )
             print(f"Wrote file {filename}")
+
+
+class MeltPoolFilter:
+    def __init__(
+        self, mu: float, sigma: float, scan_speed: float, timeseries_params: list
+    ):
+        """
+        Filtration of disparate fluctuation scales to infer a melt pool oscillations sequence.
+        Uses scipy.butter to convolve scales of fluctuations together.
+        Initialization parameters:
+        mu: mean melt pool dimension
+        sigma: target standard deviation of the fluctuations
+        scan_speed: speed in m/s
+        timeseries_params: list of [fs,duration]
+        """
+        # statistical properties
+        self.mu, self.sigma = mu, sigma
+        # process parameters
+        self.scan_speed = scan_speed
+        # timeseries related properties
+        self.fs, self.duration = timeseries_params
+        self.dt = 1 / self.fs
+        self.n_points = int(self.duration / self.dt)
+        self.t = np.arange(0, self.duration, self.dt)
+        # parametric representations of fluctuation scales
+        self.physical_effects = {}  # contains scale description and parameters
+
+    def add_effect(self, effect_name: str, effect_params: list):
+        """
+        Adds a physical effect {effect_name} with parameters
+        length_scale_m,frequency_hz,sigma_weight = effect_params
+        to the MeltPoolFiltration.physical_effects dictionary.
+        """
+        length_scale_m, frequency_hz, sigma_weight = effect_params
+        self.physical_effects[effect_name] = {
+            "length_scale_m": length_scale_m,
+            "frequency_hz": frequency_hz,
+            "sigma_weight": sigma_weight,
+        }
+
+    def initialize(self):
+        # Calculate frequencies from length scales
+        for name, params in self.physical_effects.items():
+            if params["length_scale_m"] is not None:
+                params["frequency_hz"] = self.scan_speed / params["length_scale_m"]
+
+        # Check for Nyquist limit violations
+        max_freq = max(p["frequency_hz"] for p in self.physical_effects.values())
+        if max_freq > self.fs / 2:
+            raise ValueError(
+                f"Error: Maximum frequency ({max_freq/1000:.1f} kHz) exceeds Nyquist limit ({self.fs/2000:.1f} kHz). Increase sampling rate 'fs'."
+            )
+
+        # Normalize sigma weights so the variances sum correctly
+        weights = np.array([p["sigma_weight"] for p in self.physical_effects.values()])
+        sum_of_sq_weights = np.sum(weights**2)
+        self.normalization_factor = np.sqrt(sum_of_sq_weights)
+
+        for name, params in self.physical_effects.items():
+            params["sigma_contribution"] = (
+                params["sigma_weight"] / self.normalization_factor
+            ) * self.sigma
+
+    def bandpass_filter(self, data, f0, bandwidth_fraction, fs, order=4):
+        """Applies a bandpass filter around a center frequency f0."""
+        lowcut = f0 * (1 - bandwidth_fraction / 2)
+        highcut = f0 * (1 + bandwidth_fraction / 2)
+        nyq = 0.5 * fs
+        low = lowcut / nyq
+        high = highcut / nyq
+        if high >= 1.0:
+            high = 0.999
+        if low <= 0.0001:
+            low = 0.0001
+        b, a = butter(order, [low, high], btype="band")
+        return lfilter(b, a, data)
+
+    def generate_fluctuations(self, noise_scale):
+        base_white_noise = np.random.normal(
+            loc=0, scale=noise_scale, size=self.n_points
+        )
+        final_series = np.zeros(self.n_points)
+        self.component_series = {}
+
+        # Create each component series, scale it, and add to the final series
+        for name, params in self.physical_effects.items():
+            component_noise = self.bandpass_filter(
+                data=base_white_noise,
+                f0=params["frequency_hz"],
+                bandwidth_fraction=1,
+                fs=self.fs,
+            )
+
+            std_dev = np.std(component_noise)
+            scaled_component = component_noise * (
+                params["sigma_contribution"] / std_dev
+            )
+            self.component_series[name] = scaled_component
+            final_series += scaled_component
+
+        # Adding the mean
+        final_series += self.mu
+
+        return np.column_stack([self.t, final_series])