{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "19805282", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:09.664212Z", "iopub.status.busy": "2026-10-02T14:47:09.663973Z", "iopub.status.idle": "2026-10-02T14:47:09.669299Z", "shell.execute_reply": "2026-10-02T14:47:09.668722Z" }, "papermill": { "duration": 0.008722, "end_time": "2026-10-02T14:47:09.669865+00:00", "exception": false, "start_time": "2026-10-02T14:47:09.661143+00:00", "status": "completed" }, "tags": [ "remove-input", "active-ipynb", "remove-output" ] }, "outputs": [], "source": [ "try:\n", " from openmdao.utils.notebook_utils import notebook_mode # noqa: F401\n", "except ImportError:\n", " !python -m pip install openmdao[notebooks]" ] }, { "cell_type": "markdown", "id": "381ad43b", "metadata": { "papermill": { "duration": 0.001411, "end_time": "2026-10-02T14:47:09.672954+00:00", "exception": false, "start_time": "2026-10-02T14:47:09.671543+00:00", "status": "completed" }, "tags": [] }, "source": [ "# Setting Up a Model for Efficient Linear Solves\n", "\n", "There are a number of different features that you can use to control how the linear solves are performed that will have an impact on both the speed and accuracy of the linear solution.\n", "A deeper understanding of how OpenMDAO solves the unified derivatives equations is useful in understanding when to apply certain features, and may also help you structure your model to make the most effective use of these features.\n", "The explanation of OpenMDAO's features for improving linear solver performance are broken up into three sections below:\n", "\n", "(determining_what_kind_of_linear_solver_to_use)=\n", "## Determining What Kind of Linear Solver to Use\n", "\n", "Since total derivatives are computed by solving the unified derivatives equations, there is always some kind of linear solver used by the framework whenever [compute_totals](../features/core_features/working_with_derivatives/approximating_totals.ipynb) is called.\n", "However, the specific type of linear solver that should be used will vary greatly depending on the underlying model structure.\n", "The most basic distinguishing feature of a model that governs what kind of linear solver should be used is the presence of any coupling.\n", "\n", "\n", "## Uncoupled Models\n", "\n", "If you have a completely uncoupled model, then the partial-derivative Jacobian matrix will have a lower-triangular structure.\n", "The resulting linear system can be solved using a block-forward or block-backward substitution algorithm.\n", "Alternatively you could view the solution algorithm as a single iteration of a block Gauss-Seidel algorithm.\n", "In OpenMDAO, the single-pass block Gauss-Seidel algorithm is implemented via the [LinearRunOnce](../features/building_blocks/solvers/linear_runonce.ipynb) solver.\n", "This is the default solver used by OpenMDAO on all [Groups](../features/core_features/working_with_groups/index.ipynb).\n", "\n", "![uncoupled_fwd](images/matrix_figs/uncoupled_fwd.png)\n", "\n", "
The linear system for an uncoupled system can be solved with forward substitution in forward mode
\n", "\n", "If you are using *reverse* mode, then the left-hand side of the unified derivatives equations will be the transpose-Jacobian and will have an upper-triangular structure.\n", "The upper-triangular transpose-Jacobian structure is notable, because it can also be seen in the [n2 diagram](../other_useful_docs/om_command.ipynb) that OpenMDAO can produce.\n", "\n", "![uncoupled_rev](images/matrix_figs/uncoupled_rev.png)\n", "\n", "
The linear system for an uncoupled system can be solved with back substitution in reverse mode.
\n", "\n", "\n", "## Coupled Models\n", "\n", "Coupled models will always have a non-triangular structure to their partial-derivative Jacobian.\n", "In other words, there will be nonzero entries both above and below the diagonal.\n", "\n", "![coupled_fwd](images/matrix_figs/coupled_fwd.png)\n", "\n", "
The linear system of a coupled model must be solved using either a direct factorization or an iterative linear solver.
\n", "\n", "Consequently, these linear systems **cannot** be solved with the [LinearRunOnce](../features/building_blocks/solvers/linear_runonce.ipynb) solver.\n", "There are two basic categories of linear solver that can be used in this situation:\n", "\n", "1. Direct solvers (e.g. [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb))\n", "2. Iterative solvers (e.g. [LinearBlockGS](../features/building_blocks/solvers/linear_block_gs.ipynb), [ScipyKrylov](../features/building_blocks/solvers/scipy_iter_solver.ipynb))\n", "\n", "Direct solvers make use of the Jacobian matrix, assembled in memory, in order to compute an inverse or a factorization that can be used to solve the linear system.\n", "Conversely, iterative linear solvers find the solution to the linear system without ever needing to access the Jacobian matrix directly.\n", "They search for solution vectors that drive the linear residual to 0 using only matrix-vector products.\n", "The key idea is that **some** kind of linear solver is needed when there is coupling in your model.\n", "\n", "Which type of solver is best for your model use is heavily case-dependent and sometimes can be a difficult question to answer absolutely.\n", "However, there are a few rules of thumb that can be used to guide most cases:\n", "\n", "1. Direct solvers are very simple to use, and for smaller problems, are likely to be the best option. The only downside is that the cost of computing the factorization scales is $n^3$, where $n$ is the length of your variable vector, so the compute cost can get out of control. If $n$ < 2000, try this solver first.\n", "2. Iterative solvers are more difficult to use because they do not always succeed in finding a good solution to the linear problem. Often times they require preconditioners in order to be effective. However, with adequate preconditioning, iterative solvers can dramatically outperform direct solvers for even moderate-sized problems. The trade-off you make is computational speed for complexity in getting the solver to work. Iterative solvers can also offer significant memory savings, since there isn't a need to allocate one large matrix for all the partials.\n", "\n", "```{note}\n", "There is a relationship between linear and nonlinear solvers. Any coupling in your model will affect both the linear and nonlinear solves, and thus impact which type of linear and nonlinear solvers you use.\n", "\n", "In the most basic case, an uncoupled model will use the default [NonLinearRunOnce](../features/building_blocks/solvers/nonlinear_runonce.ipynb) and the [LinearRunOnce](../features/building_blocks/solvers/linear_runonce.ipynb) solvers. These *RunOnce* solvers are a special degenerate class of Solver, which can't handle any kind of coupling or implicitness in a model. Any model with coupling will require an iterative nonlinear solver. Any model that requires an iterative nonlinear solver will also need a linear solver other than the default [LinearRunOnce](../features/building_blocks/solvers/linear_runonce.ipynb) solver.\n", "```\n", "\n", "(theory-assembled-vs-matrix-free)=\n", "## Selecting Linear Solver Architecture: Dense, Sparse, or Matrix-Free\n", "\n", "Broadly speaking, there are two classes of linear solver architecture:\n", "\n", "- Assembled Jacobian\n", "- Matrix-free\n", "\n", "At any level of the hierarchy in an OpenMDAO model, you have the option of choosing between these two options.\n", "Simple models will often just use one linear solver architecture at the top of the model hierarchy.\n", "More complex models might use both architectures at different parts of the hierarchy.\n", "At any level of the hierarchy, you can look at the aspects of the components contained within that group in order to\n", "figure out what kind of linear solver structure is needed.\n", "\n", "\n", "## Assembled-Jacobian Problems\n", "\n", "Using an assembled Jacobian means that OpenMDAO will explicitly allocate the memory for the entire\n", "Jacobian matrix up front, and then hold onto that and re-use it to perform linear solves throughout the run.\n", "This has several computational advantages, but the major one is that it helps to reduce framework\n", "overhead for models with deep system hierarchies and large numbers of variables.\n", "\n", "You should strongly consider using an assembled Jacobian if all the components in your model provide\n", "derivatives using the [compute_partials](../features/core_features/working_with_components/explicit_component.ipynb) or\n", "[linearize](../features/core_features/working_with_components/implicit_component.ipynb) methods.\n", "These methods are explicitly computing the elements of that Jacobian matrix, and so it makes sense\n", "to collect them into an actual matrix memory representation.\n", "\n", "Additionally, if your model has a very large hierarchy (i.e. many levels, many components,\n", "many variables) then an assembled Jacobian will likely offer a significant performance advantage.\n", "The reason that large models benefit is that without the assembled Jacobian, OpenMDAO must\n", "recursively loop over each level of the hierarchy, each component, and each variable in order to\n", "compute Jacobian-vector products.\n", "That triple for-loop is rather expensive, and it's much more efficient to collect the Jacobian in a\n", "single chunk of memory if possible.\n", "So even if you are using an iterative linear solver, such as [ScipyKrylov](../features/building_blocks/solvers/scipy_iter_solver.ipynb) or [PetscKrylov](../features/building_blocks/solvers/petsc_krylov.ipynb), an assembled Jacobian is generally more efficient.\n", "\n", "```{note}\n", "If you want to see how to add an assembled Jacobian to your model, check out [this feature doc](../features/core_features/working_with_derivatives/assembled_jacobian.ipynb).\n", "```\n", "\n", "## Sparse Assembled Jacobian\n", "\n", "In the majority of cases, if an assembled Jacobian is appropriate for your model, then you want to\n", "use a sparse assembled jacobian, which only allocates memory for the nonzero partial derivatives. You do this by settting `options['assembled_jac_type'] = 'sparse'` for that system. The specific sparse format for that sparse assembled jacobian will be determined by the linear solver for that system. For example, the `DirectSolver` prefers the `csc` format, while `ScipyKrylov` prefers `csr`. How does OpenMDAO know which partials are nonzero?\n", "The authors of the components in your model declared them using either a\n", "[dense](../features/core_features/working_with_derivatives/specifying_partials.ipynb) or form of `declare_partials`.\n", "\n", "Even if all of your components declared their partial derivatives as dense (or if they are all\n", "scalar variables and specifying sparsity doesn't have meaning), at the group level there is still a\n", "sparsity pattern to be taken advantage of.\n", "This sparsity arises from the way components are connected to one another, because unless there is\n", "a connection present, there is no need to allocate space associated with that portion of the Jacobian.\n", "We can see this clearly by looking at a collapsed form of the $N^2$ diagram with just the outputs shown.\n", "There are 7 scalar outputs, so we have a $7 \\times 7$ partial derivative Jacobian.\n", "Out of the possible 49 matrix entries, only 18 are actually nonzero. That makes it 63% sparse. Sellar is only a\n", "tiny toy problem, but in a real problem with thousands of variables, you will more commonly see sparsity percentages of over 90%." ] }, { "cell_type": "code", "execution_count": 2, "id": "8958cde1", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:09.676568Z", "iopub.status.busy": "2026-10-02T14:47:09.676392Z", "iopub.status.idle": "2026-10-02T14:47:11.851383Z", "shell.execute_reply": "2026-10-02T14:47:11.850728Z" }, "papermill": { "duration": 2.178142, "end_time": "2026-10-02T14:47:11.852467+00:00", "exception": false, "start_time": "2026-10-02T14:47:09.674325+00:00", "status": "completed" }, "tags": [ "remove-input", "remove-output" ] }, "outputs": [ { "data": { "application/papermill.record/text/html": "
class SellarDis1(om.ExplicitComponent):\n    """\n    Component containing Discipline 1 -- no derivatives version.\n    """\n\n    def __init__(self, units=None, scaling=None):\n        super().__init__()\n        self.execution_count = 0\n        self._units = units\n        self._do_scaling = scaling\n\n    def setup(self):\n\n        if self._units:\n            units = 'ft'\n        else:\n            units = None\n\n        if self._do_scaling:\n            ref = .1\n        else:\n            ref = 1.\n\n        # Global Design Variable\n        self.add_input('z', val=np.zeros(2), units=units)\n\n        # Local Design Variable\n        self.add_input('x', val=0., units=units)\n\n        # Coupling parameter\n        self.add_input('y2', val=1.0, units=units)\n\n        # Coupling output\n        self.add_output('y1', val=1.0, lower=0.1, upper=1000., units=units, ref=ref)\n\n    def setup_partials(self):\n        # Finite difference everything\n        self.declare_partials('*', '*', method='fd')\n\n    def compute(self, inputs, outputs):\n        """\n        Evaluates the equation\n        y1 = z1**2 + z2 + x1 - 0.2*y2\n        """\n\n        z1 = inputs['z'][0]\n        z2 = inputs['z'][1]\n        x1 = inputs['x']\n        y2 = inputs['y2']\n\n        outputs['y1'] = z1**2 + z2 + x1 - 0.2*y2\n\n        self.execution_count += 1\n
\n", "application/papermill.record/text/latex": "\\begin{Verbatim}[commandchars=\\\\\\{\\}]\n\\PY{k}{class}\\PY{+w}{ }\\PY{n+nc}{SellarDis1}\\PY{p}{(}\\PY{n}{om}\\PY{o}{.}\\PY{n}{ExplicitComponent}\\PY{p}{)}\\PY{p}{:}\n\\PY{+w}{ }\\PY{l+s+sd}{\\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\\PY{l+s+sd}{ Component containing Discipline 1 \\PYZhy{}\\PYZhy{} no derivatives version.}\n\\PY{l+s+sd}{ \\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf+fm}{\\PYZus{}\\PYZus{}init\\PYZus{}\\PYZus{}}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{k+kc}{None}\\PY{p}{,} \\PY{n}{scaling}\\PY{o}{=}\\PY{k+kc}{None}\\PY{p}{)}\\PY{p}{:}\n \\PY{n+nb}{super}\\PY{p}{(}\\PY{p}{)}\\PY{o}{.}\\PY{n+nf+fm}{\\PYZus{}\\PYZus{}init\\PYZus{}\\PYZus{}}\\PY{p}{(}\\PY{p}{)}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{execution\\PYZus{}count} \\PY{o}{=} \\PY{l+m+mi}{0}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}units} \\PY{o}{=} \\PY{n}{units}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}do\\PYZus{}scaling} \\PY{o}{=} \\PY{n}{scaling}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf}{setup}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{)}\\PY{p}{:}\n\n \\PY{k}{if} \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}units}\\PY{p}{:}\n \\PY{n}{units} \\PY{o}{=} \\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{ft}\\PY{l+s+s1}{\\PYZsq{}}\n \\PY{k}{else}\\PY{p}{:}\n \\PY{n}{units} \\PY{o}{=} \\PY{k+kc}{None}\n\n \\PY{k}{if} \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}do\\PYZus{}scaling}\\PY{p}{:}\n \\PY{n}{ref} \\PY{o}{=} \\PY{l+m+mf}{.1}\n \\PY{k}{else}\\PY{p}{:}\n \\PY{n}{ref} \\PY{o}{=} \\PY{l+m+mf}{1.}\n\n \\PY{c+c1}{\\PYZsh{} Global Design Variable}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}input}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{z}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{n}{np}\\PY{o}{.}\\PY{n}{zeros}\\PY{p}{(}\\PY{l+m+mi}{2}\\PY{p}{)}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{)}\n\n \\PY{c+c1}{\\PYZsh{} Local Design Variable}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}input}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{x}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{l+m+mf}{0.}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{)}\n\n \\PY{c+c1}{\\PYZsh{} Coupling parameter}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}input}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y2}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{l+m+mf}{1.0}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{)}\n\n \\PY{c+c1}{\\PYZsh{} Coupling output}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}output}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y1}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{l+m+mf}{1.0}\\PY{p}{,} \\PY{n}{lower}\\PY{o}{=}\\PY{l+m+mf}{0.1}\\PY{p}{,} \\PY{n}{upper}\\PY{o}{=}\\PY{l+m+mf}{1000.}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{,} \\PY{n}{ref}\\PY{o}{=}\\PY{n}{ref}\\PY{p}{)}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf}{setup\\PYZus{}partials}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{)}\\PY{p}{:}\n \\PY{c+c1}{\\PYZsh{} Finite difference everything}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{declare\\PYZus{}partials}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{*}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{*}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{method}\\PY{o}{=}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{fd}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{)}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf}{compute}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{,} \\PY{n}{inputs}\\PY{p}{,} \\PY{n}{outputs}\\PY{p}{)}\\PY{p}{:}\n\\PY{+w}{ }\\PY{l+s+sd}{\\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\\PY{l+s+sd}{ Evaluates the equation}\n\\PY{l+s+sd}{ y1 = z1**2 + z2 + x1 \\PYZhy{} 0.2*y2}\n\\PY{l+s+sd}{ \\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\n \\PY{n}{z1} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{z}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\\PY{p}{[}\\PY{l+m+mi}{0}\\PY{p}{]}\n \\PY{n}{z2} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{z}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\\PY{p}{[}\\PY{l+m+mi}{1}\\PY{p}{]}\n \\PY{n}{x1} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{x}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\n \\PY{n}{y2} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y2}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\n\n \\PY{n}{outputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y1}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]} \\PY{o}{=} \\PY{n}{z1}\\PY{o}{*}\\PY{o}{*}\\PY{l+m+mi}{2} \\PY{o}{+} \\PY{n}{z2} \\PY{o}{+} \\PY{n}{x1} \\PY{o}{\\PYZhy{}} \\PY{l+m+mf}{0.2}\\PY{o}{*}\\PY{n}{y2}\n\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{execution\\PYZus{}count} \\PY{o}{+}\\PY{o}{=} \\PY{l+m+mi}{1}\n\\end{Verbatim}\n", "application/papermill.record/text/plain": "class SellarDis1(om.ExplicitComponent):\n \"\"\"\n Component containing Discipline 1 -- no derivatives version.\n \"\"\"\n\n def __init__(self, units=None, scaling=None):\n super().__init__()\n self.execution_count = 0\n self._units = units\n self._do_scaling = scaling\n\n def setup(self):\n\n if self._units:\n units = 'ft'\n else:\n units = None\n\n if self._do_scaling:\n ref = .1\n else:\n ref = 1.\n\n # Global Design Variable\n self.add_input('z', val=np.zeros(2), units=units)\n\n # Local Design Variable\n self.add_input('x', val=0., units=units)\n\n # Coupling parameter\n self.add_input('y2', val=1.0, units=units)\n\n # Coupling output\n self.add_output('y1', val=1.0, lower=0.1, upper=1000., units=units, ref=ref)\n\n def setup_partials(self):\n # Finite difference everything\n self.declare_partials('*', '*', method='fd')\n\n def compute(self, inputs, outputs):\n \"\"\"\n Evaluates the equation\n y1 = z1**2 + z2 + x1 - 0.2*y2\n \"\"\"\n\n z1 = inputs['z'][0]\n z2 = inputs['z'][1]\n x1 = inputs['x']\n y2 = inputs['y2']\n\n outputs['y1'] = z1**2 + z2 + x1 - 0.2*y2\n\n self.execution_count += 1" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "code_src94" } }, "output_type": "display_data" } ], "source": [ "from openmdao.utils.notebook_utils import get_code\n", "from myst_nb import glue\n", "glue(\"code_src94\", get_code(\"openmdao.test_suite.components.sellar.SellarDis1\"), display=False)" ] }, { "cell_type": "markdown", "id": "4dbf1f26", "metadata": { "papermill": { "duration": 0.114677, "end_time": "2026-10-02T14:47:11.969185+00:00", "exception": false, "start_time": "2026-10-02T14:47:11.854508+00:00", "status": "completed" }, "tags": [] }, "source": [ ":::{dropdown} `SellarDis1` class definition \n", "\n", "{glue:}`code_src94`\n", ":::" ] }, { "cell_type": "code", "execution_count": 3, "id": "b5236b4b", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:11.980739Z", "iopub.status.busy": "2026-10-02T14:47:11.980496Z", "iopub.status.idle": "2026-10-02T14:47:11.993227Z", "shell.execute_reply": "2026-10-02T14:47:11.992562Z" }, "papermill": { "duration": 0.022961, "end_time": "2026-10-02T14:47:11.994051+00:00", "exception": false, "start_time": "2026-10-02T14:47:11.971090+00:00", "status": "completed" }, "tags": [ "remove-input", "remove-output" ] }, "outputs": [ { "data": { "application/papermill.record/text/html": "
class SellarDis2(om.ExplicitComponent):\n    """\n    Component containing Discipline 2 -- no derivatives version.\n    """\n\n    def __init__(self, units=None, scaling=None):\n        super().__init__()\n        self.execution_count = 0\n        self._units = units\n        self._do_scaling = scaling\n\n    def setup(self):\n        if self._units:\n            units = 'inch'\n        else:\n            units = None\n\n        if self._do_scaling:\n            ref = .18\n        else:\n            ref = 1.\n\n        # Global Design Variable\n        self.add_input('z', val=np.zeros(2), units=units)\n\n        # Coupling parameter\n        self.add_input('y1', val=1.0, units=units)\n\n        # Coupling output\n        self.add_output('y2', val=1.0, lower=0.1, upper=1000., units=units, ref=ref)\n\n    def setup_partials(self):\n        # Finite difference everything\n        self.declare_partials('*', '*', method='fd')\n\n    def compute(self, inputs, outputs):\n        """\n        Evaluates the equation\n        y2 = y1**(.5) + z1 + z2\n        """\n\n        z1 = inputs['z'][0]\n        z2 = inputs['z'][1]\n        y1 = inputs['y1']\n\n        # Note: this may cause some issues. However, y1 is constrained to be\n        # above 3.16, so lets just let it converge, and the optimizer will\n        # throw it out\n        if y1.real < 0.0:\n            y1 *= -1\n\n        outputs['y2'] = y1**.5 + z1 + z2\n\n        self.execution_count += 1\n
\n", "application/papermill.record/text/latex": "\\begin{Verbatim}[commandchars=\\\\\\{\\}]\n\\PY{k}{class}\\PY{+w}{ }\\PY{n+nc}{SellarDis2}\\PY{p}{(}\\PY{n}{om}\\PY{o}{.}\\PY{n}{ExplicitComponent}\\PY{p}{)}\\PY{p}{:}\n\\PY{+w}{ }\\PY{l+s+sd}{\\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\\PY{l+s+sd}{ Component containing Discipline 2 \\PYZhy{}\\PYZhy{} no derivatives version.}\n\\PY{l+s+sd}{ \\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf+fm}{\\PYZus{}\\PYZus{}init\\PYZus{}\\PYZus{}}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{k+kc}{None}\\PY{p}{,} \\PY{n}{scaling}\\PY{o}{=}\\PY{k+kc}{None}\\PY{p}{)}\\PY{p}{:}\n \\PY{n+nb}{super}\\PY{p}{(}\\PY{p}{)}\\PY{o}{.}\\PY{n+nf+fm}{\\PYZus{}\\PYZus{}init\\PYZus{}\\PYZus{}}\\PY{p}{(}\\PY{p}{)}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{execution\\PYZus{}count} \\PY{o}{=} \\PY{l+m+mi}{0}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}units} \\PY{o}{=} \\PY{n}{units}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}do\\PYZus{}scaling} \\PY{o}{=} \\PY{n}{scaling}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf}{setup}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{)}\\PY{p}{:}\n \\PY{k}{if} \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}units}\\PY{p}{:}\n \\PY{n}{units} \\PY{o}{=} \\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{inch}\\PY{l+s+s1}{\\PYZsq{}}\n \\PY{k}{else}\\PY{p}{:}\n \\PY{n}{units} \\PY{o}{=} \\PY{k+kc}{None}\n\n \\PY{k}{if} \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{\\PYZus{}do\\PYZus{}scaling}\\PY{p}{:}\n \\PY{n}{ref} \\PY{o}{=} \\PY{l+m+mf}{.18}\n \\PY{k}{else}\\PY{p}{:}\n \\PY{n}{ref} \\PY{o}{=} \\PY{l+m+mf}{1.}\n\n \\PY{c+c1}{\\PYZsh{} Global Design Variable}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}input}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{z}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{n}{np}\\PY{o}{.}\\PY{n}{zeros}\\PY{p}{(}\\PY{l+m+mi}{2}\\PY{p}{)}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{)}\n\n \\PY{c+c1}{\\PYZsh{} Coupling parameter}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}input}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y1}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{l+m+mf}{1.0}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{)}\n\n \\PY{c+c1}{\\PYZsh{} Coupling output}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{add\\PYZus{}output}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y2}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{val}\\PY{o}{=}\\PY{l+m+mf}{1.0}\\PY{p}{,} \\PY{n}{lower}\\PY{o}{=}\\PY{l+m+mf}{0.1}\\PY{p}{,} \\PY{n}{upper}\\PY{o}{=}\\PY{l+m+mf}{1000.}\\PY{p}{,} \\PY{n}{units}\\PY{o}{=}\\PY{n}{units}\\PY{p}{,} \\PY{n}{ref}\\PY{o}{=}\\PY{n}{ref}\\PY{p}{)}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf}{setup\\PYZus{}partials}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{)}\\PY{p}{:}\n \\PY{c+c1}{\\PYZsh{} Finite difference everything}\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{declare\\PYZus{}partials}\\PY{p}{(}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{*}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{*}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{,} \\PY{n}{method}\\PY{o}{=}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{fd}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{)}\n\n \\PY{k}{def}\\PY{+w}{ }\\PY{n+nf}{compute}\\PY{p}{(}\\PY{n+nb+bp}{self}\\PY{p}{,} \\PY{n}{inputs}\\PY{p}{,} \\PY{n}{outputs}\\PY{p}{)}\\PY{p}{:}\n\\PY{+w}{ }\\PY{l+s+sd}{\\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\\PY{l+s+sd}{ Evaluates the equation}\n\\PY{l+s+sd}{ y2 = y1**(.5) + z1 + z2}\n\\PY{l+s+sd}{ \\PYZdq{}\\PYZdq{}\\PYZdq{}}\n\n \\PY{n}{z1} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{z}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\\PY{p}{[}\\PY{l+m+mi}{0}\\PY{p}{]}\n \\PY{n}{z2} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{z}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\\PY{p}{[}\\PY{l+m+mi}{1}\\PY{p}{]}\n \\PY{n}{y1} \\PY{o}{=} \\PY{n}{inputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y1}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]}\n\n \\PY{c+c1}{\\PYZsh{} Note: this may cause some issues. However, y1 is constrained to be}\n \\PY{c+c1}{\\PYZsh{} above 3.16, so lets just let it converge, and the optimizer will}\n \\PY{c+c1}{\\PYZsh{} throw it out}\n \\PY{k}{if} \\PY{n}{y1}\\PY{o}{.}\\PY{n}{real} \\PY{o}{\\PYZlt{}} \\PY{l+m+mf}{0.0}\\PY{p}{:}\n \\PY{n}{y1} \\PY{o}{*}\\PY{o}{=} \\PY{o}{\\PYZhy{}}\\PY{l+m+mi}{1}\n\n \\PY{n}{outputs}\\PY{p}{[}\\PY{l+s+s1}{\\PYZsq{}}\\PY{l+s+s1}{y2}\\PY{l+s+s1}{\\PYZsq{}}\\PY{p}{]} \\PY{o}{=} \\PY{n}{y1}\\PY{o}{*}\\PY{o}{*}\\PY{l+m+mf}{.5} \\PY{o}{+} \\PY{n}{z1} \\PY{o}{+} \\PY{n}{z2}\n\n \\PY{n+nb+bp}{self}\\PY{o}{.}\\PY{n}{execution\\PYZus{}count} \\PY{o}{+}\\PY{o}{=} \\PY{l+m+mi}{1}\n\\end{Verbatim}\n", "application/papermill.record/text/plain": "class SellarDis2(om.ExplicitComponent):\n \"\"\"\n Component containing Discipline 2 -- no derivatives version.\n \"\"\"\n\n def __init__(self, units=None, scaling=None):\n super().__init__()\n self.execution_count = 0\n self._units = units\n self._do_scaling = scaling\n\n def setup(self):\n if self._units:\n units = 'inch'\n else:\n units = None\n\n if self._do_scaling:\n ref = .18\n else:\n ref = 1.\n\n # Global Design Variable\n self.add_input('z', val=np.zeros(2), units=units)\n\n # Coupling parameter\n self.add_input('y1', val=1.0, units=units)\n\n # Coupling output\n self.add_output('y2', val=1.0, lower=0.1, upper=1000., units=units, ref=ref)\n\n def setup_partials(self):\n # Finite difference everything\n self.declare_partials('*', '*', method='fd')\n\n def compute(self, inputs, outputs):\n \"\"\"\n Evaluates the equation\n y2 = y1**(.5) + z1 + z2\n \"\"\"\n\n z1 = inputs['z'][0]\n z2 = inputs['z'][1]\n y1 = inputs['y1']\n\n # Note: this may cause some issues. However, y1 is constrained to be\n # above 3.16, so lets just let it converge, and the optimizer will\n # throw it out\n if y1.real < 0.0:\n y1 *= -1\n\n outputs['y2'] = y1**.5 + z1 + z2\n\n self.execution_count += 1" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "code_src95" } }, "output_type": "display_data" } ], "source": [ "from openmdao.utils.notebook_utils import get_code\n", "from myst_nb import glue\n", "glue(\"code_src95\", get_code(\"openmdao.test_suite.components.sellar.SellarDis2\"), display=False)" ] }, { "cell_type": "markdown", "id": "a9d9f89a", "metadata": { "papermill": { "duration": 0.001721, "end_time": "2026-10-02T14:47:11.997580+00:00", "exception": false, "start_time": "2026-10-02T14:47:11.995859+00:00", "status": "completed" }, "tags": [] }, "source": [ ":::{dropdown} `SellarDis2` class definition \n", "\n", "{glue:}`code_src95`\n", ":::" ] }, { "cell_type": "code", "execution_count": 4, "id": "13ed919f", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:12.002416Z", "iopub.status.busy": "2026-10-02T14:47:12.002224Z", "iopub.status.idle": "2026-10-02T14:47:13.294832Z", "shell.execute_reply": "2026-10-02T14:47:13.294140Z" }, "papermill": { "duration": 1.296052, "end_time": "2026-10-02T14:47:13.295393+00:00", "exception": false, "start_time": "2026-10-02T14:47:11.999341+00:00", "status": "completed" }, "tags": [ "remove-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1790952433.242964] [runnervm8df0l:11128:0] ib_iface.c:1269 UCX ERROR mana_0: iface 0x55db535d7290 failed to create UD QP TX wr:256 sge:6 inl:64 resp:0 RX wr:4096 sge:1 resp:0 failed: Operation not supported\n", "[1790952433.243216] [runnervm8df0l:11128:0] ucp_worker.c:1412 UCX ERROR uct_iface_open(ud_verbs/mana_0:1) failed: Input/output error\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "[runnervm8df0l:11128] pml_ucx.c:313 Error: Failed to create UCP worker\n" ] }, { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import openmdao.api as om\n", "\n", "from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2\n", "\n", "class SellarMDAConnect(om.Group):\n", "\n", " def setup(self):\n", " cycle = self.add_subsystem('cycle', om.Group(), promotes_inputs=['x', 'z'])\n", " cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z'])\n", " cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z'])\n", " cycle.connect('d1.y1', 'd2.y1')\n", "\n", " ######################################\n", " # This is a \"forgotten\" connection!!\n", " ######################################\n", " #cycle.connect('d2.y2', 'd1.y2')\n", "\n", " self.set_input_defaults('x', 1.0)\n", " self.set_input_defaults('z', np.array([5.0, 2.0]))\n", "\n", " # Nonlinear Block Gauss Seidel is a gradient free solver\n", " cycle.nonlinear_solver = om.NonlinearBlockGS()\n", "\n", " self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',\n", " z=np.array([0.0, 0.0]), x=0.0),\n", " promotes_inputs=['x', 'z'])\n", "\n", " self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'))\n", " self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'))\n", "\n", " self.connect('cycle.d1.y1', ['obj_cmp.y1', 'con_cmp1.y1'])\n", " self.connect('cycle.d2.y2', ['obj_cmp.y2', 'con_cmp2.y2'])\n", "\n", "\n", "prob = om.Problem()\n", "\n", "prob.model = SellarMDAConnect()\n", "\n", "prob.driver = om.ScipyOptimizeDriver()\n", "prob.driver.options['optimizer'] = 'SLSQP'\n", "prob.driver.options['tol'] = 1e-8\n", "\n", "prob.set_solver_print(level=0)\n", "\n", "prob.model.add_design_var('x', lower=0, upper=10)\n", "prob.model.add_design_var('z', lower=0, upper=10)\n", "prob.model.add_objective('obj_cmp.obj')\n", "prob.model.add_constraint('con_cmp1.con1', upper=0)\n", "prob.model.add_constraint('con_cmp2.con2', upper=0)\n", "\n", "prob.setup()\n", "\n", "prob.set_val('x', 2.0)\n", "prob.set_val('z', [-1., -1.])\n", "\n", "om.n2(prob)" ] }, { "cell_type": "markdown", "id": "7f7efd4a", "metadata": { "papermill": { "duration": 0.001445, "end_time": "2026-10-02T14:47:13.298591+00:00", "exception": false, "start_time": "2026-10-02T14:47:13.297146+00:00", "status": "completed" }, "tags": [] }, "source": [ "
The partial-derivative Jacobian for the Sellar problem has only 18 nonzero values in it. Only 37% of the matrix elements are nonzero.
\n", "
\n", "\n", "If you chose to use the [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb), then it will use scipy's sparse [splu](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.splu.html) method to solve linear system for total derivatives.\n", "\n", "\n", "## Dense Assembled Jacobian\n", "\n", "Setting the `assembled_jac_type` option to 'dense' allocates a dense $n \\times n$ matrix, where $n$ is the sum of the sizes of all output variables in your model, to store partial derivatives in. So if you had a model that had 3 outputs of length 1000 each, then $n=3000$ and a $3000 \\times 3000$ matrix would be allocatec.\n", "\n", "Then whenever the Jacobian is needed, this dense matrix is provided.\n", "If you chose to use the [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb), then it will use scipy's [lu_factor](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html) and [lu_solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_solve.html#scipy.linalg.lu_solve)\n", "methods to solve linear system for total derivatives.\n", "\n", "If you have a very heavily-interconnected model, where many components connect to many others, then setting the `assembled_jac_type` option to 'dense' makes sense.\n", "However, the reality is that most models are very sparse at the group level even if the individual\n", "sub-Jacobians of the components are quite dense.\n", "So while there are some problems where it is appropriate, in most situations you should set `assembled_jac_type` to 'sparse' instead of 'dense'.\n", "\n", "\n", "## Matrix-Free Problems\n", "\n", "OpenMDAO is capable of solving linear systems in a matrix-free manner, to support situations where\n", "the Jacobian is too big to be fit into memory or when it's just too inefficient to do so.\n", "\n", "Practically speaking, if any components in your model use the [compute_jacvec_product](../features/core_features/working_with_components/explicit_component.ipynb) or [apply_linear](../features/core_features/working_with_components/implicit_component.ipynb) to provide derivatives, then you should be using a matrix-free linear solver architecture. These two methods provide linear operators that take in a vector and output the effect of multiplying it by a matrix. However, the underlying implementation does not actually need to assemble any matrices.\n", "\n", "Some high-fidelity PDE solvers will provide this kind of interface to get access to their partial derivatives.\n", "This kind of linear operator is also what is generally provided by algorithmic differentiations packages.\n", "\n", "Essentially, when you have problems with components that have very large array outputs (i.e. array\n", "lengths in the millions) and which run distributed across many cores, then a matrix-free linear\n", "solver architecture is something you want to consider.\n", "\n", "\n", "## Using the Model Hierarchy to Customize the Linear Solver Structure\n", "\n", "In OpenMDAO, your model is constructed via collections of Groups and Components arranged hierarchically.\n", "One of the main purposes of the hierarchy is to provide a means of sub-dividing a large and complex model into parts that can be solved using different methods.\n", "This creates a hierarchical solver architecture that is potentially both more efficient and more effective.\n", "The hierarchical solver architecture can be used for both nonlinear and linear solvers, but this section focuses specifically on the linear solver.\n", "\n", "## A Very Simple Example\n", "\n", "Consider, as an example, the [Sellar Problem](../basic_user_guide/multidisciplinary_optimization/sellar.ipynb) from the [Multidisciplinary Optimization User Guide](../basic_user_guide/basic_user_guide.ipynb).\n", "In that problem, coupling is created by a cyclic connection between the `d1` and `d2` components.\n", "You can see that coupling clearly in the n2 diagram below, because there are off-diagonal terms both above and below the diagonal inside the `cycle` group." ] }, { "cell_type": "code", "execution_count": 5, "id": "b75f329f", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:13.421945Z", "iopub.status.busy": "2026-10-02T14:47:13.421606Z", "iopub.status.idle": "2026-10-02T14:47:13.427051Z", "shell.execute_reply": "2026-10-02T14:47:13.426207Z" }, "papermill": { "duration": 0.12784, "end_time": "2026-10-02T14:47:13.427856+00:00", "exception": false, "start_time": "2026-10-02T14:47:13.300016+00:00", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "class SellarMDALinearSolver(om.Group):\n", " \"\"\"\n", " Group containing the Sellar MDA.\n", " \"\"\"\n", "\n", " def setup(self):\n", "\n", " cycle = self.add_subsystem('cycle', om.Group(), promotes=['*'])\n", " cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],\n", " promotes_outputs=['y1'])\n", " cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'],\n", " promotes_outputs=['y2'])\n", "\n", " self.set_input_defaults('x', 1.0)\n", " self.set_input_defaults('z', np.array([5.0, 2.0]))\n", "\n", " cycle.nonlinear_solver = om.NonlinearBlockGS()\n", " cycle.linear_solver = om.DirectSolver()\n", " cycle.options['assembled_jac_type'] = 'dense'\n", "\n", " self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',\n", " z=np.array([0.0, 0.0]), x=0.0),\n", " promotes=['x', 'z', 'y1', 'y2', 'obj'])\n", "\n", " self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])\n", " self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])" ] }, { "cell_type": "code", "execution_count": 6, "id": "86d80da5", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:13.433266Z", "iopub.status.busy": "2026-10-02T14:47:13.433059Z", "iopub.status.idle": "2026-10-02T14:47:13.489369Z", "shell.execute_reply": "2026-10-02T14:47:13.488529Z" }, "papermill": { "duration": 0.059665, "end_time": "2026-10-02T14:47:13.489884+00:00", "exception": false, "start_time": "2026-10-02T14:47:13.430219+00:00", "status": "completed" }, "tags": [ "remove-input" ] }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "\n", "from openmdao.test_suite.components.sellar import SellarDis1, SellarDis2\n", "\n", "class SellarMDAConnect(om.Group):\n", "\n", " def setup(self):\n", " cycle = self.add_subsystem('cycle', om.Group(), promotes_inputs=['x', 'z'])\n", " cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z'])\n", " cycle.add_subsystem('d2', SellarDis2(), promotes_inputs=['z'])\n", " cycle.connect('d1.y1', 'd2.y1')\n", "\n", " ######################################\n", " # This is a \"forgotten\" connection!!\n", " ######################################\n", " #cycle.connect('d2.y2', 'd1.y2')\n", "\n", " self.set_input_defaults('x', 1.0)\n", " self.set_input_defaults('z', np.array([5.0, 2.0]))\n", "\n", " # Nonlinear Block Gauss Seidel is a gradient free solver\n", " cycle.nonlinear_solver = om.NonlinearBlockGS()\n", "\n", " self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',\n", " z=np.array([0.0, 0.0]), x=0.0),\n", " promotes_inputs=['x', 'z'])\n", "\n", " self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'))\n", " self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'))\n", "\n", " self.connect('cycle.d1.y1', ['obj_cmp.y1', 'con_cmp1.y1'])\n", " self.connect('cycle.d2.y2', ['obj_cmp.y2', 'con_cmp2.y2'])\n", "\n", "\n", "prob = om.Problem()\n", "\n", "prob.model = SellarMDAConnect()\n", "\n", "prob.driver = om.ScipyOptimizeDriver()\n", "prob.driver.options['optimizer'] = 'SLSQP'\n", "prob.driver.options['tol'] = 1e-8\n", "\n", "prob.set_solver_print(level=0)\n", "\n", "prob.model.add_design_var('x', lower=0, upper=10)\n", "prob.model.add_design_var('z', lower=0, upper=10)\n", "prob.model.add_objective('obj_cmp.obj')\n", "prob.model.add_constraint('con_cmp1.con1', upper=0)\n", "prob.model.add_constraint('con_cmp2.con2', upper=0)\n", "\n", "prob.setup()\n", "prob.final_setup()\n", "\n", "prob.set_val('x', 2.0)\n", "prob.set_val('z', [-1., -1.])\n", "\n", "om.n2(prob)" ] }, { "cell_type": "markdown", "id": "5e96f864", "metadata": { "papermill": { "duration": 0.001754, "end_time": "2026-10-02T14:47:13.493227+00:00", "exception": false, "start_time": "2026-10-02T14:47:13.491473+00:00", "status": "completed" }, "tags": [] }, "source": [ "Since there is coupling in this model, there must also be some linear solver there to deal with it.\n", "One option would be to assign the [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb) right at the top level of the model, and have it compute an inverse of the full Jacobian.\n", "While that would certainly work, you're taking an inverse of a larger matrix than you really need to.\n", "\n", "Instead, as we've shown in the code above, you can assign the [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb) at the `cycle` level instead.\n", "The top level of the hierarchy will then be left with the default [LinearRunOnce](../features/building_blocks/solvers/linear_runonce.ipynb) solver in it.\n", "Effectively, the direct solver is being used to compute the coupled semi-total derivatives across the `cycle` group,\n", "which then makes the top level of the model have a feed-forward data path that can be solved with forward or back substitution\n", "(depending whether you select `fwd` or `rev` mode).\n", "\n", "To illustrate that visually, you can *right-click* on the cycle group in the n2 diagram above.\n", "This will collapse the cycle group to a single box, and you will see the resulting uncoupled, upper-triangular matrix structure that results.\n", "\n", "Practically speaking, for a tiny problem like [Sellar](../basic_user_guide/multidisciplinary_optimization/sellar.ipynb) there won't be any performance difference between putting\n", "the [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb) at the top, versus down in the `cycle` group. However, in larger models with hundreds or\n", "thousands of variables, the effect can be much more pronounced (e.g. if you're trying to invert a dense 10000x10000 matrix when\n", "you could be handling only a 10x10).\n", "\n", "More importantly, if you have models with high-fidelity codes like CFD or FEA in the hierarchy,\n", "you simply may not be able to use a [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb) at the top of the model, but there may still be a\n", "portion of the model where it makes sense. As you can see, understanding how to take advantage of the model hierarchy in\n", "order to customize the linear solver behavior becomes more important as your model complexity increases.\n", "\n", "\n", "## A More Realistic Example\n", "\n", "Consider an aerostructural model of an aircraft wing comprised of a Computational Fluid Dynamics (CFD) solver, a simple\n", "finite-element beam analysis, with a fuel-burn objective and a $C_l$ constraint.\n", "\n", "In OpenMDAO the model is set up as follows:" ] }, { "cell_type": "code", "execution_count": 7, "id": "8976f8e7", "metadata": { "execution": { "iopub.execute_input": "2026-10-02T14:47:13.539420Z", "iopub.status.busy": "2026-10-02T14:47:13.539272Z", "iopub.status.idle": "2026-10-02T14:47:13.586885Z", "shell.execute_reply": "2026-10-02T14:47:13.585887Z" }, "papermill": { "duration": 0.05037, "end_time": "2026-10-02T14:47:13.587383+00:00", "exception": false, "start_time": "2026-10-02T14:47:13.537013+00:00", "status": "completed" }, "tags": [ "remove-input" ] }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = om.Problem()\n", "dvs = p.model.add_subsystem('design_vars', om.IndepVarComp(), promotes=['*'])\n", "dvs.add_output('x_aero')\n", "dvs.add_output('x_struct')\n", "aerostruct = p.model.add_subsystem('aerostruct_cycle', om.Group(), promotes=['*'])\n", "aerostruct.add_subsystem('aero',\n", " om.ExecComp(['w = u+x_aero', 'Cl=u+x_aero', 'Cd = u + x_aero']),\n", " promotes=['*'])\n", "aerostruct.add_subsystem('struct', om.ExecComp(['u = w+x_struct', 'mass=x_struct']),\n", " promotes=['*'])\n", "\n", "p.model.add_subsystem('objective', om.ExecComp('f=mass+Cl/Cd'), promotes=['*'])\n", "p.model.add_subsystem('constraint', om.ExecComp('g=Cl'), promotes=['*'])\n", "\n", "p.setup()\n", "\n", "om.n2(p, outfile='aerostruct_n2.html', embeddable=True, show_browser=False)" ] }, { "cell_type": "markdown", "id": "2ede8ca3", "metadata": { "papermill": { "duration": 0.002113, "end_time": "2026-10-02T14:47:13.651127+00:00", "exception": false, "start_time": "2026-10-02T14:47:13.649014+00:00", "status": "completed" }, "tags": [] }, "source": [ "Note that this model has almost the exact same structure in its $N^2$ diagram as the sellar problem.\n", "Specifically the coupling between the aerodynamics and structural analyses can be isolated from the rest of the model.\n", "Those two are grouped together in the `aerostruct_cycle` group, giving the top level of the model a feed-forward structure.\n", "There is a subtle difference though; the Sellar problem is constructed of all explicit components but this aerostructural problem has two implicit analyses in the `aero` and `struct` components.\n", "Practically speaking, the presence of a CFD component means that the model is too big to use a [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb) at the top level of its hierarchy.\n", "\n", "Instead, based on the advice in the [Theory Manual entry on selecting which kind of linear solver to use](../theory_manual/setup_linear_solvers.ipynb),\n", "the feed-forward structure on the top level indicates that the default [LinearRunOnce](../features/building_blocks/solvers/linear_runonce.ipynb) solver is a good choice for that level of the model.\n", "\n", "So now the challenge is to select a good linear solver architecture for the `cycle` group.\n", "One possible approach is to use the [LinearBlockGS](../features/building_blocks/solvers/linear_block_gs.ipynb) solver for the `cycle`,\n", "and then assign additional solvers to the aerodynamics and structural analyses.\n", "\n", "```{note}\n", "Choosing LinearBlockGaussSeidel is analogous to solving the nonlinear system with a NonLinearBlockGaussSeidel solver.\n", "\n", "Despite the analogy, it is not required nor even advised that your linear solver architecture match your nonlinear solver architecture. It could very well be a better choice to use the [PetscKrylov](../features/building_blocks/solvers/petsc_krylov.ipynb) solver for the `cycle` level, even if the [NonlinearBlockGS](../features/building_blocks/solvers/nonlinear_block_gs.ipynb) solver was set as the nonlinear solver.\n", "```\n", "\n", "The [LinearBlockGS](../features/building_blocks/solvers/linear_block_gs.ipynb) solver requires that any implicit components underneath it have their own linear\n", "solvers to converge their part of the overall linear system. So a [PetscKrylov](../features/building_blocks/solvers/petsc_krylov.ipynb) solver is used for `aero`\n", "and a [DirectSolver](../features/building_blocks/solvers/direct_solver.ipynb) is use for `struct`. Looking back at the figure above, notice that these solvers\n", "are all called out in their respective hierarchical locations." ] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.14" }, "papermill": { "default_parameters": {}, "duration": 5.384321, "end_time": "2026-10-02T14:47:14.369037+00:00", "environment_variables": {}, "exception": null, "input_path": "/home/runner/work/OpenMDAO/OpenMDAO/openmdao/docs/openmdao_book/theory_manual/setup_linear_solvers.ipynb", "output_path": "/home/runner/work/OpenMDAO/OpenMDAO/openmdao/docs/_executed_book/theory_manual/setup_linear_solvers.ipynb", "parameters": {}, "start_time": "2026-10-02T14:47:08.984716+00:00", "version": "2.7.0" } }, "nbformat": 4, "nbformat_minor": 5 }