RustedSciThe 0.4.8

Rust framework for symbolic and numerical computing:BVP ( Newton-Raphson frozen/damped/with collocations ), IVP( BDF, Radau, Backward Euler, LSODE, LSODA, RK45, DoPri), nonlinear equations ( Levenberg, Gavin) and more
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
# BVP Damped/Frozen in RustedSciThe: user guide

This guide is written for a Rust developer who opens the BVP part of RustedSciThe for the first time and wants a practical mental model rather than a scavenger hunt through tests. We will cover what data the solver expects, how Damped and Frozen Newton differ, when Sparse or Banded is the right matrix route, what Numerical, Lambdify, and AOT actually mean, and what the current story/performance measurements tell us.

RustedSciThe's BVP solvers handle boundary value problems for systems of ordinary differential equations. You provide a first-order continuous system, a mesh, an initial guess, and boundary conditions. The solver discretizes the system, turns it into a large nonlinear algebraic problem, and solves that problem with a modified Newton method. In practice, configuring a BVP solve means configuring the whole path from continuous RHS to residual vector, Jacobian, linear solve, damping, and diagnostics.

## 1. Damped and Frozen: two Newton strategies

The main production path is Damped Newton, implemented in [`NR_Damp_solver_damped.rs`](NR_Damp_solver_damped.rs). It supports damping, bounds, Jacobian refreshes, diagnostics, generated backend handoff, and the broadest set of tested configurations. If you are not sure which strategy to use, start with Damped.

Frozen Newton, implemented in [`NR_Damp_solver_frozen.rs`](NR_Damp_solver_frozen.rs), is a specialized route with frozen linearization behavior. It is useful when you intentionally want to reuse or compare against a fixed-Jacobian style process. It is not the default "newer" or automatically faster route. A good workflow is: first make the problem converge and pass correctness checks with Damped, then use Frozen deliberately when its behavior matches the experiment.

Both solvers support symbolic problem setup through `Expr` and generated backend configuration. The pure numerical RHS-closure route is currently implemented for Damped.

## 2. What the BVP solver receives

The preferred constructor for new Damped code is `NRBVP::new_with_options(...)`:

```rust
pub fn new_with_options(
    eq_system: Vec<Expr>,
    initial_guess: DMatrix<f64>,
    values: Vec<String>,
    arg: String,
    border_conditions: HashMap<String, Vec<(usize, f64)>>,
    t0: f64,
    t_end: f64,
    n_steps: usize,
    options: DampedSolverOptions,
) -> NRBVP
```

`eq_system` is the first-order RHS system. For example, the scalar equation `y'' + y = 0` is usually rewritten as `y' = z`, `z' = -y`, so `eq_system = ["z", "-y"]` and `values = ["y", "z"]`. The order of `values` defines the state ordering and must match the RHS ordering.

`arg` is the independent variable name, usually `"x"`. The interval is given by `t0` and `t_end`, while `n_steps` is the number of mesh intervals.

`initial_guess` is the reduced unknown guess. This is an important detail in the current BVP Damped/Frozen API: the matrix is usually `values.len() x n_steps`, not `values.len() x (n_steps + 1)`. Boundary values are stored separately in `border_conditions` and are reconstructed into the full mesh state internally.

```rust
let initial_guess = DMatrix::from_vec(values.len(), n_steps, guess_values);
```

Boundary conditions use `HashMap<String, Vec<(usize, f64)>>`. The side flag `0` means the left boundary and `1` means the right boundary. For example, `y(0)=0` and `z(L)=1` are written as:

```rust
let border_conditions = HashMap::from([
    ("y".to_string(), vec![(0usize, 0.0)]),
    ("z".to_string(), vec![(1usize, 1.0)]),
]);
```

The Damped numeric route expects a well-posed first-order BVP: the total number of fixed endpoint conditions must equal `values.len()`. Those conditions do not have to be distributed one-per-variable. For example, the oscillator written as `y' = z`, `z' = -y` can be posed numerically with `y(0)=0` and `y(pi/2)=1`, leaving the derivative variable `z` free at both endpoints. Frozen does not use the numeric route.

## 3. `DampedSolverOptions`: the main configuration surface

`DampedSolverOptions` groups the values that used to be passed through long positional constructors. The explicit `new(...)` constructor is retained for compatibility, but new code should normally start from a preset and refine it with builder methods. In particular, do not pass derivative schemes as raw strings unless you are intentionally using a legacy/custom escape hatch.

In everyday code, start from one of the presets:

```rust
let options = DampedSolverOptions::sparse_damped();
let options = DampedSolverOptions::banded_damped();
let options = DampedSolverOptions::dense_damped();
```

Then refine tolerances, bounds, nonlinear strategy, and backend:

