|
| 1 | +The OptionsManager and the AppContext |
| 2 | +------------------------------------- |
| 3 | + |
| 4 | +The PETSc options provide a simple but powerful DSL for configuring composable solvers. |
| 5 | +However, their main limitation is that the values of each option is limited to primitive C types, e.g. ``str``, ``float``, ``int``, or ``complex``. |
| 6 | +Sometimes more advanced data is useful or essential for building a particular solver. |
| 7 | + |
| 8 | +:class:`petsctools.AppContext <.appctx.AppContext>` fulfils this need by providing a means of passing arbitrary Python types through to Python PETSc types (e.g. Python type PCs). |
| 9 | +In this demo we show how to use the :class:`~.appctx.AppContext` to pass data to a custom Python type PC using the variable coefficient diffusion equation as an example. |
| 10 | + |
| 11 | + |
| 12 | +Diffusion equation with variable coefficients |
| 13 | +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 14 | + |
| 15 | +The diffusion equation with coefficient :math:`\sigma(x)` depending on the spatial coordinate is: |
| 16 | + |
| 17 | +.. math:: |
| 18 | +
|
| 19 | + u - \nabla\cdot\left(\sigma(x)\nabla u\right) = b |
| 20 | +
|
| 21 | +We will solve this matrix with finite differences with the standard 3 point central stencil. |
| 22 | +The particular details of the discretisation are not essential for this demo so we will be brief in the description. |
| 23 | + |
| 24 | +If :math:`D` is the assembled matrix for the finite difference gradient stencil, and :math:`\Sigma` is a diagonal matrix with the value of the diffusion coefficient at each grid point, then the assembled matrix for the diffusion equation is: |
| 25 | + |
| 26 | +.. math:: |
| 27 | +
|
| 28 | + Au = \left(I + D^{T}\Sigma D\right)u = b, |
| 29 | + \quad |
| 30 | + \Sigma_{ii} = \sigma(x_{i}) |
| 31 | +
|
| 32 | +The following Python function takes a numpy array ``sigma`` with the value of :math:`\sigma` at each grid point and assembles a sparse (``aij``) PETSc Mat for the diffusion equation. |
| 33 | +We will use it later to build the :class:`~petsc4py.PETSc.Mat` for a :class:`~petsc4py.PETSc.KSP` to solve the diffusion equation. |
| 34 | + |
| 35 | +.. literalinclude:: ../../tests/docs/test_appctx_docs.py |
| 36 | + :language: python3 |
| 37 | + :dedent: |
| 38 | + :start-after: [appctx_docs create_mat-start] |
| 39 | + :end-before: [appctx_docs create_mat-end] |
| 40 | + |
| 41 | +A PC needing Python data |
| 42 | +~~~~~~~~~~~~~~~~~~~~~~~~ |
| 43 | + |
| 44 | +Imagine a scenario where we need to solve :math:`A` multiple times and :math:`\sigma` changes slightly each time, for example if we are solving the unsteady diffusion equation with time-varying coefficients. |
| 45 | +Rather than recomputing a preconditioner every time :math:`A` changes, we might instead find a representative :math:`\sigma_{p}` and use that to compute a preconditioner which can be reused for all solves. |
| 46 | + |
| 47 | +For simplicity, in this demo the preconditioner :math:`P` will just be the diagonal of the assembled matrix :math:`A_{p}` for the diffusion equation with :math:`\sigma_{p}` with a simple scaling factor :math:`\omega`: |
| 48 | + |
| 49 | +.. math:: |
| 50 | +
|
| 51 | + A_{p} = I + D^{T}\Sigma_{p} D, |
| 52 | + \quad |
| 53 | + P = \omega^{-1}\mathrm{diag}(A_{p}). |
| 54 | +
|
| 55 | +The diagonal of a matrix is clearly not expensive to compute, but in practice we would use a factorisation of :math:`A_{p}` which would be more expensive to compute and so would be more worthwhile reusing. |
| 56 | + |
| 57 | +The preconditioner defined above is implemented with a Python type PC called ``DiffusionJacobiPC`` in the code below. |
| 58 | +Constructing :math:`P` requires two values, :math:`\sigma_{p}` and :math:`\omega`, which must be provided by the user. |
| 59 | + |
| 60 | +1. The scaling factor :math:`\omega` is just a real number, and can therefore be passed as usual via the :class:`PETSc.Options <petsc4py.PETSc.Options>` using the ``"djacobi_scale"`` option. |
| 61 | + |
| 62 | +2. The diffusion coefficient at each grid point :math:`\sigma_{p}(x_{i})` is defined as a numpy array. |
| 63 | + This is clearly not a primitive type and so cannot be passed via the :class:`PETSc.Options <petsc4py.PETSc.Options>` directly. |
| 64 | + Instead, we access it via the :class:`~.appctx.AppContext` using the ``"djacobi_sigma"`` key. |
| 65 | + |
| 66 | +The :class:`~petsctools.appctx.AppContext` mimics the :class:`PETSc.Options <petsc4py.PETSc.Options>` very closely, but can contain arbitrary Python data. |
| 67 | +We will see below how to add ``sigma`` into the :class:`~petsctools.appctx.AppContext` so that it is available to the ``DiffusionJacobiPC``. |
| 68 | + |
| 69 | +.. literalinclude:: ../../tests/docs/test_appctx_docs.py |
| 70 | + :language: python3 |
| 71 | + :dedent: |
| 72 | + :start-after: [appctx_docs pc-start] |
| 73 | + :end-before: [appctx_docs pc-end] |
| 74 | + |
| 75 | +Assembling the system |
| 76 | +~~~~~~~~~~~~~~~~~~~~~ |
| 77 | + |
| 78 | +We specify the diffusion coefficient as some random variations :math:`\sigma'` around a constant value :math:`\overline{\sigma}`, i.e. :math:`\sigma(x) = \overline{\sigma} + \sigma'(x)`. |
| 79 | +Assuming that :math:`\sigma'` is the component that may vary from solve to solve, we use :math:`\sigma_{p}=\overline{\sigma}`. |
| 80 | + |
| 81 | +.. literalinclude:: ../../tests/docs/test_appctx_docs.py |
| 82 | + :language: python3 |
| 83 | + :dedent: |
| 84 | + :start-after: [appctx_docs create_ksp-start] |
| 85 | + :end-before: [appctx_docs create_ksp-end] |
| 86 | + |
| 87 | +The Options and the AppContext |
| 88 | +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 89 | + |
| 90 | +Now we configure ``ksp`` by passing PETSc options as key-value pairs in the ``parameters`` dictionary to :func:`petsctools.set_from_options <.options.set_from_options>`. |
| 91 | +This function will create a :class:`petsctools.OptionsManager <.options.OptionsManager>` and attach it to ``ksp``. |
| 92 | + |
| 93 | +We can see that common options, e.g. ``"ksp_type"`` are set as usual in the ``parameters`` dictionary. |
| 94 | +However, when we come to the ``"djacobi_sigma"`` value we use the :class:`petsctools.AppContextManager <.appctx.AppContextManager>` class. |
| 95 | +The :meth:`AppContextManager.add <.appctx.AppContextManager.add>` method returns a unique value which is used to associate a PETSc option to whatever data was passed to ``add``. |
| 96 | +For example, adding the key-value pair ``"djacobi_sigma": appmngr.add(sigma_p)`` to the ``parameters`` dictionary means that, during the solve, we will be able to access ``sigma_p`` via ``AppContext()["djacobi_sigma"]``. |
| 97 | + |
| 98 | +We then pass the :class:`~.appctx.AppContextManager` to :func:`~.options.set_from_options` so that it can be attached to ``ksp`` and its data can be made available later on during the solve. |
| 99 | + |
| 100 | +.. literalinclude:: ../../tests/docs/test_appctx_docs.py |
| 101 | + :language: python3 |
| 102 | + :dedent: |
| 103 | + :start-after: [appctx_docs set_from_options-start] |
| 104 | + :end-before: [appctx_docs set_from_options-end] |
| 105 | + |
| 106 | +Solving the KSP |
| 107 | +~~~~~~~~~~~~~~~ |
| 108 | + |
| 109 | +Now we come to actually solving the linear equation :math:`Au=b`. |
| 110 | +To avoid memory leaks, :func:`~.options.set_from_options` does not permanently insert the contents of ``parameters`` and the ``appmngr`` into the global :class:`PETSc.Options <petsc4py.PETSc.Options>` and :class:`~.appctx.AppContext` databases respectively. |
| 111 | +Instead, we use the :func:`petsctools.inserted_options <.options.inserted_options>` context manager. |
| 112 | +On entry, this context manager inserts the contents of ``parameters`` and ``appmngr`` into the global databases, and on exit it removes them again. |
| 113 | +This means that we need to use the :func:`~.options.inserted_options` context manager whenever these entries will be needed, for example during the solve when the KSP and PC are being set up. |
| 114 | + |
| 115 | +.. literalinclude:: ../../tests/docs/test_appctx_docs.py |
| 116 | + :language: python3 |
| 117 | + :dedent: |
| 118 | + :start-after: [appctx_docs solve-start] |
| 119 | + :end-before: [appctx_docs solve-end] |
0 commit comments