Local-Global Solver API

docs/make.jl

@@ -47,6 +47,7 @@ bibtex_plugin = CitationBibliography(
             "topics/operators.md",
             "topics/couplers.md",
             "topics/time-integration.md",
+            "topics/nonlinear-solver.md",
         ],
         "How-to guides" => [
             "Overview" => "howto/index.md",

docs/src/api-reference/solver.md

@@ -14,6 +14,7 @@ SchurComplementLinearSolver
 
 ```@docs
 NewtonRaphsonSolver
+MultiLevelNewtonRaphsonSolver
 ```
 
 

docs/src/assets/references.bib

@@ -272,6 +272,7 @@ @article{God:1959:dmn
   year={1959},
   publisher={Russian Academy of Sciences, Steklov Mathematical Institute of Russian~…}
 }
+
 @article{Ges:1989:ote,
   author={Geselowitz, D.B.},
   journal={Proceedings of the IEEE},
@@ -282,6 +283,7 @@ @article{Ges:1989:ote
   pages={857-876},
   doi={10.1109/5.29327}
 }
+
 @article{PotDubRicVinGul:2006:cmb,
   author={Potse, Mark and Dube, Bruno and Richer, Jacques and Vinet, Alain and Gulrajani, Ramesh M.},
   journal={IEEE Transactions on Biomedical Engineering}, 
@@ -292,6 +294,7 @@ @article{PotDubRicVinGul:2006:cmb
   pages={2425-2435},
   doi={10.1109/TBME.2006.880875}
 }
+
 @article{OgiBalPer:2023:seats,
   author = {Ogiermann, Dennis and Perotti, Luigi E. and Balzani, Daniel},
   journal = {International Journal for Numerical Methods in Biomedical Engineering},
@@ -302,3 +305,14 @@ @article{OgiBalPer:2023:seats
   pages = {e3670},
   doi = {https://doi.org/10.1002/cnm.3670},
 }
+
+@article{RabSanHsu:1979:mna,
+  title={A multilevel Newton algorithm with macromodeling and latency for the analysis of large-scale nonlinear circuits in the time domain},
+  author={Rabbat, N and Sangiovanni-Vincentelli, Alberto and Hsieh, Hsueh},
+  journal={IEEE Transactions on Circuits and Systems},
+  volume={26},
+  number={9},
+  pages={733--741},
+  year={1979},
+  publisher={IEEE}
+}

docs/src/topics/nonlinear-solver.md

@@ -0,0 +1,97 @@
+# Nonlinear Solver
+
+## Multi-Level Newton-Raphson
+
+A quadratically convergent Newton-Raphson scheme has been proposed by [RabSanHsu:1979:mna](@citet).
+Let us assume we have a block-nonlinear problem with unknowns $\hat{\bm{u}}$ and $\hat{\bm{q}}$ of the following form:
+
+```math
+\begin{aligned}
+    \bm{\hat{f}_{\textrm{G}}}(\hat{\bm{u}}, \hat{\bm{q}}) &= 0 \\
+    \bm{\hat{f}_{\textrm{L}}}(\hat{\bm{u}}, \hat{\bm{q}}) &= 0
+\end{aligned}
+```
+
+where solving $\bm{\hat{f}_{\textrm{L}}}(\hat{\bm{u}}, \hat{\bm{q}}) = 0$ is easy to solve for fixed $\hat{\bm{u}}$.
+If we can enforce this constraint, then we can rewrite the first equation by *implicit function theorem* as:
+```math
+\bm{\hat{f}_{\textrm{G}}}(\hat{\bm{u}}, \hat{\bm{q}}(\hat{\bm{u}})) = 0
+```
+Solving this modified problem with a Newton-Raphson algorithm requires a linearization around $\hat{\bm{u}}$, such that we have to solve at each Newton step the following linear problem
+```math
+\left( \frac{\partial  \bm{\hat{f}_{\textrm{G}}} }{\partial \hat{\bm{u}}} + \frac{\partial \bm{\hat{f}_{\textrm{G}}} }{\partial \hat{\bm{q}}} \frac{\mathrm{d} \hat{\bm{q}} }{\mathrm{d} \hat{\bm{u}}} \right) \Delta \hat{\bm{u}} = -\bm{\hat{f}_{\textrm{G}}}(\hat{\bm{u}}, \hat{\bm{q}}(\hat{\bm{u}}))
+```
+In the continuum mechanics community the system matrix is also known as the *consistent linearization*.
+
+The last missing piece the implicit function part for the system matrix, which is determined by solving an additional linear system:
+```math
+\frac{\partial \bm{\hat{f}_{\textrm{L}}}(\hat{\bm{u}}^{i}, \hat{\bm{q}}^{i}) }{\partial \hat{\bm{q}}} \frac{\mathrm{d} \hat{\bm{q}} }{\mathrm{d} \hat{\bm{u}}} = -\frac{\partial \bm{\hat{f}_{\textrm{L}}}(\hat{\bm{u}}^{i}, \hat{\bm{q}}^{i}) }{\partial \hat{\bm{u}}}
+```
+
+### Using finite-element structure
+
+Time discretization schemes applied to finite element semi-discretizations with $L_2$ variables (called *internal variables*) usually lead to block-nonlinear problems with *local-global structure*, as described above.
+This is commonly found in continuum mechanics problems.
+The local-global structure is simply a result of the algebraic decoupling of the internal variables, as they are associated with the quadrature points.
+For the resulting nonlinear form of the space-time discretization with field unknowns $u$ and internal unknowns $q = (q_1, ... q_{nqp})$, we can write the finite element discretization formally as
+
+```math
+\begin{aligned}
+    f_G(u,q,p,t) =& ... \\
+    f_{Q}(u,q,p,t) =& ... \\
+\end{aligned}
+```
+
+TODO picture with the fundamental decomposition operators
+
+The internal unknowns $q$ are located at the quadrature points which implies the following structure
+
+```math
+\begin{aligned}
+    f_{Q_1}(u,q_1,p,t) =& ... \\
+    f_{Q_2}(u,q_2,p,t) =& ... \\
+    \vdots \\
+    f_{Q_{nqp}}(u,q_{nqp},p,t) = &... \\
+\end{aligned}
+```
+
+### Example: Creep Test of 1D Linear Viscoelasticity
+
+A simple linear viscoelastic material model in 1D in weak form is:
+
+```math
+\begin{aligned}
+                0 =& \int_{\Omega} (E_0 \partial_x u(x,t) + E_1 (\partial_x u(x,t) + q(x,t))) \cdot \partial_x \delta u(x) \textrm{d}x + Neumann \\
+\partial_t q(x,t) =& \frac{E_1}{\eta_1} (\partial_x u(x,t)-q(x,t))
+\end{aligned}
+```
+
+Assuming we have a single 1D element $\Omega = [-1,1]$ with linear ansatz functions and Gauss-Legendre quadrature (i.e. 2 points) the Neumann conditon $\partial_x u(1,t) = 1$ and the Dirichlet condition $u(-1,t) = 0$, then applying a Galerkin semi-discretization yields the following linear DAE in mass matrix form:
+
+```math
+\begin{aligned}
+            0   &=  \tilde{u}_1 \\
+            0   &=  0.5\left(-(E_0 + E_1) \tilde{u}_1 + (E_0 + E_1) \tilde{u}_2 - E_1 \tilde{q}_1 - E_1 \tilde{q}_2 + 1\right) \\
+d_t \tilde{q}_1 &= \frac{E_1}{\eta_1} (-\tilde{u}_1+\tilde{u}_2-\tilde{q}_1) \\
+d_t \tilde{q}_2 &= \frac{E_1}{\eta_1} (-\tilde{u}_1+\tilde{u}_2-\tilde{q}_2)
+\end{aligned}
+```
+or, after condensing the first equation,
+```math
+\begin{aligned}
+            0   &= 0.5\left((E_0 + E_1) \tilde{u}_2(t) - E_1 \tilde{q}_1(t) - E_1 \tilde{q}_2(t) + 1\right) \\
+d_t \tilde{q}_1(t) &= \frac{E_1}{\eta_1} (-\tilde{u}_1+\tilde{u}_2(t)-\tilde{q}_1(t)) \\
+d_t \tilde{q}_2(t) &= \frac{E_1}{\eta_1} (-\tilde{u}_1+\tilde{u}_2(t)-\tilde{q}_2(t))
+\end{aligned}
+```
+
+Applying the Backward Euler in time to this system we obtain the ,,nonlinear'' system
+```math
+\begin{aligned}
+ 0 &= 0.5\left((E_0 + E_1) \hat{\tilde{u}}^n_2 - E_1 \hat{\tilde{q}}_1 - E_1 \hat{\tilde{q}}_2 + f(t_n)\right) &= f_G(\hat{\tilde{u}}_2, \hat{\tilde{q}}_1, \hat{\tilde{q}}_2) \\
+ 0 &= \hat{\tilde{q}}^n_1 - \hat{\tilde{q}}^{n-1}_1 - \Delta t_n \frac{E_1}{\eta_1} (-\hat{\tilde{u}}_1+\hat{\tilde{u}}_2-\hat{\tilde{q}}_1) &= f_{Q_1}(\hat{\tilde{u}}_2, \hat{\tilde{q}}_1) \\
+ 0 &= \hat{\tilde{q}}^n_2 - \hat{\tilde{q}}^{n-1}_2 - \Delta t_n \frac{E_1}{\eta_1} (-\hat{\tilde{u}}_1+\hat{\tilde{u}}_2-\hat{\tilde{q}}_2) &= f_{Q_2}(\hat{\tilde{u}}_2, \hat{\tilde{q}}_2) \\
+\end{aligned}
+```
+
+where we can observe that the internal problems are decoupled. The system can now be solved with the multi-level Newton-Raphson as described above.