sim_lib_femm_function/
implementation.rs

1#![forbid(unsafe_code)]
2//! Callable wrapper that evaluates a model to a quantity, field, or solution.
3//!
4//! Defines the call request, output query, evaluation result, and the callable
5//! payload that turns a FEMM model into a runtime function of its parameters.
6
7use std::{any::Any, sync::Arc};
8
9use sim_kernel::{
10    ClassId, ClassRef, Cx, DefaultFactory, Expr, Factory, Object, ObjectEncode, ObjectEncoding,
11    Result as KernelResult, Symbol, Value,
12};
13use sim_lib_femm_core::{FemmError, FemmLimits, FemmResult, ParamSet, StableId, value_as_f64};
14use sim_lib_femm_field::{Field, Projection, field_as_func};
15use sim_lib_femm_material::{BoundaryKind, Source};
16use sim_lib_femm_mesh::FemmModel;
17use sim_lib_femm_post::{Excitation, FemmSolution, QuantitySpec, quantity};
18use sim_lib_femm_solve::{GradientTrust, SolveCertificate, SteadySolve, solve_steady};
19use sim_lib_numbers_func::{Func, FuncMetadata};
20
21/// One evaluation request against a model: which parameters, output, and limits.
22///
23/// Couples a parameter binding with the [`OutputQuery`] to compute and the
24/// solver [`FemmLimits`]; `want_grad` names parameters whose sensitivities the
25/// caller also wants. See the [crate README](index.html).
26#[derive(Clone, Debug)]
27pub struct FemmCall {
28    /// The parameter binding the model is evaluated at.
29    pub params: ParamSet,
30    /// The output to compute from the solved model.
31    pub query: OutputQuery,
32    /// Parameters to also report sensitivities for, if any.
33    pub want_grad: Option<Vec<Symbol>>,
34    /// Solver budget and tolerances for this evaluation.
35    pub limits: FemmLimits,
36}
37
38/// The kind of output an evaluation produces from a solved model.
39#[derive(Clone, Debug)]
40pub enum OutputQuery {
41    /// A scalar quantity reduced from the solution (energy, flux, capacitance).
42    Quantity(QuantitySpec),
43    /// A projected field (potential or a derived component) over the mesh.
44    Field(Projection),
45    /// The full solved model solution as an opaque value.
46    Solution,
47}
48
49/// The result of evaluating a model: the output value plus optional gradient.
50#[derive(Clone, Debug)]
51pub struct FemmEval {
52    /// The computed output value (scalar, field, or solution).
53    pub value: Value,
54    /// Per-parameter sensitivities, when a gradient was requested.
55    pub gradient: Option<Vec<(Symbol, f64)>>,
56    /// Diagnostics emitted while solving and reducing the output.
57    pub diagnostics: Vec<sim_kernel::Diagnostic>,
58}
59
60/// Quantity value, certificate, and optional total gradient for a completed solve.
61#[derive(Clone, Debug)]
62pub struct QualityAnswer {
63    /// Scalar value of the requested quantity.
64    pub value: f64,
65    /// Certificate describing residual, convergence, and gradient trust.
66    pub certificate: SolveCertificate,
67    /// Gradient values and trust tag when a parameter list is supplied.
68    pub gradient: Option<(Vec<f64>, GradientTrust)>,
69}
70
71/// The opaque payload carried by a model-derived runtime function.
72///
73/// Recorded in a [`Func`]'s metadata so a differentiator can recover the model,
74/// its free variables, and the queried output to build an adjoint pass.
75#[derive(Clone)]
76pub struct FemmFuncPayload {
77    /// The model the function evaluates.
78    pub model: FemmModel,
79    /// The model inputs treated as the function's free variables.
80    pub vars: Vec<Symbol>,
81    /// The output the function returns.
82    pub query: OutputQuery,
83}
84
85impl Object for FemmFuncPayload {
86    fn display(&self, _cx: &mut Cx) -> KernelResult<String> {
87        Ok(format!(
88            "#<femm-payload model={} query={}>",
89            self.model.id.0,
90            describe_query(&self.query)
91        ))
92    }
93
94    fn as_any(&self) -> &dyn Any {
95        self
96    }
97}
98
99impl sim_kernel::ObjectCompat for FemmFuncPayload {
100    fn class(&self, cx: &mut Cx) -> KernelResult<ClassRef> {
101        if let Some(class) = cx
102            .registry()
103            .class_by_symbol(&Symbol::qualified("femm", "FuncPayload"))
104        {
105            return Ok(class.clone());
106        }
107        DefaultFactory.class_stub(ClassId(33), Symbol::qualified("femm", "FuncPayload"))
108    }
109    fn as_expr(&self, cx: &mut Cx) -> KernelResult<Expr> {
110        sim_citizen::constructor_expr(cx, self)
111    }
112    fn as_object_encoder(&self) -> Option<&dyn ObjectEncode> {
113        Some(self)
114    }
115}
116
117impl ObjectEncode for FemmFuncPayload {
118    fn object_encoding(&self, _cx: &mut Cx) -> KernelResult<ObjectEncoding> {
119        Ok(ObjectEncoding::Constructor {
120            class: func_payload_class_symbol(),
121            args: payload_constructor_args(self),
122        })
123    }
124}
125
126impl sim_citizen::Citizen for FemmFuncPayload {
127    fn citizen_symbol() -> Symbol {
128        func_payload_class_symbol()
129    }
130
131    fn citizen_version() -> u32 {
132        1
133    }
134
135    fn citizen_arity() -> usize {
136        3
137    }
138
139    fn citizen_fields() -> &'static [&'static str] {
140        &["model_id", "query", "vars"]
141    }
142}
143
144fn func_payload_class_symbol() -> Symbol {
145    Symbol::qualified("femm", "FuncPayload")
146}
147
148fn payload_constructor_args(payload: &FemmFuncPayload) -> Vec<Expr> {
149    vec![
150        Expr::Symbol(Symbol::new("v1")),
151        int_expr(payload.model.id.0),
152        Expr::String(describe_query(&payload.query)),
153        Expr::List(
154            payload
155                .vars
156                .iter()
157                .map(|name| Expr::String(name.to_string()))
158                .collect(),
159        ),
160    ]
161}
162
163fn int_expr(value: impl ToString) -> Expr {
164    Expr::Number(sim_kernel::NumberLiteral {
165        domain: Symbol::qualified("citizen", "int"),
166        canonical: value.to_string(),
167    })
168}
169
170/// Something that can be evaluated as a FEMM function of its parameters.
171pub trait FemmCallable {
172    /// Evaluates the callable for one [`FemmCall`], returning its output.
173    fn eval(&self, cx: &mut Cx, call: FemmCall) -> FemmResult<FemmEval>;
174}
175
176/// A [`FemmCallable`] that solves a concrete model on each evaluation.
177///
178/// Resolves defaults for any unbound inputs, runs the steady solve, and reduces
179/// the solution to the requested [`OutputQuery`].
180#[derive(Clone)]
181pub struct ModelCallable {
182    /// The model solved on each call.
183    pub model: FemmModel,
184}
185
186impl ModelCallable {
187    fn resolve_params(&self, params: &ParamSet) -> FemmResult<ParamSet> {
188        let mut entries = params.entries.clone();
189        for input in &self.model.inputs {
190            if entries.iter().all(|(name, _)| name != &input.name) {
191                if let Some(default) = &input.default {
192                    entries.push((input.name.clone(), default.clone()));
193                } else {
194                    return Err(FemmError::UnknownFemmParameter(input.name.to_string()));
195                }
196            }
197        }
198        Ok(ParamSet::new(entries))
199    }
200
201    fn solve_solution(
202        &self,
203        cx: &mut Cx,
204        params: &ParamSet,
205        limits: &FemmLimits,
206    ) -> FemmResult<Arc<FemmSolution>> {
207        let resolved = self.resolve_params(params)?;
208        solve_steady(cx, &self.model, &resolved, limits, None).map(|out| out.solution)
209    }
210}
211
212impl FemmCallable for ModelCallable {
213    fn eval(&self, cx: &mut Cx, call: FemmCall) -> FemmResult<FemmEval> {
214        let resolved = self.resolve_params(&call.params)?;
215        match call.query {
216            OutputQuery::Quantity(QuantitySpec::Custom { expr, .. }) => {
217                let value = sim_lib_femm_geometry::eval_expr_f64(cx, &expr, &resolved, &[])?;
218                Ok(FemmEval {
219                    value: cx
220                        .factory()
221                        .number_literal(Symbol::qualified("numbers", "f64"), value.to_string())
222                        .map_err(|err| FemmError::SensitivityUnavailable(err.to_string()))?,
223                    gradient: None,
224                    diagnostics: Vec::new(),
225                })
226            }
227            OutputQuery::Quantity(spec) => {
228                let solution = self.solve_solution(cx, &resolved, &call.limits)?;
229                let excitation = resolve_excitation(cx, &self.model, &resolved, &spec)?;
230                let scalar = quantity(&solution, &spec, &excitation)?;
231                Ok(FemmEval {
232                    value: cx
233                        .factory()
234                        .number_literal(Symbol::qualified("numbers", "f64"), scalar.to_string())
235                        .map_err(|err| FemmError::SensitivityUnavailable(err.to_string()))?,
236                    gradient: None,
237                    diagnostics: Vec::new(),
238                })
239            }
240            OutputQuery::Field(projection) => {
241                let solution = self.solve_solution(cx, &resolved, &call.limits)?;
242                let field = Field::new(solution, projection);
243                Ok(FemmEval {
244                    value: cx
245                        .factory()
246                        .opaque(Arc::new(field))
247                        .map_err(|err| FemmError::SensitivityUnavailable(err.to_string()))?,
248                    gradient: None,
249                    diagnostics: Vec::new(),
250                })
251            }
252            OutputQuery::Solution => {
253                let solution = self.solve_solution(cx, &resolved, &call.limits)?;
254                Ok(FemmEval {
255                    value: cx
256                        .factory()
257                        .opaque(solution)
258                        .map_err(|err| FemmError::SensitivityUnavailable(err.to_string()))?,
259                    gradient: None,
260                    diagnostics: Vec::new(),
261                })
262            }
263        }
264    }
265}
266
267/// Resolves the [`Excitation`] a derived quantity is evaluated against.
268///
269/// Inductance and flux linkage read the driving current of the named circuit
270/// coil source; capacitance reads the applied potential of the named conductor
271/// (its Dirichlet boundary). Quantities that do not depend on an excitation
272/// resolve to [`Excitation::none`]. A source the model does not define leaves
273/// the excitation unset, so [`quantity`] reports the precise missing-drive
274/// error rather than a silent wrong value.
275pub fn resolve_excitation(
276    cx: &mut Cx,
277    model: &FemmModel,
278    params: &ParamSet,
279    spec: &QuantitySpec,
280) -> FemmResult<Excitation> {
281    match spec {
282        QuantitySpec::Inductance { circuit } | QuantitySpec::FluxLinkage { circuit } => {
283            Ok(coil_current(cx, model, params, circuit)?
284                .map(Excitation::with_current)
285                .unwrap_or_else(Excitation::none))
286        }
287        QuantitySpec::Capacitance { conductor } => {
288            Ok(conductor_potential(cx, model, params, conductor)?
289                .map(Excitation::with_potential)
290                .unwrap_or_else(Excitation::none))
291        }
292        _ => Ok(Excitation::none()),
293    }
294}
295
296fn coil_current(
297    cx: &mut Cx,
298    model: &FemmModel,
299    params: &ParamSet,
300    circuit: &Symbol,
301) -> FemmResult<Option<f64>> {
302    for source in &model.sources {
303        if let Source::CircuitCoil { name, current, .. } = source
304            && name == circuit
305        {
306            return sim_lib_femm_geometry::eval_expr_f64(cx, current, params, &[]).map(Some);
307        }
308    }
309    Ok(None)
310}
311
312fn conductor_potential(
313    cx: &mut Cx,
314    model: &FemmModel,
315    params: &ParamSet,
316    conductor: &Symbol,
317) -> FemmResult<Option<f64>> {
318    for boundary in &model.boundaries {
319        if &boundary.name == conductor && matches!(boundary.kind, BoundaryKind::Dirichlet) {
320            return sim_lib_femm_geometry::eval_expr_f64(cx, &boundary.value, params, &[])
321                .map(Some);
322        }
323    }
324    Ok(None)
325}
326
327/// Returns the requested quantity and the certificate for a completed solve.
328///
329/// Passing `Some(params)` for `wrt` also computes a total finite-difference
330/// gradient and annotates the returned certificate with its trust level.
331/// Passing `None` skips gradient work.
332pub fn quality(
333    cx: &mut Cx,
334    solve: &SteadySolve,
335    quantity_spec: &QuantitySpec,
336    wrt: Option<&[Symbol]>,
337) -> FemmResult<QualityAnswer> {
338    let excitation = resolve_excitation(cx, &solve.model, &solve.solution.params, quantity_spec)?;
339    let value = quantity(&solve.solution, quantity_spec, &excitation)?;
340    let mut certificate = solve.certificate.clone();
341    let gradient = match wrt {
342        None => None,
343        Some(params) => {
344            let (values, trust) =
345                finite_difference_quality_gradient(cx, solve, quantity_spec, params)?;
346            certificate.set_gradient_trust(trust.clone());
347            Some((values, trust))
348        }
349    };
350    Ok(QualityAnswer {
351        value,
352        certificate,
353        gradient,
354    })
355}
356
357fn finite_difference_quality_gradient(
358    cx: &mut Cx,
359    solve: &SteadySolve,
360    quantity_spec: &QuantitySpec,
361    wrt: &[Symbol],
362) -> FemmResult<(Vec<f64>, GradientTrust)> {
363    let callable = ModelCallable {
364        model: solve.model.clone(),
365    };
366    let base_params = callable.resolve_params(&solve.solution.params)?;
367    let mut out = Vec::with_capacity(wrt.len());
368    for symbol in wrt {
369        let base_value = base_params
370            .get(symbol)
371            .ok_or_else(|| FemmError::UnknownFemmParameter(symbol.to_string()))?;
372        let x = value_as_f64(cx, base_value)?;
373        if !x.is_finite() {
374            return Err(FemmError::SensitivityUnavailable(format!(
375                "non-finite FEMM parameter {symbol}"
376            )));
377        }
378        let h = fd_step(x);
379        let plus = replace_param_value(cx, &base_params, symbol, x + h)?;
380        let minus = replace_param_value(cx, &base_params, symbol, x - h)?;
381        let plus_value = quality_at_params(cx, &solve.model, plus, quantity_spec)?;
382        let minus_value = quality_at_params(cx, &solve.model, minus, quantity_spec)?;
383        out.push((plus_value - minus_value) / (2.0 * h));
384    }
385    Ok((out, GradientTrust::FiniteDifferenceOnly))
386}
387
388fn quality_at_params(
389    cx: &mut Cx,
390    model: &FemmModel,
391    params: ParamSet,
392    quantity_spec: &QuantitySpec,
393) -> FemmResult<f64> {
394    let solved = solve_steady(cx, model, &params, &FemmLimits::default(), None)?;
395    let excitation = resolve_excitation(cx, model, &params, quantity_spec)?;
396    quantity(&solved.solution, quantity_spec, &excitation)
397}
398
399fn replace_param_value(
400    cx: &mut Cx,
401    params: &ParamSet,
402    name: &Symbol,
403    value: f64,
404) -> FemmResult<ParamSet> {
405    let mut found = false;
406    let mut entries = params.entries.clone();
407    for (symbol, slot) in &mut entries {
408        if symbol == name {
409            *slot = cx
410                .factory()
411                .number_literal(Symbol::qualified("numbers", "f64"), value.to_string())
412                .map_err(|err| FemmError::SensitivityUnavailable(err.to_string()))?;
413            found = true;
414        }
415    }
416    if found {
417        Ok(ParamSet::new(entries))
418    } else {
419        Err(FemmError::UnknownFemmParameter(name.to_string()))
420    }
421}
422
423fn fd_step(value: f64) -> f64 {
424    1.0e-6 * value.abs().max(1.0)
425}
426
427/// Wraps a model as a sim-numbers [`Func`] of the named variables.
428///
429/// The returned function solves the model on call and reduces it to `query`;
430/// its metadata carries a [`FemmFuncPayload`] and an adjoint differentiator
431/// hint so sensitivity analysis can recover the model.
432///
433/// # Examples
434///
435/// ```
436/// use sim_kernel::Symbol;
437/// use sim_lib_femm_fixtures::parallel_plate_capacitor;
438/// use sim_lib_femm_function::{femm_as_func, OutputQuery};
439/// use sim_lib_femm_post::QuantitySpec;
440///
441/// let vars = vec![Symbol::new("gap-mm")];
442/// let func = femm_as_func(
443///     parallel_plate_capacitor(),
444///     vars.clone(),
445///     OutputQuery::Quantity(QuantitySpec::Energy { region: None }),
446/// );
447/// assert_eq!(func.vars, vars);
448/// ```
449pub fn femm_as_func(model: FemmModel, vars: Vec<Symbol>, query: OutputQuery) -> Func {
450    let callable = ModelCallable {
451        model: model.clone(),
452    };
453    let closure_vars = vars.clone();
454    let payload_vars = closure_vars.clone();
455    let closure_query = query.clone();
456    Func {
457        vars,
458        body_cas: None,
459        body_native: Some(Arc::new(move |cx, args| {
460            let params = ParamSet::new(
461                closure_vars
462                    .iter()
463                    .cloned()
464                    .zip(args.iter().cloned())
465                    .collect::<Vec<_>>(),
466            );
467            callable
468                .eval(
469                    cx,
470                    FemmCall {
471                        params,
472                        query: closure_query.clone(),
473                        want_grad: None,
474                        limits: FemmLimits::default(),
475                    },
476                )
477                .map(|out| out.value)
478                .map_err(sim_kernel::Error::from)
479        })),
480        metadata: FuncMetadata {
481            source: Some(Symbol::qualified("femm", "model")),
482            differentiator_hint: Some(Symbol::new("femm-adjoint")),
483            payload: DefaultFactory
484                .opaque(Arc::new(FemmFuncPayload {
485                    model: model.clone(),
486                    vars: payload_vars,
487                    query: query.clone(),
488                }))
489                .ok(),
490        },
491    }
492}
493
494/// Wraps a model's potential field as a sim-numbers [`Func`] over position.
495///
496/// Builds a trivial single-element solution for `model` and exposes its
497/// potential projection as a spatial function, used where a model is consumed
498/// as a field-valued function rather than a parameter-to-scalar map.
499pub fn femm_field_func(model: FemmModel) -> Func {
500    let field = Arc::new(FemmSolution {
501        id: StableId(model.id.0 + 1),
502        model_id: model.id,
503        physics: model.physics.clone(),
504        formulation: model.formulation.clone(),
505        params: ParamSet::default(),
506        mesh: sim_lib_femm_mesh::FemMesh2 {
507            xy: vec![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
508            tri: vec![[0, 1, 2]],
509            elem_region: vec![Symbol::new("air")],
510            edge_boundary: Vec::new(),
511        },
512        u: vec![0.0, 1.0, 1.0],
513        diagnostics: sim_lib_femm_flow::SolveDiagnostics {
514            method: Symbol::new("femm-ptc"),
515            converged: true,
516            iterations: 1,
517            final_residual: 0.0,
518            events: Vec::new(),
519            diagnostics: Vec::new(),
520        },
521    });
522    field_as_func(Field::new(field, Projection::Potential))
523}
524
525pub(crate) fn describe_query(query: &OutputQuery) -> String {
526    match query {
527        OutputQuery::Quantity(QuantitySpec::Custom { name, .. }) => format!("quantity:{name}"),
528        OutputQuery::Quantity(_) => "quantity".to_owned(),
529        OutputQuery::Field(projection) => format!("field:{projection:?}"),
530        OutputQuery::Solution => "solution".to_owned(),
531    }
532}
sim_lib_femm_function/implementation.rs

sim_lib_femm_function/
implementation.rs