```rust
let strategy = SolverParams {
    max_jac: Some(6),
    max_damp_iter: Some(6),
    damp_factor: Some(0.5),
    adaptive: None,
};

let options = DampedSolverOptions::banded_damped()
    .trapezoid_derivative()
    .with_strategy_params(Some(strategy))
    .with_abs_tolerance(1e-8)
    .with_rel_tolerance(HashMap::from([
        ("y".to_string(), 1e-8),
        ("z".to_string(), 1e-8),
    ]))
    .with_bounds(HashMap::from([
        ("y".to_string(), (-2.0, 2.0)),
        ("z".to_string(), (-2.0, 2.0)),
    ]))
    .with_max_iterations(40)
    .with_loglevel(Some("none".to_string()));
```

The derivative-scheme builder selects the discretization scheme. `strategy` is usually `"Damped"`. `linear_sys_method` can normally stay `None` when the matrix route is selected through `method` and generated backend configuration. `method` is the broad matrix route: `"Dense"`, `"Sparse"`, or `"Banded"`.

`FrozenSolverOptions` is similar but smaller. It does not expose the same bounds and per-variable relative tolerance surface as Damped:

```rust
let options = FrozenSolverOptions::banded_frozen()
    .with_tolerance(1e-8)
    .with_max_iterations(40);
```

The derivative scheme builder is shared conceptually by Damped and Frozen:

```rust
use RustedSciThe::numerical::BVP_Damp::BvpDerivativeScheme;

let options = DampedSolverOptions::sparse_damped()
    .with_scheme(BvpDerivativeScheme::Trapezoid);

let options = DampedSolverOptions::sparse_damped()
    .forward_derivative();

let options = FrozenSolverOptions::banded_frozen()
    .trapezoid_derivative();
```

`with_scheme_name("...")` exists only as a compatibility escape hatch for legacy or experimental scheme names.

## 4. Dense, Sparse, Banded: a cost model, not a style preference

`Dense` is the small/debug path. It is useful for tiny systems and tests where simplicity matters more than memory or asymptotic cost.

`Sparse` is the general production route for large systems with sparse but not necessarily narrow Jacobians. It is a safe default for chemistry-like systems, irregular coupling, and cases where you do not want to assert banded structure too early.

`Banded` is the right route when nonzeros are concentrated near the diagonal. This is common after finite-difference discretization of first-order BVP/PDE-like systems: each node is coupled mostly to its neighbors. RustedSciThe's Banded route uses a LAPACK-style banded LU solver, so both storage and factorization exploit the bandwidth. On current combustion BVP story tests, Banded consistently reduces linear-solve time compared with Sparse when the problem really is narrow-band.

The warning is simple: Banded is excellent when the structure is truly banded. If the problem has long-range or irregular coupling, forcing Banded can hurt both robustness and performance.

## 5. Three backend routes: Numerical, Lambdify, AOT

For BVP Damped/Frozen, it helps to separate the Newton method from the evaluator backend. The nonlinear algorithm can remain the same while residual and Jacobian evaluation take different routes.

### Numerical path

The Damped Numerical path bypasses symbolic assembly. You provide a Rust RHS closure `f(x, y, params)`, and [`numeric_discretization.rs`](numeric_discretization.rs) builds the discrete residual directly from that closure. This route is meant for models whose source of truth is already Rust code, not `Expr`.

There are now two explicit Damped numerical entry points. `NRBVP::new_numeric_fd_with_options(...)` accepts the RHS closure and asks the solver to approximate the Newton Jacobian by finite differences. `NRBVP::new_numeric_with_jacobian_options(...)` accepts both the RHS closure and the small continuous Jacobian closure `df/dy`; the large discretized BVP Jacobian is assembled internally from that local Jacobian, the mesh, and the boundary-condition layout. This is the important ergonomic difference from the older low-level pattern: users no longer need to pass an empty `eq_system` just to say "this is a pure numeric problem".

Boundary handling follows the same first-order BVP rule as the symbolic route: provide exactly as many fixed endpoint values as there are state variables. They may be attached to different variables or to the same variable at two endpoints, as long as the resulting reduced system is square. A common two-point setup is:

```rust
let border_conditions = HashMap::from([(
    "y".to_string(),
    vec![(0usize, 0.0), (1usize, 1.0)],
)]);
```

Here `z = y'` is still part of the state vector and still appears in the RHS and Jacobian closures, but it has no direct endpoint value.

For the common sparse case there are shorter wrappers, `NRBVP::new_numeric_fd(...)` and `NRBVP::new_numeric_with_jacobian(...)`. They still ask for `bounds` and per-variable `rel_tolerance`, because those are not optional decoration in the damped Newton algorithm: bounds limit the admissible step and tolerances define the weighted convergence test. If you need Banded, custom nonlinear strategy parameters, or custom logging, use the `*_with_options` constructors.

