conspire/domain/partition/feti/pcg/
mod.rs

1use crate::{
2    domain::feti::{
3        THREADS,
4        dual::{
5            dual_action, dual_operator, dual_precondition, dual_precondition_dirichlet,
6            dual_precondition_scaled_dirichlet,
7        },
8        dual_primal::{coarse::Coarse, rigid_projector::RigidProjector},
9        subdomain::Subdomain,
10    },
11    math::{
12        Scalar, Tensor, Vector,
13        optimize::{Krylov, KrylovError, KrylovMethod},
14    },
15};
16
17/// Solves the dual (interface) problem `(F + C S_pp^-1 C^T) . lambda = rhs`
18/// by conjugate gradients, Dirichlet-preconditioned. For a symmetric positive
19/// definite tangent the operator is SPD given the corners are pinned, so this
20/// needs no projection against a rigid-body null space the way plain FETI
21/// would; a nonsymmetric tangent needs GMRES instead. Dirichlet is the default over lumped: its
22/// condition-number bound is near mesh-independent (`1 + log(H/h)^2`) where
23/// lumped's degrades with the subdomain-to-mesh-size ratio, at the price of
24/// one local interior solve per subdomain per iteration, so it wins as soon
25/// as a subdomain has more than a handful of elements per edge, the
26/// realistic regime. `dual_precondition` (lumped) stays available for the
27/// pathologically-small-subdomain case.
28pub(crate) fn projected_pcg<B>(
29    subdomains: &[Subdomain<B>],
30    coarse: &Coarse,
31    rhs: &Vector,
32) -> Result<Vector, KrylovError>
33where
34    B: Sync,
35{
36    projected_pcg_with(
37        subdomains,
38        coarse,
39        rhs,
40        Preconditioner::Dirichlet,
41        Krylov::default().rel_tol,
42        KrylovMethod::default(),
43    )
44}
45
46/// Which FETI-DP preconditioner the dual PCG applies.
47///
48/// `Dirichlet` is the default: its condition-number bound is near
49/// mesh-independent, at the price of one interior solve per subdomain per
50/// iteration over `Lumped`'s cheap matvec — see `projected_pcg`'s own doc
51/// for the tradeoff in full.
52#[derive(Clone, Copy, Debug, PartialEq, Eq)]
53pub enum Preconditioner {
54    /// `sum_s B_s K_dd,s B_s^T` — a local matvec in place of a local solve.
55    Lumped,
56    /// `sum_s B_b,s S_s B_b,s^T` — near mesh-independent convergence.
57    Dirichlet,
58}
59
60pub(crate) fn projected_pcg_with<B>(
61    subdomains: &[Subdomain<B>],
62    coarse: &Coarse,
63    rhs: &Vector,
64    preconditioner: Preconditioner,
65    rel_tol: Scalar,
66    method: KrylovMethod,
67) -> Result<Vector, KrylovError>
68where
69    B: Sync,
70{
71    Krylov {
72        rel_tol,
73        method,
74        ..Krylov::default()
75    }
76    .solve(
77        |lambda| dual_operator(subdomains, lambda, coarse, THREADS),
78        |lambda: &Vector| match preconditioner {
79            Preconditioner::Lumped => dual_precondition(subdomains, lambda, THREADS),
80            Preconditioner::Dirichlet => dual_precondition_dirichlet(subdomains, lambda, THREADS),
81        },
82        rhs,
83    )
84}
85
86/// Solves the dual problem of classical FETI, `F lambda - G alpha = d` with
87/// `G^T lambda = e`, for the multipliers and the rigid-body amplitudes.
88///
89/// The multipliers are `lambda_0 + mu`, with `lambda_0 = G (G^T G)^-1 e`
90/// satisfying the constraint and `mu` found by the same Krylov solve as
91/// FETI-DP on `P F P mu = P (d - F lambda_0)`, preconditioned by `P M P`, with
92/// the Dirichlet preconditioner scaled by multiplicity. The
93/// projector `P` keeps every iterate in the space where `G^T mu = 0`, which is
94/// what makes the floating subdomains' local solves consistent.
95pub(crate) fn rigid_projected_pcg<B>(
96    subdomains: &[Subdomain<B>],
97    projector: &RigidProjector,
98    d: &Vector,
99    e: &Vector,
100    preconditioner: Preconditioner,
101    rel_tol: Scalar,
102    method: KrylovMethod,
103) -> Result<(Vector, Vector), KrylovError>
104where
105    B: Sync,
106{
107    let particular = projector.particular(e, d.len());
108    let rhs = projector.project(&(d.clone() - dual_action(subdomains, &particular, THREADS)));
109    let correction = Krylov {
110        rel_tol,
111        method,
112        ..Krylov::default()
113    }
114    .solve(
115        |mu: &Vector| projector.project(&dual_action(subdomains, &projector.project(mu), THREADS)),
116        |residual: &Vector| {
117            let residual = projector.project(residual);
118            projector.project(&match preconditioner {
119                Preconditioner::Lumped => dual_precondition(subdomains, &residual, THREADS),
120                Preconditioner::Dirichlet => {
121                    dual_precondition_scaled_dirichlet(subdomains, &residual, THREADS)
122                }
123            })
124        },
125        &rhs,
126    )?;
127    let lambda = particular + projector.project(&correction);
128    let alpha = projector.amplitudes(&(dual_action(subdomains, &lambda, THREADS) - d.clone()));
129    Ok((lambda, alpha))
130}
131
132/// Recovers each subdomain's full local solution (corner and dual DOFs
133/// together) once the coarse (corner) solution and multipliers are known:
134/// `u_d,s = K_dd,s^-1 (f_d,s - K_dp,s . u_p - B_s^T . lambda)`, with `u_p`
135/// gathered into this subdomain's raw local corner positions and the
136/// coupling term `K_dd,s^-1 K_dp,s . u_p = dual_map_s . u_p` already
137/// available from `condense()`.
138pub(crate) fn primal_recovery<B>(
139    subdomains: &[Subdomain<B>],
140    local_forces: &[Vector],
141    corner_solution: &Vector,
142    lambda: &Vector,
143) -> Vec<Vector> {
144    subdomains
145        .iter()
146        .zip(local_forces.iter())
147        .map(|(subdomain, local_force)| {
148            let multiplier_rhs = subdomain
149                .interface()
150                .apply_transpose(lambda, subdomain.num_local());
151            let combined_rhs = local_force - &multiplier_rhs;
152            let dual_solution = subdomain.local_solve(&combined_rhs);
153            let primal_local = subdomain.gather_primal(corner_solution);
154            let coupling_correction =
155                subdomain.scatter_dual(&(subdomain.dual_map() * &primal_local));
156            dual_solution - coupling_correction + subdomain.scatter_primal(&primal_local)
157        })
158        .collect()
159}
conspire/domain/partition/feti/pcg/mod.rs

conspire/domain/partition/feti/pcg/
mod.rs
