From 9dca44da06adb6d61db80f10157118130f7b5ed3 Mon Sep 17 00:00:00 2001 From: "Thomas H. Gibson" Date: Thu, 13 Jan 2022 13:43:11 -0600 Subject: [PATCH] Fix metric computation for affine + overintegrated (#210) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add overintegration option for wave-op-mpi * Add OI+lazy to examples CI testing * Draft: Use P^0 discr on affine groups even if quad grid requested * Refactor geoderiv connection/interpolation steps * Fix quadrature interpolation * Add pytato pytest array context factory * Update metric unit tests * Clarifying comment Co-authored-by: Andreas Klöckner * Add nonaffine option for wave op example * Run lazy + nonaffine + overintegration in examples CI testing * Add quadrature case for metric unit tests * Add optional import for debugging outside of pytest * Add docstrings for _use_geoderiv_connection Co-authored-by: Andreas Kloeckner --- .github/workflows/ci.yml | 1 + examples/wave/wave-op-mpi.py | 88 +++++++++++---- grudge/array_context.py | 8 ++ grudge/geometry/metrics.py | 208 ++++++++++++++++++++++------------- test/test_grudge.py | 44 +------- test/test_metrics.py | 125 +++++++++++++++++++++ 6 files changed, 335 insertions(+), 139 deletions(-) create mode 100644 test/test_metrics.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 30224593..3cef42c2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -65,6 +65,7 @@ jobs: run_examples python wave/wave-op-mpi.py --lazy + python wave/wave-op-mpi.py --lazy --quad --nonaffine docs: name: Documentation diff --git a/examples/wave/wave-op-mpi.py b/examples/wave/wave-op-mpi.py index 63b052ed..3be4035b 100644 --- a/examples/wave/wave-op-mpi.py +++ b/examples/wave/wave-op-mpi.py @@ -45,6 +45,8 @@ from pytools.obj_array import flat_obj_array, make_obj_array from meshmode.dof_array import DOFArray from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.dof_desc import as_dofdesc, DOFDesc, DISCR_TAG_BASE, DISCR_TAG_QUAD +from grudge.trace_pair import TracePair from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer @@ -76,8 +78,9 @@ class WaveState: def wave_flux(dcoll, c, w_tpair): u = w_tpair.u v = w_tpair.v + dd = w_tpair.dd - normal = thaw(dcoll.normal(w_tpair.dd), u.int.array_context) + normal = thaw(dcoll.normal(dd), u.int.array_context) flux_weak = WaveState( u=v.avg @ normal, @@ -91,14 +94,29 @@ def wave_flux(dcoll, c, w_tpair): v=0.5 * v_jump * normal, ) - return op.project(dcoll, w_tpair.dd, "all_faces", c*flux_weak) + return op.project(dcoll, dd, dd.with_dtag("all_faces"), c*flux_weak) -def wave_operator(dcoll, c, w): - u = w.u - v = w.v +def wave_operator(dcoll, c, w, quad_tag=None): + dd_base = as_dofdesc("vol") + dd_vol = DOFDesc("vol", quad_tag) + dd_faces = DOFDesc("all_faces", quad_tag) + dd_btag = as_dofdesc(BTAG_ALL).with_discr_tag(quad_tag) - dir_w = op.project(dcoll, "vol", BTAG_ALL, w) + def interp_to_surf_quad(utpair): + local_dd = utpair.dd + local_dd_quad = local_dd.with_discr_tag(quad_tag) + return TracePair( + local_dd_quad, + interior=op.project(dcoll, local_dd, local_dd_quad, utpair.int), + exterior=op.project(dcoll, local_dd, local_dd_quad, utpair.ext) + ) + + w_quad = op.project(dcoll, dd_base, dd_vol, w) + u = w_quad.u + v = w_quad.v + + dir_w = op.project(dcoll, dd_base, dd_btag, w) dir_u = dir_w.u dir_v = dir_w.v dir_bval = WaveState(u=dir_u, v=dir_v) @@ -108,19 +126,20 @@ def wave_operator(dcoll, c, w): op.inverse_mass( dcoll, WaveState( - u=-c*op.weak_local_div(dcoll, v), - v=-c*op.weak_local_grad(dcoll, u) + u=-c*op.weak_local_div(dcoll, dd_vol, v), + v=-c*op.weak_local_grad(dcoll, dd_vol, u) ) + op.face_mass( dcoll, + dd_faces, wave_flux( dcoll, c=c, w_tpair=op.bdry_trace_pair(dcoll, - BTAG_ALL, + dd_btag, interior=dir_bval, exterior=dir_bc) ) + sum( - wave_flux(dcoll, c=c, w_tpair=tpair) + wave_flux(dcoll, c=c, w_tpair=interp_to_surf_quad(tpair)) for tpair in op.interior_trace_pairs(dcoll, w) ) ) @@ -164,7 +183,8 @@ def bump(actx, dcoll, t=0): / source_width**2)) -def main(ctx_factory, dim=2, order=3, visualize=False, lazy=False): +def main(ctx_factory, dim=2, order=3, + visualize=False, lazy=False, use_quad=False, use_nonaffine_mesh=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -186,11 +206,20 @@ def main(ctx_factory, dim=2, order=3, visualize=False, lazy=False): nel_1d = 16 if mesh_dist.is_mananger_rank(): - from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh( - a=(-0.5,)*dim, - b=(0.5,)*dim, - nelements_per_axis=(nel_1d,)*dim) + if use_nonaffine_mesh: + from meshmode.mesh.generation import generate_warped_rect_mesh + # FIXME: *generate_warped_rect_mesh* in meshmode warps a + # rectangle domain with hard-coded vertices at (-0.5,)*dim + # and (0.5,)*dim. Should extend the function interface to allow + # for specifying the end points directly. + mesh = generate_warped_rect_mesh(dim=dim, order=order, + nelements_side=nel_1d) + else: + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dim, + b=(0.5,)*dim, + nelements_per_axis=(nel_1d,)*dim) logger.info("%d elements", mesh.nelements) @@ -203,8 +232,23 @@ def main(ctx_factory, dim=2, order=3, visualize=False, lazy=False): else: local_mesh = mesh_dist.receive_mesh_part() - dcoll = DiscretizationCollection(actx, local_mesh, order=order, - mpi_communicator=comm) + from meshmode.discretization.poly_element import \ + QuadratureSimplexGroupFactory, \ + default_simplex_group_factory + dcoll = DiscretizationCollection( + actx, local_mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(base_dim=dim, order=order), + # High-order quadrature to integrate inner products of polynomials + # on warped geometry exactly (metric terms are also polynomial) + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(3*order), + } + ) + + if use_quad: + quad_tag = DISCR_TAG_QUAD + else: + quad_tag = None fields = WaveState( u=bump(actx, dcoll), @@ -217,7 +261,7 @@ def main(ctx_factory, dim=2, order=3, visualize=False, lazy=False): vis = make_visualizer(dcoll) def rhs(t, w): - return wave_operator(dcoll, c=c, w=w) + return wave_operator(dcoll, c=c, w=w, quad_tag=quad_tag) compiled_rhs = actx.compile(rhs) @@ -278,6 +322,8 @@ if __name__ == "__main__": parser.add_argument("--visualize", action="store_true") parser.add_argument("--lazy", action="store_true", help="switch to a lazy computation mode") + parser.add_argument("--quad", action="store_true") + parser.add_argument("--nonaffine", action="store_true") args = parser.parse_args() @@ -286,6 +332,8 @@ if __name__ == "__main__": dim=args.dim, order=args.order, visualize=args.visualize, - lazy=args.lazy) + lazy=args.lazy, + use_quad=args.quad, + use_nonaffine_mesh=args.nonaffine) # vim: foldmethod=marker diff --git a/grudge/array_context.py b/grudge/array_context.py index 2af7ece0..0f85a62b 100644 --- a/grudge/array_context.py +++ b/grudge/array_context.py @@ -32,6 +32,7 @@ from meshmode.array_context import ( ) from arraycontext.pytest import ( _PytestPyOpenCLArrayContextFactoryWithClass, + _PytestPytatoPyOpenCLArrayContextFactory, register_pytest_array_context_factory) @@ -56,6 +57,11 @@ class PytestPyOpenCLArrayContextFactory( actx_class = PyOpenCLArrayContext +class PytestPytatoPyOpenCLArrayContextFactory( + _PytestPytatoPyOpenCLArrayContextFactory): + actx_class = PytatoPyOpenCLArrayContext + + # deprecated class PytestPyOpenCLArrayContextFactoryWithHostScalars( _PytestPyOpenCLArrayContextFactoryWithClass): @@ -65,6 +71,8 @@ class PytestPyOpenCLArrayContextFactoryWithHostScalars( register_pytest_array_context_factory("grudge.pyopencl", PytestPyOpenCLArrayContextFactory) +register_pytest_array_context_factory("grudge.pytato-pyopencl", + PytestPytatoPyOpenCLArrayContextFactory) # }}} diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py index 34e84e51..492d8b7f 100644 --- a/grudge/geometry/metrics.py +++ b/grudge/geometry/metrics.py @@ -80,11 +80,38 @@ from arraycontext import register_multivector_as_array_container register_multivector_as_array_container() +def _geometry_to_quad_if_requested( + dcoll, inner_dd, dd, vec, _use_geoderiv_connection): + + def to_quad(vec): + if not dd.uses_quadrature(): + return vec + return dcoll.connection_from_dds(inner_dd, dd)(vec) + + # FIXME: At least for eager evaluation, this is somewhat inefficient, as + # all element groups vectors are upsampled to the quadrature grid, but then + # only the data for the non-affinely-mapped ones is used. + all_quad_vec = to_quad(vec) + + if not _use_geoderiv_connection: + return all_quad_vec + + return DOFArray( + vec.array_context, + tuple( + geoderiv_vec_i if megrp.is_affine else all_quad_vec_i + for megrp, geoderiv_vec_i, all_quad_vec_i in zip( + dcoll.discr_from_dd(inner_dd).mesh.groups, + dcoll._base_to_geoderiv_connection(inner_dd)(vec), + all_quad_vec))) + + # {{{ Metric computations def forward_metric_nth_derivative( actx: ArrayContext, dcoll: DiscretizationCollection, - xyz_axis, ref_axes, dd=None) -> DOFArray: + xyz_axis, ref_axes, dd=None, + *, _use_geoderiv_connection=False) -> DOFArray: r"""Pointwise metric derivatives representing repeated derivatives of the physical coordinate enumerated by *xyz_axis*: :math:`x_{\mathrm{xyz\_axis}}` with respect to the coordiantes on the reference element :math:`\xi_i`: @@ -108,6 +135,12 @@ def forward_metric_nth_derivative( by the axis index. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: If *True*, process returned + :class:`~meshmode.dof_array.DOFArray`\ s through + :meth:`~grudge.DiscretizationCollection._base_to_geoderiv_connection`. + This should be set based on whether the code using the result of this + function is able to make use of these arrays. (This is an internal + argument and not intended for use outside :mod:`grudge`.) :returns: a :class:`~meshmode.dof_array.DOFArray` containing the pointwise metric derivative at each nodal coordinate. """ @@ -139,15 +172,13 @@ def forward_metric_nth_derivative( thaw(dcoll.discr_from_dd(inner_dd).nodes(), actx)[xyz_axis] ) - if dd.uses_quadrature(): - vec = dcoll.connection_from_dds(inner_dd, dd)(vec) - - return vec + return _geometry_to_quad_if_requested( + dcoll, inner_dd, dd, vec, _use_geoderiv_connection) def forward_metric_derivative_vector( - actx: ArrayContext, dcoll: DiscretizationCollection, rst_axis, dd=None - ) -> np.ndarray: + actx: ArrayContext, dcoll: DiscretizationCollection, rst_axis, dd=None, + *, _use_geoderiv_connection=False) -> np.ndarray: r"""Computes an object array containing the forward metric derivatives of each physical coordinate. @@ -156,19 +187,23 @@ def forward_metric_derivative_vector( which will be taken. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s containing the pointwise metric derivatives at each nodal coordinate. """ return make_obj_array([ - forward_metric_nth_derivative(actx, dcoll, i, rst_axis, dd=dd) + forward_metric_nth_derivative( + actx, dcoll, i, rst_axis, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) for i in range(dcoll.ambient_dim) ] ) def forward_metric_derivative_mv( - actx: ArrayContext, dcoll: DiscretizationCollection, rst_axis, dd=None - ) -> MultiVector: + actx: ArrayContext, dcoll: DiscretizationCollection, rst_axis, dd=None, + *, _use_geoderiv_connection=False) -> MultiVector: r"""Computes a :class:`pymbolic.geometric_algebra.MultiVector` containing the forward metric derivatives of each physical coordinate. @@ -177,17 +212,21 @@ def forward_metric_derivative_mv( which will be taken. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a :class:`pymbolic.geometric_algebra.MultiVector` containing the forward metric derivatives in each physical coordinate. """ return MultiVector( - forward_metric_derivative_vector(actx, dcoll, rst_axis, dd=dd) + forward_metric_derivative_vector( + actx, dcoll, rst_axis, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) ) def forward_metric_derivative_mat( - actx: ArrayContext, dcoll: DiscretizationCollection, dd=None - ) -> np.ndarray: + actx: ArrayContext, dcoll: DiscretizationCollection, dd=None, + *, _use_geoderiv_connection=False) -> np.ndarray: r"""Computes the forward metric derivative matrix, also commonly called the Jacobian matrix, with entries defined as the forward metric derivatives: @@ -205,6 +244,8 @@ def forward_metric_derivative_mat( :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a matrix containing the evaluated forward metric derivatives of each physical coordinate, with respect to each reference coordinate. """ @@ -217,13 +258,15 @@ def forward_metric_derivative_mat( result = np.zeros((ambient_dim, dim), dtype=object) for j in range(dim): - result[:, j] = forward_metric_derivative_vector(actx, dcoll, j, dd=dd) + result[:, j] = forward_metric_derivative_vector( + actx, dcoll, j, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) return result def first_fundamental_form(actx: ArrayContext, dcoll: DiscretizationCollection, - dd=None) -> np.ndarray: + dd=None, *, _use_geoderiv_connection=False) -> np.ndarray: r"""Computes the first fundamental form using the Jacobian matrix: .. math:: @@ -243,25 +286,30 @@ def first_fundamental_form(actx: ArrayContext, dcoll: DiscretizationCollection, :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a matrix containing coefficients of the first fundamental form. """ if dd is None: dd = DD_VOLUME - mder = forward_metric_derivative_mat(actx, dcoll, dd=dd) + mder = forward_metric_derivative_mat( + actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection) return mder.T.dot(mder) def inverse_metric_derivative_mat( - actx: ArrayContext, dcoll: DiscretizationCollection, dd=None - ) -> np.ndarray: + actx: ArrayContext, dcoll: DiscretizationCollection, dd=None, + *, _use_geoderiv_connection=False) -> np.ndarray: r"""Computes the inverse metric derivative matrix, which is the inverse of the Jacobian (forward metric derivative) matrix. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a matrix containing the evaluated inverse metric derivatives. """ ambient_dim = dcoll.ambient_dim @@ -275,15 +323,16 @@ def inverse_metric_derivative_mat( for i in range(dim): for j in range(ambient_dim): result[i, j] = inverse_metric_derivative( - actx, dcoll, i, j, dd=dd + actx, dcoll, i, j, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection ) return result def inverse_first_fundamental_form( - actx: ArrayContext, dcoll: DiscretizationCollection, dd=None - ) -> np.ndarray: + actx: ArrayContext, dcoll: DiscretizationCollection, dd=None, + *, _use_geoderiv_connection=False) -> np.ndarray: r"""Computes the inverse of the first fundamental form: .. math:: @@ -300,6 +349,8 @@ def inverse_first_fundamental_form( :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a matrix containing coefficients of the inverse of the first fundamental form. """ @@ -309,10 +360,12 @@ def inverse_first_fundamental_form( dim = dcoll.discr_from_dd(dd).dim if dcoll.ambient_dim == dim: - inv_mder = inverse_metric_derivative_mat(actx, dcoll, dd=dd) + inv_mder = inverse_metric_derivative_mat( + actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection) inv_form1 = inv_mder.dot(inv_mder.T) else: - form1 = first_fundamental_form(actx, dcoll, dd=dd) + form1 = first_fundamental_form( + actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection) if dim == 1: inv_form1 = 1.0 / form1 @@ -329,7 +382,8 @@ def inverse_first_fundamental_form( def inverse_metric_derivative( - actx: ArrayContext, dcoll: DiscretizationCollection, rst_axis, xyz_axis, dd + actx: ArrayContext, dcoll: DiscretizationCollection, rst_axis, xyz_axis, dd, + *, _use_geoderiv_connection=False ) -> DOFArray: r"""Computes the inverse metric derivative of the physical coordinate enumerated by *xyz_axis* with respect to the @@ -339,6 +393,8 @@ def inverse_metric_derivative( :arg xyz_axis: an integer denoting the physical coordinate axis. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a :class:`~meshmode.dof_array.DOFArray` containing the inverse metric derivative at each nodal coordinate. """ @@ -350,8 +406,11 @@ def inverse_metric_derivative( "the derivative matrix is not square!" ) - par_vecs = [forward_metric_derivative_mv(actx, dcoll, rst, dd) - for rst in range(dim)] + par_vecs = [ + forward_metric_derivative_mv( + actx, dcoll, rst, dd, + _use_geoderiv_connection=_use_geoderiv_connection) + for rst in range(dim)] # Yay Cramer's rule! from functools import reduce, partial @@ -368,7 +427,9 @@ def inverse_metric_derivative( return outerprod(vecs) volume_pseudoscalar_inv = outerprod( - forward_metric_derivative_mv(actx, dcoll, rst_axis, dd) + forward_metric_derivative_mv( + actx, dcoll, rst_axis, dd, + _use_geoderiv_connection=_use_geoderiv_connection) for rst_axis in range(dim) ).inv() @@ -393,12 +454,8 @@ def inverse_surface_metric_derivative( :arg xyz_axis: an integer denoting the physical coordinate axis. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. - :arg _use_geoderiv_connection: If *True*, process returned - :class:`~meshmode.dof_array.DOFArray`\ s through - :meth:`~grudge.DiscretizationCollection._base_to_geoderiv_connection`. - This should be set based on whether the code using the result of this - function is able to make use of these arrays. (This is an internal - argument and not intended for use outside :mod:`grudge`.) + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a :class:`~meshmode.dof_array.DOFArray` containing the inverse metric derivative at each nodal coordinate. """ @@ -411,18 +468,17 @@ def inverse_surface_metric_derivative( if ambient_dim == dim: result = inverse_metric_derivative( - actx, dcoll, rst_axis, xyz_axis, dd=dd + actx, dcoll, rst_axis, xyz_axis, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection ) else: inv_form1 = inverse_first_fundamental_form(actx, dcoll, dd=dd) result = sum( inv_form1[rst_axis, d]*forward_metric_nth_derivative( - actx, dcoll, xyz_axis, d, dd=dd + actx, dcoll, xyz_axis, d, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection ) for d in range(dim)) - if _use_geoderiv_connection: - result = dcoll._base_to_geoderiv_connection(dd)(result) - return result @@ -440,12 +496,8 @@ def inverse_surface_metric_derivative_mat( :arg times_area_element: If *True*, each entry of the matrix is premultiplied with the value of :func:`area_element`, reflecting the typical use of the matrix in integrals evaluating weak derivatives. - :arg _use_geoderiv_connection: If *True*, process returned - :class:`~meshmode.dof_array.DOFArray`\ s through - :meth:`~grudge.DiscretizationCollection._base_to_geoderiv_connection`. - This should be set based on whether the code using the result of this - function is able to make use of these arrays. (This is an internal - argument and not intended for use outside :mod:`grudge`.) + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a :class:`~meshmode.dof_array.DOFArray` containing the inverse metric derivatives in per-group arrays of shape ``(xyz_dimension, rst_dimension, nelements, ndof)``. @@ -485,13 +537,13 @@ def _signed_face_ones( all_faces_conn = dcoll.connection_from_dds( DD_VOLUME, DOFDesc(dd.domain_tag) ) - signed_face_ones = dcoll.discr_from_dd(dd).zeros( + signed_ones = dcoll.discr_from_dd(dd.with_discr_tag(DISCR_TAG_BASE)).zeros( actx, dtype=dcoll.real_dtype ) + 1 from arraycontext import to_numpy, from_numpy, thaw - _signed_face_ones_numpy = to_numpy(signed_face_ones, actx) + _signed_face_ones_numpy = to_numpy(signed_ones, actx) for igrp, grp in enumerate(all_faces_conn.groups): for batch in grp.batches: @@ -504,13 +556,15 @@ def _signed_face_ones( def parametrization_derivative( - actx: ArrayContext, dcoll: DiscretizationCollection, dd - ) -> MultiVector: + actx: ArrayContext, dcoll: DiscretizationCollection, dd, + *, _use_geoderiv_connection=False) -> MultiVector: r"""Computes the product of forward metric derivatives spanning the tangent space with topological dimension *dim*. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a :class:`pymbolic.geometric_algebra.MultiVector` containing the product of metric derivatives. """ @@ -529,25 +583,31 @@ def parametrization_derivative( from pytools import product return product( - forward_metric_derivative_mv(actx, dcoll, rst_axis, dd) + forward_metric_derivative_mv( + actx, dcoll, rst_axis, dd, + _use_geoderiv_connection=_use_geoderiv_connection) for rst_axis in range(dim) ) def pseudoscalar(actx: ArrayContext, dcoll: DiscretizationCollection, - dd=None) -> MultiVector: + dd=None, *, _use_geoderiv_connection=False) -> MultiVector: r"""Computes the field of pseudoscalars for the domain/discretization identified by *dd*. :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: A :class:`~pymbolic.geometric_algebra.MultiVector` of :class:`~meshmode.dof_array.DOFArray`\ s. """ if dd is None: dd = DD_VOLUME - return parametrization_derivative(actx, dcoll, dd).project_max_grade() + return parametrization_derivative( + actx, dcoll, dd, + _use_geoderiv_connection=_use_geoderiv_connection).project_max_grade() def area_element( @@ -561,12 +621,8 @@ def area_element( :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization. - :arg _use_geoderiv_connection: If *True*, process returned - :class:`~meshmode.dof_array.DOFArray`\ s through - :meth:`~grudge.DiscretizationCollection._base_to_geoderiv_connection`. - This should be set based on whether the code using the result of this - function is able to make use of these arrays. (This is an internal - argument and not intended for use outside :mod:`grudge`.) + :arg _use_geoderiv_connection: For internal use. See + :func:`forward_metric_nth_derivative` for an explanation. :returns: a :class:`~meshmode.dof_array.DOFArray` containing the transformed volumes for each element. """ @@ -576,14 +632,9 @@ def area_element( @memoize_in(dcoll, (area_element, dd, _use_geoderiv_connection)) def _area_elements(): result = actx.np.sqrt( - pseudoscalar(actx, dcoll, dd=dd).norm_squared()) - - # NOTE: Don't interpolate 0-dim geometric factors to the - # "geometry discretization" - # (the metrics are just single scalars for facets in 1-D) - if dcoll.discr_from_dd(dd).dim != 0: - if _use_geoderiv_connection: - result = dcoll._base_to_geoderiv_connection(dd)(result) + pseudoscalar( + actx, dcoll, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection).norm_squared()) return freeze(result, actx) @@ -595,7 +646,8 @@ def area_element( # {{{ Surface normal vectors def rel_mv_normal( - actx: ArrayContext, dcoll: DiscretizationCollection, dd=None) -> MultiVector: + actx: ArrayContext, dcoll: DiscretizationCollection, dd=None, + *, _use_geoderiv_connection=False) -> MultiVector: r"""Computes surface normals at each nodal location as a :class:`~pymbolic.geometric_algebra.MultiVector` relative to the pseudoscalar of the discretization described by *dd*. @@ -608,7 +660,12 @@ def rel_mv_normal( # NOTE: Don't be tempted to add a sign here. As it is, it produces # exterior normals for positively oriented curves. - pder = pseudoscalar(actx, dcoll, dd=dd) / area_element(actx, dcoll, dd=dd) + pder = pseudoscalar( + actx, dcoll, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) \ + / area_element( + actx, dcoll, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) # Dorst Section 3.7.2 return pder << pder.I.inv() @@ -654,7 +711,9 @@ def mv_normal( f"their ambient dimension ({ambient_dim})") if dim == ambient_dim - 1: - result = rel_mv_normal(actx, dcoll, dd=dd) + result = rel_mv_normal( + actx, dcoll, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) else: # NOTE: In the case of (d - 2)-dimensional curves, we don't really have # enough information on the face to decide what an "exterior face normal" @@ -672,21 +731,18 @@ def mv_normal( project(dcoll, dof_desc.DD_VOLUME, dd, rel_mv_normal( actx, dcoll, - dd=dof_desc.DD_VOLUME + dd=dof_desc.DD_VOLUME, + _use_geoderiv_connection=_use_geoderiv_connection ).as_vector(dtype=object)) ) - pder = pseudoscalar(actx, dcoll, dd=dd) + pder = pseudoscalar( + actx, dcoll, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) mv = -(volm_normal ^ pder) << volm_normal.I.inv() result = mv / actx.np.sqrt(mv.norm_squared()) - # NOTE: Don't interpolate normals to the "geometry discretization" - # in 1-D (just a single value for every facet: +1 or -1) - if dim != 0: - if _use_geoderiv_connection: - result = dcoll._base_to_geoderiv_connection(dd)(result) - return freeze(result, actx) n = _normal() diff --git a/test/test_grudge.py b/test/test_grudge.py index 77566d13..06f0710c 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -34,6 +34,7 @@ pytest_generate_tests = pytest_generate_tests_for_array_contexts( from arraycontext.container.traversal import thaw from meshmode import _acf # noqa: F401 + from meshmode.dof_array import flat_norm import meshmode.mesh.generation as mgen @@ -52,49 +53,6 @@ import logging logger = logging.getLogger(__name__) -# {{{ inverse metric - -@pytest.mark.parametrize("dim", [2, 3]) -def test_inverse_metric(actx_factory, dim): - actx = actx_factory() - - mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, - nelements_per_axis=(6,)*dim, order=4) - - def m(x): - result = np.empty_like(x) - result[0] = ( - 1.5*x[0] + np.cos(x[0]) - + 0.1*np.sin(10*x[1])) - result[1] = ( - 0.05*np.cos(10*x[0]) - + 1.3*x[1] + np.sin(x[1])) - if len(x) == 3: - result[2] = x[2] - return result - - from meshmode.mesh.processing import map_mesh - mesh = map_mesh(mesh, m) - - dcoll = DiscretizationCollection(actx, mesh, order=4) - - from grudge.geometry import \ - forward_metric_derivative_mat, inverse_metric_derivative_mat - - mat = forward_metric_derivative_mat(actx, dcoll).dot( - inverse_metric_derivative_mat(actx, dcoll)) - - for i in range(mesh.dim): - for j in range(mesh.dim): - tgt = 1 if i == j else 0 - - err = flat_norm(mat[i, j] - tgt, ord=np.inf) - logger.info("error[%d, %d]: %.5e", i, j, err) - assert err < 1.0e-12, (i, j, err) - -# }}} - - # {{{ mass operator trig integration @pytest.mark.parametrize("ambient_dim", [1, 2, 3]) diff --git a/test/test_metrics.py b/test/test_metrics.py new file mode 100644 index 00000000..586b093a --- /dev/null +++ b/test/test_metrics.py @@ -0,0 +1,125 @@ +__copyright__ = """ +Copyright (C) 2015 Andreas Kloeckner +Copyright (C) 2021 University of Illinois Board of Trustees +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + +from grudge.array_context import ( + PytestPyOpenCLArrayContextFactory, + PytestPytatoPyOpenCLArrayContextFactory +) +from arraycontext import pytest_generate_tests_for_array_contexts +pytest_generate_tests = pytest_generate_tests_for_array_contexts( + [PytestPyOpenCLArrayContextFactory, + PytestPytatoPyOpenCLArrayContextFactory]) + +from meshmode.dof_array import flat_norm +import meshmode.mesh.generation as mgen + +from grudge import DiscretizationCollection + +import pytest + +import logging + +logger = logging.getLogger(__name__) + + +# {{{ inverse metric + +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("nonaffine", [False, True]) +@pytest.mark.parametrize("use_quad", [False, True]) +def test_inverse_metric(actx_factory, dim, nonaffine, use_quad): + actx = actx_factory() + + order = 3 + mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, + nelements_per_axis=(6,)*dim, order=order) + + if nonaffine: + def m(x): + result = np.empty_like(x) + result[0] = ( + 1.5*x[0] + np.cos(x[0]) + + 0.1*np.sin(10*x[1])) + result[1] = ( + 0.05*np.cos(10*x[0]) + + 1.3*x[1] + np.sin(x[1])) + if len(x) == 3: + result[2] = x[2] + return result + + from meshmode.mesh.processing import map_mesh + mesh = map_mesh(mesh, m) + + from grudge.dof_desc import as_dofdesc, DISCR_TAG_BASE, DISCR_TAG_QUAD + from meshmode.discretization.poly_element import \ + QuadratureSimplexGroupFactory, \ + default_simplex_group_factory + + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(base_dim=dim, order=order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order + 1), + } + ) + + from grudge.geometry import \ + forward_metric_derivative_mat, inverse_metric_derivative_mat + + dd = as_dofdesc("vol") + if use_quad: + dd = dd.with_discr_tag(DISCR_TAG_QUAD) + + mat = forward_metric_derivative_mat( + actx, dcoll, dd, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting).dot( + inverse_metric_derivative_mat( + actx, dcoll, dd, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting)) + + for i in range(mesh.dim): + for j in range(mesh.dim): + tgt = 1 if i == j else 0 + + err = actx.to_numpy(flat_norm(mat[i, j] - tgt, ord=np.inf)) + logger.info("error[%d, %d]: %.5e", i, j, err) + assert err < 1.0e-12, (i, j, err) + +# }}} + + +# You can test individual routines by typing +# $ python test_metrics.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) + +# vim: fdm=marker -- GitLab