Frozen is different. `BackendSelectionPolicy::NumericOnly` is intentionally rejected for the Frozen solver rather than silently falling back to Lambdify or finite differences. Frozen currently expects a symbolic Lambdify/AOT route with a prepared Jacobian callback.

### Lambdify

Lambdify means symbolic equations (`Expr`) are assembled into the BVP residual/Jacobian and lowered into executable Rust callbacks without an external compiler. It is the best correctness baseline: low friction, no dynamic artifact build, no `gcc`, `tcc`, or `zig` dependency.

Lambdify remains the best first run for a new formulation: at that stage the priority is checking equations, boundary conditions, and the initial guess without depending on an external toolchain. For small one-off solves it often wins end-to-end. After the `AtomView` optimizations, however, this should not be generalized to large problems: in the measured `n_steps = 1000` combustion BVP, cold `tcc` AOT beat Lambdify on both Sparse and Banded routes.

### AOT

AOT means ahead-of-time generated backend. Symbolic expressions are lowered through an intermediate representation, code is generated, a separate artifact is compiled, and the solver links that artifact as runtime callbacks. The philosophy is to pay preparation cost once and reduce per-call cost later.

RustedSciThe supports several AOT toolchains: C through `gcc` or `tcc`, Zig, and a Rust backend in the codegen layer. The corresponding compiler must be installed and visible in `PATH`. In current Windows release story runs, `tcc` is more than a cheap artifact compiler: on large `AtomView` problems it is competitive with Lambdify in full time-to-solution and, in several cold measurements, already wins. `gcc` can be attractive for runtime throughput but usually pays a higher startup price, while Zig sometimes has high or unstable build cost. Treat this as measured guidance, not a universal law; recheck ranking for another machine or equation structure.

AOT is unquestionably useful when you solve the same large symbolic problem repeatedly: parameter sweeps, continuation, batch runs, or production workflows where artifacts are reused. The boundary is friendlier than expected, though: with `AtomView` and `tcc`, AOT can be a reasonable choice even for one large solve, especially on the Banded route. For a small BVP, artifact preparation will still usually not pay for itself.

## 6. ExprLegacy and AtomView

The symbolic BVP pipeline has two assembly backends.

`ExprLegacy` is the older, established expression route. It remains valuable as a compatibility baseline.

`AtomView` is the newer route designed for efficient assembly and code generation. Separate process-isolated release measurements on the combustion BVP family show that it removes a particularly expensive symbolic-Jacobian construction cost in both Sparse and Banded routes while preserving the solution to the expected numerical tolerance. For that reason, the Sparse and Banded production presets for both Damped and Frozen solvers now select `AtomView` by default.

For a banded Lambdify solve the production spelling is now intentionally short:

```rust
use RustedSciThe::symbolic::symbolic_functions_BVP::BvpSymbolicAssemblyBackend;

let options = DampedSolverOptions::banded_damped()
    .with_banded_lambdify();

// Compatibility/control route when reproducing an ExprLegacy run:
let legacy_options = DampedSolverOptions::banded_damped()
    .with_banded_lambdify()
    .with_symbolic_assembly_backend(BvpSymbolicAssemblyBackend::ExprLegacy);
```

The same default is inherited by `FrozenSolverOptions::{sparse_frozen, banded_frozen}()` and by the Sparse/Banded AOT presets. `ExprLegacy` remains a supported explicit override for validation, compatibility runs, and investigation of frontend-sensitive behavior.

## 7. Generated backend config: Lambdify, AOT, and build policy

`GeneratedBackendConfig` is the central generated-backend configuration object. It stores backend selection policy, AOT build policy, codegen backend, C compiler, symbolic assembly backend, matrix backend override, and chunking policy.

High-level modes:

```rust
SparseGeneratedBackendMode::Defaults
SparseGeneratedBackendMode::RequirePrebuilt
SparseGeneratedBackendMode::BuildIfMissingRelease

BandedGeneratedBackendMode::Defaults
BandedGeneratedBackendMode::Lambdify
BandedGeneratedBackendMode::BuildIfMissingRelease
```

`Defaults` prefers compiled AOT when an artifact is already available and otherwise falls back to Lambdify. For both Sparse and Banded production modes it now also means `AtomView` symbolic assembly; Banded additionally selects the faithful LAPACK-style banded linear solver. `RequirePrebuilt` is useful in production because a missing artifact becomes a clear typed error instead of an unexpected build. `BuildIfMissingRelease` is convenient during development and experiments: the first run builds the artifact if needed.

Repeated-solve presets:

```rust
let sparse_aot = GeneratedBackendConfig::sparse_atomview_for_repeated_solves();
let banded_aot = GeneratedBackendConfig::banded_atomview_for_repeated_solves();
```

Explicit toolchain choices:

```rust
let cfg = GeneratedBackendConfig::banded_atomview_build_if_missing_release_gcc();
let cfg = GeneratedBackendConfig::banded_atomview_build_if_missing_release_tcc();
let cfg = GeneratedBackendConfig::banded_atomview_build_if_missing_release_zig();
```

`BuildIfMissing` exists so artifacts can be reused. The first run may pay for generation and compilation; later runs can use the compiled backend. In CI or deployment, use `RequirePrebuilt` when implicit compilation would be a bug.

A useful workflow for a new large symbolic problem has two passes. First run `AtomView + Lambdify` to establish a simple solution and matrix-structure control. Then compare it with `AtomView + tcc AOT` on the same `Sparse` or `Banded` route. In the recorded honest cold `combustion-1000` run, AOT `tcc/whole` completed in about `2.46 s` versus `2.84 s` for Banded Lambdify, and in about `2.68 s` versus `4.33 s` for Sparse Lambdify. This is strong enough evidence to test TCC AOT early, rather than reserving it only for late repeated-solve optimization.

There is now also a controlled warm result for the production lifecycle. In the
Damped Banded AtomView `combustion-1000` story, one setup solve built the `tcc`
artifact, and five later strict `RequirePrebuilt` solves were compared with five
Lambdify solves using alternating order and an equal five-second cooldown. Warm
AOT took `431.4 +/- 10.9 ms` versus `468.3 +/- 12.5 ms` for Lambdify, a modest
but repeatable improvement of about `7.9%`. The setup build cost means the
measured benefit pays back after roughly eight subsequent solves. This is the
intended meaning of the repeated-solve preset: not an unconditional one-shot
victory, but a compiled route that becomes useful when the same large model is
solved again.

### AOT artifact lifecycle and cleanup

AOT artifacts are generated under the crate `target` directory and are tracked
by the solver's `AotResolver` through a problem key and manifest metadata. Treat
that registry as the source of truth: do not hand-edit generated directories or
try to infer correctness from a filename alone. If the mesh/discretized problem,
symbolic frontend, matrix route, toolchain, or chunking policy changes, treat it
as a distinct artifact lifecycle and let the manifest/resolver contract decide
whether an existing artifact can be reused.

The recommended lifecycle is explicit:

1. Use `BuildIfMissing` for a setup/warmup run when the artifact is allowed to
   be created.
2. Switch to `RequirePrebuilt` for the measured or production loop, so a missing
   artifact is a typed error rather than hidden compilation.
3. If this is a cold-build benchmark or a cleanup/debug workflow, call
   `cleanup_registered_aot_artifacts()` only after the compiled callbacks from
   that solver are no longer needed.

```rust
use RustedSciThe::numerical::BVP_Damp::generated_solver_handoff::AotBuildPolicy;

// First pass: BuildIfMissing preset, for example
// GeneratedBackendConfig::banded_atomview_build_if_missing_release_tcc().
solver.try_eq_generate()?;
solver.solve();

// Strict repeated-solve pass: reuse the already built artifact.
let strict_cfg = solver
    .generated_backend_config()
    .clone()
    .with_aot_build_policy(AotBuildPolicy::RequirePrebuilt);
solver.set_generated_backend_config(strict_cfg);

solver.try_eq_generate()?;
solver.solve();

// Optional explicit cleanup for story/debug/cold-build workflows.
// Do this only when no live callback from this solver is still needed.
let removed = solver.cleanup_registered_aot_artifacts()?;
println!("removed {removed} generated AOT artifact trees");
```

Cleanup is deliberately not automatic on `drop`. Compiled callbacks may still be
alive in the process, and on Windows a loaded dynamic library can also keep file
handles open. The cleanup method only removes registered generated trees owned
by the current resolver snapshot; it does not unregister linked callbacks or
unload dynamic libraries. That conservative behavior is less magical, but it is
the safer production contract.

## 8. Chunking and parallel execution: the honest measurement story

AOT can split residual and sparse Jacobian value evaluation into chunks. The relevant user-level types are `AotChunkingPolicy`, `ResidualChunkingStrategy`, `SparseChunkingStrategy`, `AotExecutionPolicy`, and `ParallelExecutorConfig`.

