\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", "A [DenseJacobian](../_srcdocs/packages/jacobians/assembled_jacobian.ipynb) 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\n", "[DenseJacobian](../_srcdocs/packages/jacobians/assembled_jacobian.ipynb) would allocate a $3000 \\times 3000$ matrix.\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 a\n", "[DenseJacobian](../_srcdocs/packages/jacobians/assembled_jacobian.ipynb) 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 don't want to use\n", "the [DenseJacobian](../_srcdocs/packages/jacobians/assembled_jacobian.ipynb).\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.md).\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": null, "metadata": {}, "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", " d1 = cycle.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],\n", " promotes_outputs=['y1'])\n", " d2 = 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", "\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": null, "metadata": { "tags": [ "remove-input" ] }, "outputs": [], "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", " cycle.set_input_defaults('x', 1.0)\n", " cycle.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", "metadata": {}, "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": null, "metadata": { "tags": [ "remove-input" ] }, "outputs": [], "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", "metadata": {}, "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.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }