1use serde::{Deserialize, Serialize};
2use std::collections::{HashMap, VecDeque};
3use std::time::Instant;
4
5use crate::{
6 assembly::AssemblySummary,
7 diagnostics::{FeaDiagnostic, FeaDiagnosticSeverity},
8 physics::{
9 coupling::{electro_thermal, thermo_mechanical},
10 structural,
11 },
12 progress::{emit_phase, is_cancelled, FeaProgressPhase, FeaProgressStatus},
13 solve::runtime_tensor_solver::RuntimeTensorPreparedLinearSystem,
14 ComputeBackend, FeaElectroThermalContext, FeaPrepContext, FeaThermoMechanicalContext,
15};
16
17mod diagnostics;
18mod linear_step;
19
20use diagnostics::{push_transient_quality_diagnostics, TransientQualityDiagnosticInputs};
21use linear_step::{
22 build_step_rhs, solve_implicit_step_system, strain_energy, LinearStepStats,
23 RuntimeTensorStepCache,
24};
25
26#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
27pub struct TransientSolveOptions {
28 pub time_step_s: f64,
29 pub min_time_step_s: f64,
30 pub max_time_step_s: f64,
31 pub step_count: usize,
32 pub max_linear_iters: usize,
33 pub tolerance: f64,
34 pub residual_target: f64,
35 pub adaptive_time_step: bool,
36 pub max_step_retries: usize,
37 pub adapt_min_scale: f64,
38 pub adapt_max_scale: f64,
39 pub adapt_growth_exponent: f64,
40 pub adapt_retry_growth_cap: f64,
41 pub adapt_nonconverged_shrink: f64,
42 pub dt_bucket_rel_tolerance: f64,
43 #[serde(default = "default_progress_operation")]
44 pub progress_operation: String,
45 pub prep_context: Option<FeaPrepContext>,
46 pub thermo_mechanical_context: Option<FeaThermoMechanicalContext>,
47 pub electro_thermal_context: Option<FeaElectroThermalContext>,
48}
49
50fn default_progress_operation() -> String {
51 "fea.run_transient".to_string()
52}
53
54impl Default for TransientSolveOptions {
55 fn default() -> Self {
56 Self {
57 time_step_s: 1.0e-3,
58 min_time_step_s: 1.0e-6,
59 max_time_step_s: 2.0e-2,
60 step_count: 10,
61 max_linear_iters: 128,
62 tolerance: 1.0e-8,
63 residual_target: 1.0e-6,
64 adaptive_time_step: true,
65 max_step_retries: 4,
66 adapt_min_scale: 0.8,
67 adapt_max_scale: 1.25,
68 adapt_growth_exponent: 0.35,
69 adapt_retry_growth_cap: 1.05,
70 adapt_nonconverged_shrink: 0.75,
71 dt_bucket_rel_tolerance: 0.0,
72 progress_operation: default_progress_operation(),
73 prep_context: None,
74 thermo_mechanical_context: None,
75 electro_thermal_context: None,
76 }
77 }
78}
79
80#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
81pub struct TransientSolveResult {
82 pub converged_steps: usize,
83 pub total_steps: usize,
84 pub time_points_s: Vec<f64>,
85 pub displacement_snapshots: Vec<Vec<f64>>,
86 pub residual_norms: Vec<f64>,
87 pub accepted_time_steps_s: Vec<f64>,
88 pub diagnostics: Vec<FeaDiagnostic>,
89 pub solver_method: String,
90 pub solver_backend: String,
91 pub solver_host_sync_count: u32,
92 pub device_apply_k_count: u32,
93 pub device_apply_k_attempt_count: u32,
94 pub preconditioner: String,
95}
96
97pub fn solve_transient_system(
98 summary: &AssemblySummary,
99 options: TransientSolveOptions,
100 backend: ComputeBackend,
101) -> TransientSolveResult {
102 if summary.dof_count == 0 || options.step_count == 0 {
103 return TransientSolveResult {
104 converged_steps: 0,
105 total_steps: options.step_count,
106 time_points_s: vec![0.0],
107 displacement_snapshots: vec![vec![0.0; summary.dof_count]],
108 residual_norms: Vec::new(),
109 accepted_time_steps_s: Vec::new(),
110 diagnostics: vec![FeaDiagnostic {
111 code: "FEA_TRANSIENT_EMPTY_SYSTEM".to_string(),
112 severity: FeaDiagnosticSeverity::Warning,
113 message: "transient solve skipped because assembled system has zero DOFs or step_count is zero"
114 .to_string(),
115 }],
116 solver_method: "implicit_euler_pcg".to_string(),
117 solver_backend: "cpu_reference".to_string(),
118 solver_host_sync_count: 0,
119 device_apply_k_count: 0,
120 device_apply_k_attempt_count: 0,
121 preconditioner: "none".to_string(),
122 };
123 }
124
125 let use_runtime_tensor = backend == ComputeBackend::Gpu;
126 let thermo_context = options.thermo_mechanical_context.as_ref();
127 let electro_context = options.electro_thermal_context.as_ref();
128 let thermo_severity_base = thermo_mechanical::severity(thermo_context);
129 let thermo_temporal_variation = thermo_mechanical::temporal_profile_variation(thermo_context);
130 let electro_severity_base = electro_thermal::severity(electro_context);
131 let electro_temporal_variation = electro_thermal::temporal_profile_variation(electro_context);
132 let mut thermo_severity_sum = 0.0_f64;
133 let mut thermo_time_scale_sum = 0.0_f64;
134 let mut thermo_time_extrapolated = 0usize;
135 let mut thermo_time_clamped = 0usize;
136 let mut electro_severity_sum = 0.0_f64;
137 let mut electro_time_scale_sum = 0.0_f64;
138 let mut electro_severity_peak = 0.0_f64;
139 let mut thermo_severity_peak = 0.0_f64;
140 let mut effective_residual_target_peak = options.residual_target;
141 let mut thermo_growth_limit_min = 1.0_f64;
142 let mut thermo_nonconverged_shrink_min = 1.0_f64;
143 let min_dt = options.min_time_step_s.max(1.0e-9);
144 let max_dt = options.max_time_step_s.max(min_dt);
145 let mut dt = options.time_step_s.clamp(min_dt, max_dt);
146 let mut x = vec![0.0; summary.dof_count];
147 let mut time_points_s = vec![0.0];
148 let mut displacement_snapshots = vec![x.clone()];
149 let mut residual_norms = Vec::with_capacity(options.step_count);
150 let mut accepted_time_steps_s = Vec::with_capacity(options.step_count);
151 let mut converged_steps = 0usize;
152 let mut retry_budget_hits = 0usize;
153 let mut energies = Vec::with_capacity(options.step_count + 1);
154 energies.push(strain_energy(summary, &x));
155 let mut solver_backend = "cpu_reference".to_string();
156 let mut solver_host_sync_count = 0u32;
157 let mut device_apply_k_count = 0u32;
158 let mut device_apply_k_attempt_count = 0u32;
159 let mut selected_preconditioner = "none".to_string();
160 let mut prepared_runtime_systems_by_dt: HashMap<u64, RuntimeTensorPreparedLinearSystem> =
161 HashMap::new();
162 let mut prepared_runtime_system_lru = VecDeque::new();
163 let mut prepared_runtime_cache_hits = 0usize;
164 let mut prepared_runtime_cache_misses = 0usize;
165 let mut prepared_build_ms = 0.0_f64;
166 let mut solve_ms = 0.0_f64;
167 let mut fallback_apply_count = 0u32;
168 let mut adapt_increase_steps = 0usize;
169 let mut adapt_decrease_steps = 0usize;
170 let mut adapt_hold_steps = 0usize;
171 let mut adapt_scale_sum = 0.0_f64;
172 let mut adapt_scale_min = f64::INFINITY;
173 let mut adapt_scale_max = 0.0_f64;
174 let dt_bucket_rel_tolerance = options.dt_bucket_rel_tolerance.max(0.0);
175 let progress_operation = options.progress_operation.as_str();
176
177 for step_index in 0..options.step_count {
178 if is_cancelled() {
179 emit_phase(
180 progress_operation,
181 FeaProgressPhase::Solve,
182 FeaProgressStatus::Cancelled,
183 "transient solve cancelled",
184 Some(step_index as u64),
185 Some(options.step_count as u64),
186 );
187 break;
188 }
189 emit_phase(
190 progress_operation,
191 FeaProgressPhase::Solve,
192 FeaProgressStatus::Advanced,
193 format!("solving transient step {}", step_index + 1),
194 Some(step_index as u64),
195 Some(options.step_count as u64),
196 );
197 let step_progress = if options.step_count <= 1 {
198 1.0
199 } else {
200 step_index as f64 / (options.step_count - 1) as f64
201 };
202 let thermo_time_sample =
203 thermo_mechanical::sample_time_profile(thermo_context, step_progress);
204 let thermo_time_scale = thermo_time_sample.scale;
205 if thermo_time_sample.extrapolated {
206 thermo_time_extrapolated = thermo_time_extrapolated.saturating_add(1);
207 }
208 if thermo_time_sample.clamped {
209 thermo_time_clamped = thermo_time_clamped.saturating_add(1);
210 }
211 let electro_time_scale = electro_thermal::time_scale(electro_context, step_progress);
212 let thermo_severity = (thermo_severity_base * thermo_time_scale).clamp(0.0, 1.0);
213 let electro_severity = (electro_severity_base * electro_time_scale).clamp(0.0, 1.0);
214 thermo_severity_sum += thermo_severity;
215 thermo_time_scale_sum += thermo_time_scale;
216 thermo_severity_peak = thermo_severity_peak.max(thermo_severity);
217 electro_severity_sum += electro_severity;
218 electro_time_scale_sum += electro_time_scale;
219 electro_severity_peak = electro_severity_peak.max(electro_severity);
220 let thermo_policy =
221 thermo_mechanical::transient_policy(options.residual_target, thermo_severity);
222 let effective_residual_target = thermo_policy.effective_residual_target;
223 let thermo_growth_limit = thermo_policy.growth_limit;
224 let thermo_nonconverged_shrink = thermo_policy.nonconverged_shrink;
225 effective_residual_target_peak =
226 effective_residual_target_peak.max(effective_residual_target);
227 thermo_growth_limit_min = thermo_growth_limit_min.min(thermo_growth_limit);
228 thermo_nonconverged_shrink_min =
229 thermo_nonconverged_shrink_min.min(thermo_nonconverged_shrink);
230 let mut step_dt = dt;
231 let mut retries = 0usize;
232 let (next_x, residual_norm, converged, step_stats) = loop {
233 let rhs = build_step_rhs(summary, &x, step_dt);
234 let solve_start = Instant::now();
235 let solved = solve_implicit_step_system(
236 summary,
237 &rhs,
238 step_dt,
239 &options,
240 use_runtime_tensor,
241 RuntimeTensorStepCache {
242 prepared_systems_by_dt: &mut prepared_runtime_systems_by_dt,
243 prepared_lru: &mut prepared_runtime_system_lru,
244 cache_hits: &mut prepared_runtime_cache_hits,
245 cache_misses: &mut prepared_runtime_cache_misses,
246 prepared_build_ms: &mut prepared_build_ms,
247 dt_bucket_rel_tolerance,
248 },
249 );
250 solve_ms += solve_start.elapsed().as_secs_f64() * 1_000.0;
251 if !options.adaptive_time_step {
252 break solved;
253 }
254 let (candidate_x, candidate_residual, candidate_converged, candidate_stats) = solved;
255 if candidate_converged && candidate_residual <= effective_residual_target * 4.0 {
256 break (
257 candidate_x,
258 candidate_residual,
259 candidate_converged,
260 candidate_stats,
261 );
262 }
263 if retries >= options.max_step_retries || step_dt <= min_dt * 1.01 {
264 retry_budget_hits += 1;
265 break (
266 candidate_x,
267 candidate_residual,
268 candidate_converged,
269 candidate_stats,
270 );
271 }
272 step_dt = (step_dt * 0.5).clamp(min_dt, max_dt);
273 retries += 1;
274 };
275
276 if let Some(LinearStepStats {
277 solver_backend: step_solver_backend,
278 host_sync_count,
279 device_apply_k_count: step_device_apply_k_count,
280 device_apply_k_attempt_count: step_device_apply_k_attempt_count,
281 preconditioner,
282 }) = step_stats
283 {
284 solver_backend = step_solver_backend;
285 solver_host_sync_count = solver_host_sync_count.saturating_add(host_sync_count);
286 device_apply_k_count = device_apply_k_count.saturating_add(step_device_apply_k_count);
287 device_apply_k_attempt_count =
288 device_apply_k_attempt_count.saturating_add(step_device_apply_k_attempt_count);
289 fallback_apply_count = fallback_apply_count.saturating_add(
290 step_device_apply_k_attempt_count.saturating_sub(step_device_apply_k_count),
291 );
292 selected_preconditioner = preconditioner;
293 }
294
295 x = next_x;
296 let next_time = time_points_s.last().copied().unwrap_or(0.0) + step_dt;
297 time_points_s.push(next_time);
298 displacement_snapshots.push(x.clone());
299 residual_norms.push(residual_norm);
300 accepted_time_steps_s.push(step_dt);
301 energies.push(strain_energy(summary, &x));
302 if converged {
303 converged_steps += 1;
304 }
305
306 if options.adaptive_time_step {
307 let next_dt = recommend_next_time_step(structural::TransientAdaptivityInput {
308 step_dt,
309 residual_norm,
310 residual_target: effective_residual_target,
311 min_dt,
312 max_dt,
313 converged,
314 retries,
315 adapt_nonconverged_shrink: options.adapt_nonconverged_shrink,
316 adapt_growth_exponent: options.adapt_growth_exponent,
317 adapt_min_scale: options.adapt_min_scale,
318 adapt_max_scale: options.adapt_max_scale,
319 adapt_retry_growth_cap: options.adapt_retry_growth_cap,
320 thermo_growth_limit,
321 thermo_nonconverged_shrink,
322 });
323 let scale = next_dt / step_dt.max(1.0e-12);
324 adapt_scale_sum += scale;
325 adapt_scale_min = adapt_scale_min.min(scale);
326 adapt_scale_max = adapt_scale_max.max(scale);
327 if scale > 1.01 {
328 adapt_increase_steps += 1;
329 } else if scale < 0.99 {
330 adapt_decrease_steps += 1;
331 } else {
332 adapt_hold_steps += 1;
333 }
334 dt = next_dt;
335 } else {
336 dt = step_dt;
337 }
338 }
339
340 emit_phase(
341 progress_operation,
342 FeaProgressPhase::Solve,
343 FeaProgressStatus::Advanced,
344 "transient step solve loop complete",
345 Some(accepted_time_steps_s.len() as u64),
346 Some(options.step_count as u64),
347 );
348
349 let mut max_step_l2_jump_ratio = 0.0_f64;
350 let mut nonfinite_displacement_count = 0usize;
351 for window in displacement_snapshots.windows(2) {
352 let prev = &window[0];
353 let next = &window[1];
354 let prev_norm = prev.iter().map(|value| value * value).sum::<f64>().sqrt();
355 let next_norm = next.iter().map(|value| value * value).sum::<f64>().sqrt();
356 let mut jump_norm_sq = 0.0_f64;
357 for (a, b) in prev.iter().zip(next.iter()) {
358 let d = b - a;
359 jump_norm_sq += d * d;
360 if !b.is_finite() {
361 nonfinite_displacement_count += 1;
362 }
363 }
364 let jump_norm = jump_norm_sq.sqrt();
365 let jump_ratio = jump_norm / prev_norm.max(next_norm).max(1.0);
366 max_step_l2_jump_ratio = max_step_l2_jump_ratio.max(jump_ratio);
367 }
368
369 let mut diagnostics = vec![FeaDiagnostic {
370 code: "FEA_TRANSIENT_METHOD".to_string(),
371 severity: FeaDiagnosticSeverity::Info,
372 message: "solver=implicit_euler_pcg matrix_free=true".to_string(),
373 }];
374 push_transient_quality_diagnostics(
375 &mut diagnostics,
376 TransientQualityDiagnosticInputs {
377 options: &options,
378 dt_final: dt,
379 converged_steps,
380 retry_budget_hits,
381 accepted_time_steps_s: &accepted_time_steps_s,
382 residual_norms: &residual_norms,
383 energies: &energies,
384 use_runtime_tensor,
385 prepared_cache_entries: prepared_runtime_systems_by_dt.len(),
386 prepared_cache_hits: prepared_runtime_cache_hits,
387 prepared_cache_misses: prepared_runtime_cache_misses,
388 prepared_build_ms,
389 solve_ms,
390 fallback_apply_count,
391 adapt_increase_steps,
392 adapt_decrease_steps,
393 adapt_hold_steps,
394 adapt_scale_mean: if accepted_time_steps_s.is_empty() {
395 1.0
396 } else {
397 adapt_scale_sum / accepted_time_steps_s.len() as f64
398 },
399 adapt_scale_min: if adapt_scale_min.is_finite() {
400 adapt_scale_min
401 } else {
402 1.0
403 },
404 adapt_scale_max: if adapt_scale_max > 0.0 {
405 adapt_scale_max
406 } else {
407 1.0
408 },
409 dt_bucket_rel_tolerance,
410 max_step_l2_jump_ratio,
411 nonfinite_displacement_count,
412 thermo_severity_mean: if options.step_count == 0 {
413 0.0
414 } else {
415 thermo_severity_sum / options.step_count as f64
416 },
417 thermo_time_scale_mean: if options.step_count == 0 {
418 1.0
419 } else {
420 thermo_time_scale_sum / options.step_count as f64
421 },
422 thermo_severity_peak,
423 thermo_temporal_variation,
424 thermo_time_extrapolated,
425 thermo_time_clamped,
426 effective_residual_target_peak,
427 thermo_growth_limit_min,
428 thermo_nonconverged_shrink_min,
429 },
430 );
431 if electro_severity_peak > 0.0 {
432 diagnostics.push(FeaDiagnostic {
433 code: "FEA_ET_TRANSIENT".to_string(),
434 severity: if electro_severity_peak <= 0.6 && electro_temporal_variation <= 0.5 {
435 FeaDiagnosticSeverity::Info
436 } else {
437 FeaDiagnosticSeverity::Warning
438 },
439 message: format!(
440 "severity_mean={} time_scale_mean={} severity_peak={} temporal_variation={}",
441 if options.step_count == 0 {
442 0.0
443 } else {
444 electro_severity_sum / options.step_count as f64
445 },
446 if options.step_count == 0 {
447 1.0
448 } else {
449 electro_time_scale_sum / options.step_count as f64
450 },
451 electro_severity_peak,
452 electro_temporal_variation,
453 ),
454 });
455 }
456
457 TransientSolveResult {
458 converged_steps,
459 total_steps: options.step_count,
460 time_points_s,
461 displacement_snapshots,
462 residual_norms,
463 accepted_time_steps_s,
464 diagnostics,
465 solver_method: "implicit_euler_pcg".to_string(),
466 solver_backend,
467 solver_host_sync_count,
468 device_apply_k_count,
469 device_apply_k_attempt_count,
470 preconditioner: selected_preconditioner,
471 }
472}
473
474fn recommend_next_time_step(input: structural::TransientAdaptivityInput) -> f64 {
475 structural::recommend_next_time_step(input)
476}