```rust
use RustedSciThe::numerical::BVP_Damp::generated_solver_handoff::{
    AotChunkingPolicy, AotExecutionPolicy, GeneratedBackendConfig,
};
use RustedSciThe::symbolic::codegen::codegen_orchestrator::{
    ParallelExecutorConfig, ParallelFallbackPolicy,
};
use RustedSciThe::symbolic::codegen::codegen_runtime_api::ResidualChunkingStrategy;
use RustedSciThe::symbolic::codegen::codegen_tasks::SparseChunkingStrategy;

let chunking = AotChunkingPolicy::with_parts(
    Some(ResidualChunkingStrategy::ByTargetChunkCount { target_chunks: 4 }),
    Some(SparseChunkingStrategy::ByTargetChunkCount { target_chunks: 4 }),
);

let exec = ParallelExecutorConfig {
    jobs_per_worker: 1,
    max_residual_jobs: None,
    max_sparse_jobs: None,
    fallback_policy: ParallelFallbackPolicy::Auto,
};

let cfg = GeneratedBackendConfig::sparse_atomview_build_if_missing_release_tcc()
    .with_aot_chunking_policy(chunking)
    .with_aot_execution_policy(AotExecutionPolicy::Parallel(exec));
```

Current BVP combustion story tests show a nuanced result: chunking is correct and the true parallel path exists, but its value depends on scale and on the dominant stage. For medium `n_steps = 200` and `n_steps = 1000` problems, callback value evaluation can take only a few milliseconds, so parallel execution may barely move full wall clock. On the large Banded `AtomView` problem with `n_steps = 3000`, `tcc/chunk4` did execute four residual and four Jacobian jobs without fallback and reduced the hot callback intervals: residual values from about `13.8` to `7.1 ms`, and Jacobian values from about `6.7` to `1.7 ms`.

Even that large run does not make forced chunking a universal default: full cold times for `tcc/whole` and `tcc/chunk4` were statistically close, because hot callback work is only one part of the solve. This is why `Auto` fallback matters. If the chunk/job workload is too small, the runtime may intentionally use the sequential path to avoid scheduling overhead. That is not a hidden failure; it is the desired production behavior. If you want to study break-even points on your machine, use the diagnostic tests and inspect `actual_jobs`, `fallback`, `residual_ms`, `jacobian_ms`, and `linear_ms`. In production, leave chunking on `Auto` unless a measured model-specific configuration earns an explicit override.

The most direct production spelling is to leave chunking at its default and ask
the compiled backend to decide at runtime:

```rust
let generated = GeneratedBackendConfig::banded_atomview_build_if_missing_release_tcc()
    .with_aot_execution_policy(AotExecutionPolicy::Auto);
```

After `solver.try_solver()?`, inspect the diagnostics map rather than guessing
from the requested configuration:

```rust
let stats = solver.get_statistics();
let d = &stats.diagnostics;

let auto_mode = d.get("aot.auto.execution_mode");
let auto_residual_reason = d.get("aot.auto.residual.reason");
let auto_jacobian_reason = d.get("aot.auto.sparse_jacobian.reason");
let residual_jobs = d.get("aot.runtime.residual.actual_jobs");
let jacobian_jobs = d.get("aot.runtime.sparse_jacobian.actual_jobs");
let residual_fallback = d.get("aot.runtime.residual.fallback_reason");
let jacobian_fallback = d.get("aot.runtime.sparse_jacobian.fallback_reason");
let parallel_requested = d.get("aot.runtime.parallel_requested");

println!(
    "Auto={auto_mode:?}, residual={auto_residual_reason:?}/{residual_jobs:?}/{residual_fallback:?}, \
     jacobian={auto_jacobian_reason:?}/{jacobian_jobs:?}/{jacobian_fallback:?}, \
     parallel_requested={parallel_requested:?}"
);
```

The important point is semantic, not cosmetic: `aot.auto.*` describes the plan,
while `aot.runtime.*` describes what the linked callback actually did. In the
12-core `combustion-1000` AtomView+tcc story, `Auto` chose `Sequential` for both
Sparse and Banded. The plan said `work_per_chunk_too_small` for the Jacobian
stage and `single_chunk_or_job` for the residual stage; runtime diagnostics
matched that with `actual_jobs = 1`, `parallel_requested = false`, and
`single_requested_job`. That is a successful Auto decision. Forced chunking
should only be used when a larger model, or a break-even test, proves that there
is enough callback work per job to pay for scheduling.

### Comparing cold and warm runs

`total_ms` is meaningful only when the lifecycle included in a run is clear. A cold wall-clock test starts a fresh process and includes symbolic frontend construction, code generation, compilation, linking, and the Newton solve; it answers "how long from pressing the button until a result?" A warm or `RequirePrebuilt` run starts with an existing artifact and answers a different question: "how expensive is one more solve in a series?"

Both views matter. Cold measurements protect against an AOT pipeline that has excellent callbacks but takes too long before the first result. Warm or prebuilt measurements show its value in parameter sweeps, continuation, and batch workflows. Do not compare cold Lambdify with warm AOT, or draw conclusions from a table in which one row reused an artifact while another built it.

`RequirePrebuilt` is also a correctness contract, not merely a performance
switch. In the current Damped Sparse/Banded and Frozen Sparse/Banded lifecycle
stories, strict reuse remains on `AotCompiled`, reports no compiler/linker
interval, and preserves the Lambdify solution to roundoff. A prebuilt row that
compiles again or silently falls back to Lambdify should be treated as a
regression. In particular, Frozen combustion-1000 now has release-confirmed
artifact reuse on both matrix routes; the Sparse reuse rows agree with
Lambdify to `4.44e-16`.

## 9. End-to-end example: Symbolic Lambdify + Banded

This example solves `y' = z`, `z' = -y` with `y(0)=0`, `z(0)=1`. The exact solution is `y=sin(x)`, `z=cos(x)`, which makes it a good sanity check.

```rust
use std::collections::HashMap;
use nalgebra::{DMatrix, DVector};
use RustedSciThe::numerical::BVP_Damp::NR_Damp_solver_damped::{
    DampedSolverOptions, NRBVP,
};
use RustedSciThe::symbolic::symbolic_engine::Expr;

let n_steps = 80;
let values = vec!["y".to_string(), "z".to_string()];
let eq_system = vec![
    Expr::parse_expression("z"),
    Expr::parse_expression("-y"),
];

let t0 = 0.0;
let t_end = std::f64::consts::FRAC_PI_2;
let h = (t_end - t0) / n_steps as f64;
let mut guess = Vec::with_capacity(values.len() * n_steps);
for i in 0..n_steps {
    let x = t0 + i as f64 * h;
    guess.push(x.sin());
    guess.push(x.cos());
}
let initial_guess = DMatrix::from_column_slice(
    values.len(),
    n_steps,
    DVector::from_vec(guess).as_slice(),
);

let border_conditions = HashMap::from([
    ("y".to_string(), vec![(0usize, 0.0)]),
    ("z".to_string(), vec![(0usize, 1.0)]),
]);

let bounds = HashMap::from([
    ("y".to_string(), (-2.0, 2.0)),
    ("z".to_string(), (-2.0, 2.0)),
]);

let rel_tol = HashMap::from([
    ("y".to_string(), 1e-8),
    ("z".to_string(), 1e-8),
]);

let options = DampedSolverOptions::banded_damped()
    .with_banded_lambdify()
    .with_abs_tolerance(1e-8)
    .with_rel_tolerance(rel_tol)
    .with_bounds(bounds)
    .with_max_iterations(40)
    .with_loglevel(Some("none".to_string()));

let mut solver = NRBVP::new_with_options(
    eq_system,
    initial_guess,
    values,
    "x".to_string(),
    border_conditions,
    t0,
    t_end,
    n_steps,
    options,
);

solver.dont_save_log(true);
solver.try_solve()?;
let result = solver.get_result().expect("BVP result should be available");
let stats = solver.get_statistics();

println!("solution shape = {} x {}", result.nrows(), result.ncols());
println!("timers = {:?}", stats.timers);
# Ok::<(), Box<dyn std::error::Error>>(())
```

This is a strong first example because it uses symbolic equations but does not require an external compiler.

## 10. End-to-end example: Numerical path

The Numerical path is useful when the RHS closure is the source of truth. The constructor below does not receive symbolic equations at all. The discrete residual is built from the Rust closure, and the Jacobian choice is explicit: either finite differences through `new_numeric_fd_with_options`, or a user-provided continuous Jacobian through `new_numeric_with_jacobian_options`.

```rust
use std::collections::HashMap;
use nalgebra::{DMatrix, DVector};
use RustedSciThe::numerical::BVP_Damp::NR_Damp_solver_damped::{
    DampedSolverOptions, NRBVP,
};

let n_steps = 80;
let values = vec!["y".to_string(), "z".to_string()];
let t0 = 0.0;
let t_end = std::f64::consts::FRAC_PI_2;
let h = (t_end - t0) / n_steps as f64;

let mut guess = Vec::with_capacity(values.len() * n_steps);
for i in 0..n_steps {
    let x = t0 + i as f64 * h;
    guess.push(x.sin());
    guess.push(x.cos());
}
let initial_guess = DMatrix::from_column_slice(
    values.len(),
    n_steps,
    DVector::from_vec(guess).as_slice(),
);

let border_conditions = HashMap::from([
    ("y".to_string(), vec![(0usize, 0.0)]),
    ("z".to_string(), vec![(0usize, 1.0)]),
]);

let bounds = HashMap::from([
    ("y".to_string(), (-2.0, 2.0)),
    ("z".to_string(), (-2.0, 2.0)),
]);

let rel_tol = HashMap::from([
    ("y".to_string(), 1e-8),
    ("z".to_string(), 1e-8),
]);

let options = DampedSolverOptions::sparse_damped()
    .with_abs_tolerance(1e-8)
    .with_rel_tolerance(rel_tol)
    .with_bounds(bounds)
    .with_loglevel(Some("none".to_string()));

let mut solver = NRBVP::new_numeric_fd_with_options(
    initial_guess,
    values,
    "x".to_string(),
    border_conditions,
    t0,
    t_end,
    n_steps,
    options,
    |_x, y, _params| {
    DVector::from_vec(vec![y[1], -y[0]])
    },
);

solver.dont_save_log(true);
solver.try_solve()?;
# Ok::<(), Box<dyn std::error::Error>>(())
```

