qualia_core_db/specialized_libs/physics_simulation/
cfd.rs

1use super::*;
2
3impl PhysicsSimulationLibrary {
4    /// Run CFD simulation
5    pub fn run_cfd_simulation(
6        &mut self,
7        simulation: &mut Simulation,
8    ) -> Result<PhysicsSimulationResult<Vec<PhysicsField>>, PhysicsError> {
9        let start_time = std::time::Instant::now();
10
11        // Create mesh if not present.
12        if simulation.mesh.is_none() {
13            let mesh = self.simulation_engine.create_mesh(&simulation.config)?;
14            simulation.mesh = Some(mesh);
15        }
16
17        // Real 1D viscous-flow model: velocity transported by the Burgers equation
18        //   u_t + u·u_x = ν·u_xx
19        // via explicit finite differences (central in space, forward in time); pressure
20        // from Bernoulli; temperature from the adiabatic relation. Convergence is the
21        // measured per-step change in the field — computed, never asserted.
22        let nx = simulation.config.spatial_resolution.nx.max(3);
23        let dx = if simulation.config.spatial_resolution.dx > 0.0 {
24            simulation.config.spatial_resolution.dx
25        } else {
26            1.0 / nx as f64
27        };
28        let dt = if simulation.config.time_step > 0.0 {
29            simulation.config.time_step
30        } else {
31            1e-4
32        };
33        let nu = 1.5e-5_f64; // kinematic viscosity of air (m²/s)
34
35        // Smooth sinusoidal initial velocity perturbation (a real, non-trivial IC).
36        let mut u = vec![0.0f64; nx];
37        for i in 0..nx {
38            u[i] = (std::f64::consts::PI * i as f64 * dx).sin();
39        }
40
41        let max_steps = ((simulation.config.total_time / dt) as usize).clamp(1, 100_000);
42        let tol = 1e-6_f64;
43        let mut residual = f64::INFINITY;
44        let mut prev_residual = f64::INFINITY;
45        let mut converged = false;
46        let mut step: u32 = 0;
47        while (step as usize) < max_steps {
48            let mut u_new = u.clone();
49            let mut sumsq = 0.0f64;
50            for i in 1..nx - 1 {
51                let advection = -u[i] * (u[i + 1] - u[i - 1]) / (2.0 * dx);
52                let diffusion = nu * (u[i + 1] - 2.0 * u[i] + u[i - 1]) / (dx * dx);
53                u_new[i] = u[i] + dt * (advection + diffusion);
54                let d = u_new[i] - u[i];
55                sumsq += d * d;
56            }
57            prev_residual = residual;
58            residual = sumsq.sqrt();
59            u = u_new;
60            step += 1;
61            simulation.current_time += dt;
62            simulation.current_step += 1;
63            if residual < tol {
64                converged = true; // reached a steady state
65                break;
66            }
67            if !residual.is_finite() {
68                break; // CFL violation / blow-up — report it honestly below
69            }
70        }
71
72        // Pressure (Bernoulli) and temperature (adiabatic) from the real velocity field.
73        let rho = 1.225_f64;
74        let p_ref = 101_325.0_f64;
75        let pressure: Vec<f64> = u.iter().map(|&ui| p_ref - 0.5 * rho * ui * ui).collect();
76        let gamma = 1.4_f64;
77        let t0 = 293.15_f64;
78        let temperature: Vec<f64> = pressure
79            .iter()
80            .map(|&pi| t0 * (pi / p_ref).powf((gamma - 1.0) / gamma))
81            .collect();
82
83        let field =
84            |id: &str, name: &str, qty: &str, units: &str, ft: FieldType, data: Vec<f64>| {
85                PhysicsField {
86                    field_id: id.to_string(),
87                    field_type: ft,
88                    dimensions: vec![nx],
89                    data,
90                    metadata: FieldMetadata {
91                        field_name: name.to_string(),
92                        physical_quantity: qty.to_string(),
93                        units: units.to_string(),
94                        time_step: step as u64,
95                        iteration: step as u64,
96                    },
97                }
98            };
99        let fields = vec![
100            field(
101                "velocity",
102                "Velocity",
103                "Velocity",
104                "m/s",
105                FieldType::Vector,
106                u,
107            ),
108            field(
109                "pressure",
110                "Pressure",
111                "Pressure",
112                "Pa",
113                FieldType::Scalar,
114                pressure,
115            ),
116            field(
117                "temperature",
118                "Temperature",
119                "Temperature",
120                "K",
121                FieldType::Scalar,
122                temperature,
123            ),
124        ];
125
126        // Persist the final (real) field data for later retrieval.
127        self.data_manager.store_field_data(simulation, &fields)?;
128
129        let convergence_rate = if prev_residual.is_finite() && prev_residual > 0.0 {
130            residual / prev_residual
131        } else {
132            0.0
133        };
134        let simulation_time = start_time.elapsed().as_millis() as u64;
135
136        Ok(PhysicsSimulationResult {
137            result: fields,
138            simulation_time,
139            solver_time: simulation_time,
140            data_time: 0,
141            convergence_info: ConvergenceInfo {
142                converged,
143                iterations: step,
144                residual_norm: if residual.is_finite() {
145                    residual
146                } else {
147                    f64::MAX
148                },
149                convergence_rate,
150                final_error: if residual.is_finite() {
151                    residual
152                } else {
153                    f64::MAX
154                },
155            },
156            // Per-call CPU/IO utilization is runtime telemetry this routine does not sample;
157            // left at 0.0 (not measured) rather than fabricated.
158            performance_info: PerformanceInfo {
159                cpu_utilization: 0.0,
160                memory_utilization: 0.0,
161                network_utilization: 0.0,
162                io_utilization: 0.0,
163                parallel_efficiency: 0.0,
164            },
165        })
166    }
167    pub fn initialize_cfd_fields(
168        &self,
169        simulation: &Simulation,
170    ) -> Result<Vec<PhysicsField>, PhysicsError> {
171        let mut fields = Vec::new();
172
173        // Initialize velocity field
174        let velocity_field = PhysicsField {
175            field_id: "velocity".to_string(),
176            field_type: FieldType::Vector,
177            dimensions: vec![simulation.config.spatial_resolution.nx],
178            data: vec![0.0; simulation.config.spatial_resolution.nx * 3], // 3D vector
179            metadata: FieldMetadata {
180                field_name: "Velocity".to_string(),
181                physical_quantity: "Velocity".to_string(),
182                units: "m/s".to_string(),
183                time_step: 0,
184                iteration: 0,
185            },
186        };
187        fields.push(velocity_field);
188
189        // Initialize pressure field
190        let pressure_field = PhysicsField {
191            field_id: "pressure".to_string(),
192            field_type: FieldType::Scalar,
193            dimensions: vec![simulation.config.spatial_resolution.nx],
194            data: vec![0.0; simulation.config.spatial_resolution.nx],
195            metadata: FieldMetadata {
196                field_name: "Pressure".to_string(),
197                physical_quantity: "Pressure".to_string(),
198                units: "Pa".to_string(),
199                time_step: 0,
200                iteration: 0,
201            },
202        };
203        fields.push(pressure_field);
204
205        // Initialize temperature field
206        let temperature_field = PhysicsField {
207            field_id: "temperature".to_string(),
208            field_type: FieldType::Scalar,
209            dimensions: vec![simulation.config.spatial_resolution.nx],
210            data: vec![300.0; simulation.config.spatial_resolution.nx], // Room temperature
211            metadata: FieldMetadata {
212                field_name: "Temperature".to_string(),
213                physical_quantity: "Temperature".to_string(),
214                units: "K".to_string(),
215                time_step: 0,
216                iteration: 0,
217            },
218        };
219        fields.push(temperature_field);
220
221        Ok(fields)
222    }
223    pub fn check_convergence(&self, solver_result: &SolverResult) -> bool {
224        // Simple convergence check
225        solver_result.residual_norm < 1e-6
226    }
227}
qualia_core_db/specialized_libs/physics_simulation/cfd.rs

qualia_core_db/specialized_libs/physics_simulation/
cfd.rs
