From 794d48a09481fa79c9095ceacf302dca2501ece6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 20 Mar 2018 20:17:38 -0500 Subject: [PATCH 1/4] Doc improvements, more docs to code files --- doc/tools.rst | 40 ++-------------------------------------- modepy/modal_decay.py | 4 ++++ 2 files changed, 6 insertions(+), 38 deletions(-) diff --git a/doc/tools.rst b/doc/tools.rst index 1d64d3f..b5b0006 100644 --- a/doc/tools.rst +++ b/doc/tools.rst @@ -8,7 +8,7 @@ Version query .. data:: VERSION - A tuple like *(2013,1,1)* indicating modepy's version. + A tuple like *(2013, 1, 1)* indicating modepy's version. .. data:: VERSION_TEXT @@ -17,49 +17,13 @@ Version query Interpolation matrices ---------------------- -.. currentmodule:: modepy - -.. autofunction:: vandermonde - -Vandermonde matrices are very useful to concisely express interpolation. For -instance, one may use the inverse :math:`V^{-1}` of a Vandermonde matrix -:math:`V` to map nodal values (i.e. point values at the nodes corresponding to -:math:`V`) into modal (i.e. basis) coefficients. :math:`V` itself maps modal -coefficients to nodal values. This allows easy resampling: - -.. autofunction:: resampling_matrix - -Vandermonde matrices also obey the following relationship useful for obtaining -point interpolants: - -.. math:: - - V^T [\text{interpolation coefficents to point $x$}] = \phi_i(x), - -where :math:`(\phi_i)_i` is the basis of functions underlying :math:`V`. - -.. autofunction:: inverse_mass_matrix - -.. autofunction:: mass_matrix - -.. autofunction:: modal_face_mass_matrix - -.. autofunction:: nodal_face_mass_matrix - -Differentiation is also convenient to express by using :math:`V^{-1}` to -obtain modal values and then using a Vandermonde matrix for the derivatives -of the basis to return to nodal values. - -.. autofunction:: differentiation_matrices +.. automodule:: modepy.matrices Modal decay/residual -------------------- .. automodule:: modepy.modal_decay -.. autofunction:: fit_modal_decay -.. autofunction:: estimate_relative_expansion_residual - Interpolation quality --------------------- diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py index 9004123..a16208a 100644 --- a/modepy/modal_decay.py +++ b/modepy/modal_decay.py @@ -39,6 +39,10 @@ The method implemented in this module follows this article: http://arxiv.org/abs/1102.3190 .. versionadded:: 2013.2 + +.. autofunction:: simplex_interp_error_coefficient_estimator_matrix +.. autofunction:: fit_modal_decay +.. autofunction:: estimate_relative_expansion_residual """ -- GitLab From d8f3d59c3ecbbfeeb0d86f1f04125233870f206b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 20 Mar 2018 20:18:56 -0500 Subject: [PATCH 2/4] Second half of matrices doc move --- modepy/matrices.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/modepy/matrices.py b/modepy/matrices.py index 4e53001..34aaaa3 100644 --- a/modepy/matrices.py +++ b/modepy/matrices.py @@ -27,6 +27,44 @@ import numpy as np import numpy.linalg as la +__doc__ = r""" +.. currentmodule:: modepy + +.. autofunction:: vandermonde + +Vandermonde matrices are very useful to concisely express interpolation. For +instance, one may use the inverse :math:`V^{-1}` of a Vandermonde matrix +:math:`V` to map nodal values (i.e. point values at the nodes corresponding to +:math:`V`) into modal (i.e. basis) coefficients. :math:`V` itself maps modal +coefficients to nodal values. This allows easy resampling: + +.. autofunction:: resampling_matrix + +Vandermonde matrices also obey the following relationship useful for obtaining +point interpolants: + +.. math:: + + V^T [\text{interpolation coefficents to point $x$}] = \phi_i(x), + +where :math:`(\phi_i)_i` is the basis of functions underlying :math:`V`. + +.. autofunction:: inverse_mass_matrix + +.. autofunction:: mass_matrix + +.. autofunction:: modal_face_mass_matrix + +.. autofunction:: nodal_face_mass_matrix + +Differentiation is also convenient to express by using :math:`V^{-1}` to +obtain modal values and then using a Vandermonde matrix for the derivatives +of the basis to return to nodal values. + +.. autofunction:: differentiation_matrices +""" + + def vandermonde(functions, nodes): """Return a (generalized) Vandermonde matrix. -- GitLab From 87129fab3fbfc7d77fa560e25f7e199e8f4c3311 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 20 Mar 2018 20:19:15 -0500 Subject: [PATCH 3/4] Add versions of the simplex bases that return mode IDs --- modepy/modes.py | 85 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 67 insertions(+), 18 deletions(-) diff --git a/modepy/modes.py b/modepy/modes.py index e6dfcc1..f841629 100644 --- a/modepy/modes.py +++ b/modepy/modes.py @@ -61,10 +61,14 @@ Dimension-independent basis getters for simplices Scientific Computing 6, no. 4 (December 1, 1991): 345–390. http://dx.doi.org/10.1007/BF01060030 +.. autofunction:: simplex_onb_with_mode_ids + .. autofunction:: simplex_onb .. autofunction:: grad_simplex_onb +.. autofunction:: simplex_monomial_basis_with_mode_ids + .. autofunction:: simplex_monomial_basis .. autofunction:: grad_simplex_monomial_basis @@ -438,15 +442,17 @@ def grad_monomial(order, rst): # {{{ dimension-independent interface for simplices -def simplex_onb(dims, n): +def simplex_onb_with_mode_ids(dims, n): """Return a list of orthonormal basis functions in dimension *dims* of maximal total degree *n*. - :returns: a class:`tuple` of functions, each of which - accepts arrays of shape *(dims, npts)* - and return the function values as an array of size *npts*. - 'Scalar' evaluation, by passing just one vector of length *dims*, - is also supported. + :returns: a tuple ``(mode_ids, basis)``, where *basis* is a class:`tuple` + of functions, each of which accepts arrays of shape *(dims, npts)* and + return the function values as an array of size *npts*. 'Scalar' + evaluation, by passing just one vector of length *dims*, is also supported. + *mode_ids* is a tuple of the same length as *basis*, where each entry + is a tuple of integers describing the order of the mode along the coordinate + axes. See the following publications: @@ -454,9 +460,7 @@ def simplex_onb(dims, n): * |koornwinder-ref| * |dubiner-ref| - .. versionchanged:: 2013.2 - - Made return value a tuple, to make bases hashable. + ... versionadded: 2018.1 """ from functools import partial from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ @@ -469,18 +473,45 @@ def simplex_onb(dims, n): else: return np.ones(x.shape[1]) - return (zerod_basis,) + return ((0,),), (zerod_basis,) elif dims == 1: - return tuple(partial(jacobi, 0, 0, i) for i in range(n+1)) + mode_ids = tuple(range(n+1)) + return mode_ids, tuple(partial(jacobi, 0, 0, i) for i in mode_ids) elif dims == 2: - return tuple(partial(pkdo_2d, order) for order in gnitstam(n, dims)) + mode_ids = tuple(gnitstam(n, dims)) + return mode_ids, tuple(partial(pkdo_2d, order) for order in mode_ids) elif dims == 3: - return tuple(partial(pkdo_3d, order) for order in gnitstam(n, dims)) + mode_ids = tuple(gnitstam(n, dims)) + return mode_ids, tuple(partial(pkdo_3d, order) for order in mode_ids) else: raise NotImplementedError("%d-dimensional bases" % dims) +def simplex_onb(dims, n): + """Return a list of orthonormal basis functions in dimension *dims* of maximal + total degree *n*. + + :returns: a class:`tuple` of functions, each of which + accepts arrays of shape *(dims, npts)* + and return the function values as an array of size *npts*. + 'Scalar' evaluation, by passing just one vector of length *dims*, + is also supported. + + See the following publications: + + * |proriol-ref| + * |koornwinder-ref| + * |dubiner-ref| + + .. versionchanged:: 2013.2 + + Made return value a tuple, to make bases hashable. + """ + mode_ids, basis = simplex_onb_with_mode_ids(dims, n) + return basis + + def grad_simplex_onb(dims, n): """Return the gradients of the functions returned by :func:`simplex_onb`. @@ -515,6 +546,27 @@ def grad_simplex_onb(dims, n): raise NotImplementedError("%d-dimensional bases" % dims) +def simplex_monomial_basis_with_mode_ids(dims, n): + """Return a list of monomial basis functions in dimension *dims* of maximal + total degree *n*. + + :returns: a tuple ``(mode_ids, basis)``, where *basis* is a class:`tuple` + of functions, each of which accepts arrays of shape *(dims, npts)* and + return the function values as an array of size *npts*. 'Scalar' + evaluation, by passing just one vector of length *dims*, is also supported. + *mode_ids* is a tuple of the same length as *basis*, where each entry + is a tuple of integers describing the order of the mode along the coordinate + axes. + + .. versionadded:: 2018.1 + """ + from functools import partial + from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ + as gnitstam + mode_ids = tuple(gnitstam(n, dims)) + return mode_ids, tuple(partial(monomial, order) for order in mode_ids) + + def simplex_monomial_basis(dims, n): """Return a list of monomial basis functions in dimension *dims* of maximal total degree *n*. @@ -527,11 +579,8 @@ def simplex_monomial_basis(dims, n): .. versionadded:: 2016.1 """ - - from functools import partial - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam - return tuple(partial(monomial, order) for order in gnitstam(n, dims)) + mode_ids, basis = simplex_monomial_basis_with_mode_ids(dims, n) + return basis def grad_simplex_monomial_basis(dims, n): -- GitLab From 0d5419d4312620f89d6d284eb136d9739ab56a09 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 20 Mar 2018 20:19:33 -0500 Subject: [PATCH 4/4] Add simplex_interp_error_coefficient_estimator_matrix --- modepy/modal_decay.py | 35 ++++++++++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 3 deletions(-) diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py index a16208a..dee45fe 100644 --- a/modepy/modal_decay.py +++ b/modepy/modal_decay.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2010-2012 Andreas Kloeckner" @@ -24,6 +22,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range + import numpy as np import numpy.linalg as la @@ -46,6 +46,35 @@ The method implemented in this module follows this article: """ +def simplex_interp_error_coefficient_estimator_matrix( + unit_nodes, order, n_tail_orders): + """Return a matrix :math:`C` that, when multiplied by a vector of nodal values, + yields the coeffiicients belonging to the basis functions of the *n_tail_orders* + highest orders. + + The 2-norm of the resulting coefficents can be used as an estimate of the + interpolation error. + + .. versionadded:: 2018.1 + """ + + from modepy.matrices import vandermonde + from modepy.modes import simplex_onb_with_mode_ids + + dim, nunit_nodes = unit_nodes.shape + + mode_ids, basis = simplex_onb_with_mode_ids(dim, order) + vdm = vandermonde(basis, unit_nodes) + vdm_inv = la.inv(vdm) + + order_vector = np.array([sum(mode_id) for mode_id in mode_ids]) + + max_order = np.max(order_vector) + assert max_order == order + + return vdm_inv[order_vector > max_order-n_tail_orders] + + def make_mode_number_vector(mode_order_tuples, ignored_modes): node_cnt = len(mode_order_tuples) -- GitLab