If you have an analytical local Jacobian, use the second numerical constructor. The closure returns the `2 x 2` continuous Jacobian of `f(x, y)` at one mesh node; the solver handles the global BVP matrix.

```rust
let mut solver = NRBVP::new_numeric_with_jacobian_options(
    initial_guess,
    values,
    "x".to_string(),
    border_conditions,
    t0,
    t_end,
    n_steps,
    options,
    |_x, y, _params| DVector::from_vec(vec![y[1], -y[0]]),
    |_x, _y, _params| DMatrix::from_row_slice(2, 2, &[0.0, 1.0, -1.0, 0.0]),
);
```

This route is intentionally free of symbolic semantics. If it fails, the likely source is numeric discretization, Newton choreography, or linear algebra rather than expression parsing or codegen.

## 11. End-to-end example: AtomView + AOT + Banded

For a large structured BVP that will be solved repeatedly, the typical production candidate is Banded + AtomView + AOT. The first run may be expensive because it builds the artifact, but repeated runs reuse compiled callbacks and the LAPACK-style banded linear solver.

```rust
use RustedSciThe::numerical::BVP_Damp::generated_solver_handoff::GeneratedBackendConfig;
use RustedSciThe::numerical::BVP_Damp::NR_Damp_solver_damped::DampedSolverOptions;

let generated = GeneratedBackendConfig::banded_atomview_for_repeated_solves();

let options = DampedSolverOptions::banded_damped()
    .with_generated_backend_config(generated)
    .with_abs_tolerance(1e-8)
    .with_max_iterations(60)
    .with_loglevel(Some("none".to_string()));
```

Explicit toolchain selection:

```rust
let generated = GeneratedBackendConfig::banded_atomview_build_if_missing_release_tcc();
let generated = GeneratedBackendConfig::banded_atomview_build_if_missing_release_gcc();
let generated = GeneratedBackendConfig::banded_atomview_build_if_missing_release_zig();
```

Production runs can require a prebuilt artifact:

```rust
use RustedSciThe::numerical::BVP_Damp::generated_solver_handoff::{
    AotBuildPolicy, GeneratedBackendConfig,
};

let generated = GeneratedBackendConfig::banded_atomview_for_repeated_solves()
    .with_aot_build_policy(AotBuildPolicy::RequirePrebuilt);
```

`RequirePrebuilt` surfaces a typed error when the artifact is missing. That is usually preferable to an unexpected compilation step in production.

## 12. Practical backend recommendations

When writing a new BVP model, start with Lambdify. It is the simplest way to check equations, boundary conditions, bounds, and initial guesses.

If your model already exists as Rust code, use the Numerical path. Choose `new_numeric_fd_with_options` when maintaining an analytical Jacobian is not worth it, and choose `new_numeric_with_jacobian_options` when the continuous `df/dy` is available and you want to avoid finite-difference noise. For very stiff or badly scaled symbolic models, Lambdify/AOT may still be preferable because they give a mechanically derived Jacobian from the `Expr` source of truth.

If the problem is large and symbolic, test AOT with `tcc` alongside Lambdify rather than waiting until a repeated-solve use case appears. On the corrected `AtomView` path, TCC AOT can already match or beat Lambdify in honest cold end-to-end measurements. For local-stencil BVP/PDE-like problems, test Banded early; for irregular sparse coupling, Sparse is safer. For a single small solve, still begin with Lambdify because build/bootstrap can easily consume any runtime saving.

If the model will be solved repeatedly, distinguish artifact creation from
artifact use explicitly: create or validate it with `BuildIfMissing`, and use
`RequirePrebuilt` in the measured/production loop. On the measured Damped
Banded combustion-1000 case this warm route is about `1.09x` faster per solve
than Lambdify, but requires several reuses before the first build is amortized.

When selecting an AOT toolchain, `tcc` is now the sensible first candidate: it delivered the best practical cold profile in our large BVP runs. `gcc`, Zig, and the Rust backend remain important platform-dependent alternatives and controls, but should not silently be treated as the performance default.

For chunking, do not assume parallel is better. Measure. Inspect callback value time, matrix assembly time, linear solve time, total solver time, actual jobs, and fallback decisions. Current measurements establish both facts: real parallel execution exists, and on Banded `n_steps = 3000` it materially accelerates callback values; full cold wall clock can nevertheless remain nearly unchanged. Therefore `Auto`, rather than forced `chunk4`, is the appropriate library default.

## 13. Reading statistics and story tables

`solver.get_statistics()` returns counters and timers. This interface is available on both the Damped and Frozen solvers. Frozen deliberately remains a different nonlinear strategy, but it now reports the same generated-backend facts needed for honest Lambdify/AOT comparisons: selected backend, symbolic frontend, handoff preparation stages, callback diagnostics, Newton iterations, Jacobian rebuilds, and linear solves. In BVP story tests, the most useful entries are:

`linear_ms` for Newton linear solves, `jac_ms` for Jacobian preparation/evaluation at solver level, `fun_ms` for residual evaluation, `Callback Residual Values` and `Callback Jacobian Values` for low-level callback arithmetic, `Callback Jacobian Matrix Assembly` for matrix construction from values, and counters such as `iters`, `linsys`, and `jac_re`.

For generated AOT routes the same statistics object also carries a `diagnostics` map. The most useful keys are `generated.selected_backend`, `generated.handoff.initial_generate_wall_ms`, `generated.handoff.build_policy_wall_ms`, `generated.handoff.post_build_rebind_wall_ms`, `aot.runtime.execution_policy`, `aot.runtime.parallel_requested`, `aot.runtime.residual.actual_jobs`, `aot.runtime.residual.fallback_reason`, `aot.runtime.sparse_jacobian.actual_jobs`, `aot.runtime.sparse_jacobian.fallback_reason`, and the corresponding `work_per_job` / `work_per_chunk` entries. These fields are intentionally runtime facts, not just requested configuration: if `Auto` or a threshold policy falls back to sequential execution, the statistics will say so directly.

Older tables may still show `Symbolic Operations`; newer code also exports the same timer as `Backend Preparation`. The alias is deliberate. In a Lambdify/AOT route the stage is genuinely symbolic/codegen preparation, while in Numerical path it is numeric discretization and backend handoff rather than symbolic work.

The current combustion BVP measurements suggest the following working picture. For `n_steps = 200`, cold AOT can still lose time to preparation. For `n_steps = 1000`, after switching to `AtomView`, cold TCC AOT beat Lambdify in the recorded comparison on both Sparse and Banded routes; in a separate cooldown-controlled warm Damped Banded run, strict prebuilt TCC also beat Lambdify by about `7.9%`. For Banded `n_steps = 3000`, the production `AtomView` frontend preserved the solution to roundoff and placed TCC AOT at roughly equal or better full cold wall-clock time, although the Lambdify row in that experiment was noisy. Banded remains the natural candidate for a narrow Jacobian band, while chunking should be selected by `Auto` or by a dedicated measurement: it improves hot callback work but does not guarantee proportional whole-solve acceleration.

These are engineering conclusions, not marketing promises. Re-run the story tests on your own hardware when changing toolchain, thread count, mesh size, or equation structure.

## 14. Runnable guide examples

The `examples/` directory contains small executable companions to this chapter. Each example solves a boundary-value problem with an analytical reference solution, so it demonstrates not only the configuration syntax but also an immediate correctness check. The Damped examples cover the pure numerical route, symbolic Lambdify and generated AOT; the Frozen examples cover its supported symbolic Lambdify and AOT routes. The separate Frozen numerical-route note explains why a closure-defined problem should be routed through Damped rather than represented by an empty symbolic placeholder.

```powershell
cargo run --example bvp_damped_numerical_guide
cargo run --example bvp_damped_lambdify_guide
cargo run --example bvp_damped_aot_guide
cargo run --example bvp_frozen_lambdify_guide
cargo run --example bvp_frozen_aot_guide
cargo run --example bvp_frozen_numerical_route_guide
```

The AOT examples deliberately demonstrate both phases of the artifact lifecycle: a first `BuildIfMissing` solve followed by a `RequirePrebuilt` solve. They require `tcc` on `PATH`; if it is not installed, they report that the AOT demonstration was skipped instead of disguising a toolchain failure as a solver failure.

## 15. Where to look next

Basic correctness tests live in `BVP_Damp_tests.rs`. Heavier story and performance runs live in `BVP_Damp_tests3.rs` and `BVP_Damp_tests4.rs`. The running story-test ledger is [`BVP_DAMP_STORY_TESTS.md`](BVP_DAMP_STORY_TESTS.md).

The lower-level codegen/performance layer is documented in [`../../symbolic/codegen/tests/BVP_CODEGEN_STORY_TESTS.md`](../../symbolic/codegen/tests/BVP_CODEGEN_STORY_TESTS.md). Use it when you need to distinguish solver-loop cost from generated callback cost.

Prefer `try_solve()` in production code. The older `solve()` wrapper is kept for compatibility and panics on failure, while `try_solve()` returns typed errors that can be logged, retried, or surfaced to users.