diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 72a2932d0..d5fbc5670 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -32,11 +32,11 @@ import pyopencl as cl import pyopencl.tools as cl_tools from arraycontext import flatten -from meshmode.mesh import BTAG_ALL +from meshmode.mesh import BTAG_ALL, TensorProductElementGroup import grudge.dof_desc as dof_desc import grudge.op as op -from grudge.array_context import PyOpenCLArrayContext +from grudge.array_context import NumpyArrayContext, PytatoPyOpenCLArrayContext logger = logging.getLogger(__name__) @@ -96,24 +96,27 @@ def __call__(self, evt, basename, overwrite=True): # }}} -def main(ctx_factory, dim=2, order=4, visualize=False): +def main(ctx_factory, dim=1, order=2, lazy=False, + visualize=False, group_cls=TensorProductElementGroup): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext( - queue, - allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)), - force_device_scalars=True, - ) + + if lazy is False: + actx = NumpyArrayContext() + else: + actx = PytatoPyOpenCLArrayContext( + queue, + allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)),) # {{{ parameters # domain [-d/2, d/2]^dim d = 1.0 # number of points in each dimension - npoints = 20 + npoints = 10 # final time - final_time = 1.0 + final_time = 0.5 # velocity field c = np.array([0.5] * dim) @@ -129,7 +132,8 @@ def main(ctx_factory, dim=2, order=4, visualize=False): from meshmode.mesh.generation import generate_box_mesh mesh = generate_box_mesh( [np.linspace(-d/2, d/2, npoints) for _ in range(dim)], - order=order) + order=order, + group_cls=group_cls) from grudge.discretization import make_discretization_collection @@ -163,7 +167,10 @@ def u_analytic(x, t=0): def rhs(t, u): return adv_operator.operator(t, u) - dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u)) + rhs_compiled = actx.compile(rhs) + + # dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u)) + dt = 0.01 logger.info("Timestep size: %g", dt) @@ -172,7 +179,7 @@ def rhs(t, u): # {{{ time stepping from grudge.shortcuts import set_up_rk4 - dt_stepper = set_up_rk4("u", float(dt), u, rhs) + dt_stepper = set_up_rk4("u", float(dt), u, rhs_compiled) plot = Plotter(actx, dcoll, order, visualize=visualize, ylim=[-1.1, 1.1]) @@ -200,13 +207,16 @@ def rhs(t, u): import argparse parser = argparse.ArgumentParser() - parser.add_argument("--dim", default=2, type=int) - parser.add_argument("--order", default=4, type=int) + parser.add_argument("--dim", default=1, type=int) + parser.add_argument("--order", default=2, type=int) parser.add_argument("--visualize", action="store_true") + parser.add_argument("--lazy", action="store_true") + parser.add_argument("--tp-elements", action="store_true") args = parser.parse_args() logging.basicConfig(level=logging.INFO) main(cl.create_some_context, dim=args.dim, order=args.order, - visualize=args.visualize) + visualize=args.visualize, + lazy=args.lazy) diff --git a/examples/euler/acoustic_pulse.py b/examples/euler/acoustic_pulse.py index c44a6322e..0388c695c 100644 --- a/examples/euler/acoustic_pulse.py +++ b/examples/euler/acoustic_pulse.py @@ -29,12 +29,16 @@ import pyopencl as cl import pyopencl.tools as cl_tools -from arraycontext import ArrayContext -from meshmode.mesh import BTAG_ALL +from arraycontext import ArrayContext, NumpyArrayContext +from meshmode.discretization.poly_element import ( + InterpolatoryEdgeClusteredGroupFactory, + QuadratureGroupFactory, +) +from meshmode.mesh import BTAG_ALL, SimplexElementGroup, TensorProductElementGroup from pytools.obj_array import make_obj_array import grudge.op as op -from grudge.array_context import PyOpenCLArrayContext, PytatoPyOpenCLArrayContext +from grudge.array_context import PytatoPyOpenCLArrayContext from grudge.models.euler import ConservedEulerField, EulerOperator, InviscidWallBC from grudge.shortcuts import rk4_step @@ -106,7 +110,8 @@ def run_acoustic_pulse(actx, final_time=1, resolution=16, overintegration=False, - visualize=False): + visualize=False, + tensor_product_elements=False): # eos-related parameters gamma = 1.4 @@ -115,18 +120,19 @@ def run_acoustic_pulse(actx, from meshmode.mesh.generation import generate_regular_rect_mesh + if tensor_product_elements: + group_cls = TensorProductElementGroup + else: + group_cls = SimplexElementGroup + dim = 2 box_ll = -0.5 box_ur = 0.5 mesh = generate_regular_rect_mesh( a=(box_ll,)*dim, b=(box_ur,)*dim, - nelements_per_axis=(resolution,)*dim) - - from meshmode.discretization.poly_element import ( - QuadratureSimplexGroupFactory, - default_simplex_group_factory, - ) + nelements_per_axis=(resolution,)*dim, + group_cls=group_cls) from grudge.discretization import make_discretization_collection from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD @@ -141,9 +147,8 @@ def run_acoustic_pulse(actx, dcoll = make_discretization_collection( actx, mesh, discr_tag_to_group_factory={ - DISCR_TAG_BASE: default_simplex_group_factory( - base_dim=mesh.dim, order=order), - DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) + DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order=order), + DISCR_TAG_QUAD: QuadratureGroupFactory(2*order) } ) @@ -182,12 +187,20 @@ def rhs(t, q): # {{{ time stepping + import time + step = 0 t = 0.0 + elapsed = 0.0 while t < final_time: if step % 10 == 0: norm_q = actx.to_numpy(op.norm(dcoll, fields, 2)) - logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q) + if step != 0: + logger.info("[%04d] t = %.5f |q| = %.5e time per step = %.5f", + step, t, norm_q, elapsed / step) + else: + logger.info("[%04d] t = %.5f |q| = %.5e time per step = %.5f", + step, t, norm_q, 0) if visualize: vis.write_vtk_file( f"{exp_name}-{step:04d}.vtu", @@ -199,8 +212,10 @@ def rhs(t, q): ) assert norm_q < 5 + start = time.time() fields = actx.thaw(actx.freeze(fields)) fields = rk4_step(fields, t, dt, compiled_rhs) + elapsed += time.time() - start t += dt step += 1 @@ -208,7 +223,8 @@ def rhs(t, q): def main(ctx_factory, order=3, final_time=1, resolution=16, - overintegration=False, visualize=False, lazy=False): + overintegration=False, visualize=False, lazy=False, + tensor_product_elements=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -218,11 +234,7 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)), ) else: - actx = PyOpenCLArrayContext( - queue, - allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)), - force_device_scalars=True, - ) + actx = NumpyArrayContext() run_acoustic_pulse( actx, @@ -230,7 +242,8 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, resolution=resolution, overintegration=overintegration, final_time=final_time, - visualize=visualize + visualize=visualize, + tensor_product_elements=tensor_product_elements ) @@ -238,9 +251,10 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, import argparse parser = argparse.ArgumentParser() - parser.add_argument("--order", default=3, type=int) + parser.add_argument("--tpe", action="store_true") + parser.add_argument("--order", default=2, type=int) parser.add_argument("--tfinal", default=0.1, type=float) - parser.add_argument("--resolution", default=16, type=int) + parser.add_argument("--resolution", default=4, type=int) parser.add_argument("--oi", action="store_true", help="use overintegration") parser.add_argument("--visualize", action="store_true", @@ -256,4 +270,5 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, resolution=args.resolution, overintegration=args.oi, visualize=args.visualize, - lazy=args.lazy) + lazy=args.lazy, + tensor_product_elements=args.tpe) diff --git a/grudge/array_context.py b/grudge/array_context.py index aca1edc08..c6d1d05f0 100644 --- a/grudge/array_context.py +++ b/grudge/array_context.py @@ -38,13 +38,30 @@ from typing import TYPE_CHECKING, Any, Optional from warnings import warn +import pytato as pt +from pytato.analysis import get_num_nodes + +import loopy as lp from meshmode.array_context import ( PyOpenCLArrayContext as _PyOpenCLArrayContextBase, PytatoPyOpenCLArrayContext as _PytatoPyOpenCLArrayContextBase, ) +from meshmode.transform_metadata import ( + DiscretizationDOFAxisTag, + DiscretizationElementAxisTag, + DiscretizationEntityAxisTag, +) from pytools import to_identifier from pytools.tag import Tag +from grudge.transform.mappers import ( + tensor_product_algebraic_transforms, +) +from grudge.transform.metadata import ( + OutputIsTensorProductDOFArrayOrdered, + TensorProductDOFAxisTag, +) + logger = logging.getLogger(__name__) @@ -132,16 +149,101 @@ def __init__(self, queue: "pyopencl.CommandQueue", super().__init__(queue, allocator, wait_event_queue_length, force_device_scalars) + def transform_loopy_program(self, t_unit): + knl = t_unit.default_entrypoint + + # {{{ process tensor product specific metadata + + # NOTE: This differs from the lazy case b/c we don't have access to axis + # tags that can be manipulated pre-execution. In eager, we update + # strides/loop nest ordering for the output array + if knl.tags_of_type(OutputIsTensorProductDOFArrayOrdered): + new_args = [] + for arg in knl.args: + if arg.is_output: + arg = arg.copy(dim_tags=( + f"N{len(arg.shape)-1}," + + ",".join(f"N{i}" + for i in range(len(arg.shape)-1)) + )) + + new_args.append(arg) + + knl = knl.copy(args=new_args) + t_unit = t_unit.with_kernel(knl) + + # }}} + + return super().transform_loopy_program(t_unit) + # }}} # {{{ pytato +def deduplicate_data_wrappers(dag): + import pytato as pt + data_wrapper_cache = {} + data_wrappers_encountered = 0 + + def cached_data_wrapper_if_present(ary): + nonlocal data_wrappers_encountered + + if isinstance(ary, pt.DataWrapper): + + data_wrappers_encountered += 1 + cache_key = (ary.data.base_data.int_ptr, ary.data.offset, + ary.shape, ary.data.strides) + try: + result = data_wrapper_cache[cache_key] + except KeyError: + result = ary + data_wrapper_cache[cache_key] = result + + return result + else: + return ary + + dag = pt.transform.map_and_copy(dag, cached_data_wrapper_if_present) + + if data_wrappers_encountered: + logger.info("data wrapper de-duplication: " + "%d encountered, %d kept, %d eliminated", + data_wrappers_encountered, + len(data_wrapper_cache), + data_wrappers_encountered - len(data_wrapper_cache)) + + return dag + + +def remove_redundant_tensor_product_reshapes(ary): + # FIXME: variable names can be more clear + if isinstance(ary, pt.Reshape): + if isinstance(ary.array, pt.Reshape): + if ary.array.array.shape == ary.shape: + return ary.array.array + + return ary + + +def remove_redundant_index_lambda_expressions(ary): + # FIXME: this can be made much more robust + if isinstance(ary, pt.IndexLambda): + if len(ary.shape) >= 3: + if 0.0 in ary.expr.children: + return list(ary.bindings.values())[0] + + return ary + + class PytatoPyOpenCLArrayContext(_PytatoPyOpenCLArrayContextBase): """Inherits from :class:`meshmode.array_context.PytatoPyOpenCLArrayContext`. Extends it to understand :mod:`grudge`-specific transform metadata. (Of which there isn't any, for now.) """ + + dot_codes = [] + def __init__(self, queue, allocator=None, *, compile_trace_callback: Callable[[Any, str, Any], None] | None @@ -162,6 +264,70 @@ def __init__(self, queue, allocator=None, super().__init__(queue, allocator, compile_trace_callback=compile_trace_callback) + def transform_dag(self, dag): + + if 1: + dag = tensor_product_algebraic_transforms(dag) + + dag = pt.transform.map_and_copy( + dag, remove_redundant_tensor_product_reshapes) + + dag = pt.transform.map_and_copy( + dag, remove_redundant_index_lambda_expressions) + + dag = deduplicate_data_wrappers(dag) + + from grudge.transform.metadata import unify_discretization_entity_tags + logger.info("transform_dag.infer_axes_tags") + dag = unify_discretization_entity_tags(dag) + + def eliminate_reshapes_of_data_wrappers(ary): + if (isinstance(ary, pt.Reshape) + and isinstance(ary.array, pt.DataWrapper)): + return pt.make_data_wrapper(ary.array.data.reshape(ary.shape), + tags=ary.tags, + axes=ary.axes) + else: + return ary + + dag = pt.transform.map_and_copy(dag, + eliminate_reshapes_of_data_wrappers) + + logger.info("transform_dag.infer_axes_tags") + from grudge.transform.metadata import unify_discretization_entity_tags + dag = unify_discretization_entity_tags(dag) + + def _untag_impl_stored(expr): + if isinstance(expr, pt.InputArgumentBase): + return expr + else: + return expr.without_tags(pt.tags.ImplStored(), + verify_existence=False) + + dag = pt.make_dict_of_named_arrays({ + name: _untag_impl_stored(named_ary.expr) + for name, named_ary in dag.items()}) + + self.dot_codes.append(pt.get_dot_graph(dag)) + + return dag + + def transform_loopy_program(self, t_unit): + knl = t_unit.default_entrypoint + + redn_inames = { + iname + for insn in knl.instructions + for iname in insn.reduction_inames() + } + + for iname in knl.inames: + if knl.inames[iname].tags_of_type(DiscretizationElementAxisTag): + if iname not in redn_inames: + knl = lp.tag_inames(knl, {iname: "g.0"}) + + return t_unit.with_kernel(knl) + # }}} @@ -226,7 +392,7 @@ def _dag_to_compiled_func(self, dict_of_named_arrays, "transform_dag.infer_axes_tags[pre-partition]"): from meshmode.transform_metadata import DiscretizationEntityAxisTag dict_of_named_arrays = pt.unify_axes_tags( - dict_of_named_arrays, + dict_of_named_arrays, tag_t=DiscretizationEntityAxisTag, ) @@ -550,6 +716,7 @@ def __call__(self): class PytestNumpyArrayContextFactory(_PytestNumpyArrayContextFactory): + from arraycontext import NumpyArrayContext actx_class = NumpyArrayContext def __call__(self): diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py index e7b726e2b..4c4833726 100644 --- a/grudge/geometry/metrics.py +++ b/grudge/geometry/metrics.py @@ -520,21 +520,33 @@ def inverse_surface_metric_derivative_mat( @memoize_in(dcoll, (inverse_surface_metric_derivative_mat, dd, times_area_element, _use_geoderiv_connection)) def _inv_surf_metric_deriv(): + # FIXME: two branches to avoid multiplication by 1 (for now) if times_area_element: + from warnings import warn + warn("Passing `times_area_element` to " + "`inverse_surface_metric_derivative_mat` is deprecated and " + "will stop working in 2025", stacklevel=1) + multiplier = area_element(actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection) + mat = actx.np.stack([ + actx.np.stack([ + multiplier * inverse_surface_metric_derivative( + actx, dcoll, + rst_axis, xyz_axis, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) + for rst_axis in range(dcoll.dim)]) + for xyz_axis in range(dcoll.ambient_dim)]) + else: - multiplier = 1 - - mat = actx.np.stack([ - actx.np.stack([ - multiplier - * inverse_surface_metric_derivative( - actx, dcoll, - rst_axis, xyz_axis, dd=dd, - _use_geoderiv_connection=_use_geoderiv_connection) - for rst_axis in range(dcoll.dim)]) - for xyz_axis in range(dcoll.ambient_dim)]) + mat = actx.np.stack([ + actx.np.stack([ + inverse_surface_metric_derivative( + actx, dcoll, + rst_axis, xyz_axis, dd=dd, + _use_geoderiv_connection=_use_geoderiv_connection) + for rst_axis in range(dcoll.dim)]) + for xyz_axis in range(dcoll.ambient_dim)]) return actx.freeze( actx.tag(NameHint(f"inv_metric_deriv_{dd.as_identifier()}"), diff --git a/grudge/matrices.py b/grudge/matrices.py new file mode 100644 index 000000000..d25a2af43 --- /dev/null +++ b/grudge/matrices.py @@ -0,0 +1,466 @@ +""" +Core DG Matrices +^^^^^^^^^^^^^^^^ + +Strong differentiation +---------------------- + +.. autofunction:: reference_derivative_matrices + +Weak differentiation +-------------------- + +.. autofunction:: reference_stiffness_transpose_matrices + +Mass, inverse mass, and face mass matrices +------------------------------------------ + +.. autofunction:: reference_mass_matrix +.. autofunction:: reference_inverse_mass_matrix +.. autofunction:: reference_face_mass_matrix +""" + +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2024 Addison Alvey-Blanco +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. +""" + +from collections.abc import Mapping + +import numpy as np +import numpy.linalg as la + +import modepy as mp +from arraycontext import ArrayContext, ArrayOrContainer, tag_axes +from arraycontext.metadata import NameHint +from meshmode.discretization import ( + DiscretizationDOFAxisTag, + InterpolatoryElementGroupBase, + NodalElementGroupBase, +) +from meshmode.transform_metadata import ( + DiscretizationAmbientDimAxisTag, + DiscretizationFaceAxisTag, +) +from pytools import keyed_memoize_on_first_arg +from pytools.tag import Tag + +from grudge.tools import ( + get_accurate_quadrature_rule, + get_basis_for_face_group, + get_element_group_basis, + get_element_group_nodes, + get_faces_for_volume_group, + get_quadrature_for_face, +) +from grudge.transform.metadata import ( + TensorProductMassOperatorInverseTag, + TensorProductMassOperatorTag, + TensorProductStiffnessOperatorTag, + get_dof_axis_tag_type +) + + +@keyed_memoize_on_first_arg( + lambda output_group, input_group, use_tensor_product_fast_eval: ( + input_group.discretization_key(), + output_group.discretization_key(), + use_tensor_product_fast_eval + ) +) +def reference_derivative_matrices( + actx: ArrayContext, + input_group: InterpolatoryElementGroupBase, + output_group: NodalElementGroupBase, + use_tensor_product_fast_eval: bool = True, + axis_tags: Mapping[int, tuple[Tag]] | None = None, + ary_tags: tuple[Tag] | None = None, + ) -> ArrayOrContainer: + """ + Computes all reference derivative matrices. See + :func:`~grudge.matrices.reference_derivative_matrix` for more information. + """ + + if axis_tags is None: + dof_axis_tag = get_dof_axis_tag_type(output_group, + use_tensor_product_fast_eval) + + if not use_tensor_product_fast_eval: + axis_tags = {0: (DiscretizationAmbientDimAxisTag(),)} + axis_tags.update({ + i+1: (dof_axis_tag(),) for i in range(2) + }) # type: ignore + else: + axis_tags = ({ + i: (dof_axis_tag(),) for i in range(2) + }) # type: ignore + if ary_tags is None: + ary_tags = (NameHint("ref_deriv_mat"),) + + basis = get_element_group_basis( + input_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + input_nodes = get_element_group_nodes( + input_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + output_nodes = get_element_group_nodes( + output_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy( + np.asarray( + mp.diff_matrices(basis, + output_nodes, + from_nodes=input_nodes), + order="C") + ) + ) + ) + + +@keyed_memoize_on_first_arg( + lambda output_group, input_group, use_tensor_product_fast_eval: ( + input_group.discretization_key(), + output_group.discretization_key(), + use_tensor_product_fast_eval + ) +) +def reference_stiffness_transpose_matrices( + actx: ArrayContext, + input_group: NodalElementGroupBase, + output_group: InterpolatoryElementGroupBase, + use_tensor_product_fast_eval: bool = True, + axis_tags: Mapping[int, tuple[Tag]] | None = None, + ary_tags: tuple[Tag] | None = None, + ) -> ArrayOrContainer: + + if axis_tags is None: + dof_axis_tag = get_dof_axis_tag_type(output_group, + use_tensor_product_fast_eval) + if not use_tensor_product_fast_eval: + axis_tags = {0: (DiscretizationAmbientDimAxisTag(),)} + axis_tags.update({ + i+1: (dof_axis_tag(),) for i in range(2) + }) # type: ignore + else: + axis_tags = { + i: (dof_axis_tag(),) for i in range(2) + } + if ary_tags is None: + ary_tags = (NameHint("stiff_t"),) + if use_tensor_product_fast_eval: + ary_tags += (TensorProductStiffnessOperatorTag(),) # type: ignore + + basis = get_element_group_basis( + output_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + nodes = get_element_group_nodes( + output_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + quadrature = get_accurate_quadrature_rule( + input_group, + required_exactness=2*output_group.order, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + if use_tensor_product_fast_eval: + num_matrices = 1 + else: + num_matrices = input_group.dim + + if input_group == output_group: + stiffness_t = np.asarray([ + mp.nodal_quad_bilinear_form( + quadrature=quadrature, + test_basis=basis, + trial_basis=basis, + input_nodes=nodes, + test_derivative_ax=rst_axis + ) + for rst_axis in range(num_matrices) + ], order="C" + ) + + if use_tensor_product_fast_eval: + stiffness_t = stiffness_t[0] + + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy(stiffness_t) + ) + ) + + stiffness_t = np.asarray([ + mp.nodal_quad_operator( + quadrature=quadrature, + test_basis=basis, + nodes=nodes, + test_derivative_ax=rst_axis + ) + for rst_axis in range(num_matrices) + ], order="C") + + if use_tensor_product_fast_eval: + stiffness_t = stiffness_t[0] + + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy(stiffness_t) + ) + ) + + +@keyed_memoize_on_first_arg( + lambda output_group, input_group, use_tensor_product_fast_eval: ( + input_group.discretization_key(), + output_group.discretization_key(), + use_tensor_product_fast_eval + ) +) +def reference_mass_matrix( + actx: ArrayContext, + input_group: NodalElementGroupBase, + output_group: InterpolatoryElementGroupBase, + use_tensor_product_fast_eval: bool = False, + axis_tags: Mapping[int, tuple[Tag]] | None = None, + ary_tags: tuple[Tag] | None = None + ) -> ArrayOrContainer: + + if axis_tags is None: + dof_axis_tag = get_dof_axis_tag_type(output_group, + use_tensor_product_fast_eval) + axis_tags = {i: (dof_axis_tag(),) for i in range(2)} + if ary_tags is None: + ary_tags = (NameHint("ref_mass"),) + if use_tensor_product_fast_eval: + ary_tags += (TensorProductMassOperatorTag(),) # type: ignore + + basis = get_element_group_basis( + output_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + nodes = get_element_group_nodes( + output_group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + quadrature = get_accurate_quadrature_rule( + input_group, + required_exactness=2*output_group.order, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + if input_group == output_group: + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy( + mp.nodal_quad_bilinear_form( + quadrature=quadrature, + test_basis=basis, + trial_basis=basis, + input_nodes=nodes + ) + ) + ) + ) + + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy( + mp.nodal_quad_operator( + quadrature=quadrature, + test_basis=basis, + nodes=nodes + ) + ) + ) + ) + + +@keyed_memoize_on_first_arg( + lambda group, use_tensor_product_fast_eval: ( + group.discretization_key(), + use_tensor_product_fast_eval + ) +) +def reference_inverse_mass_matrix( + actx: ArrayContext, + group: InterpolatoryElementGroupBase, + use_tensor_product_fast_eval: bool = False, + axis_tags: Mapping[int, tuple[Tag]] | None = None, + ary_tags: tuple[Tag] | None = None + ) -> ArrayOrContainer: + + if axis_tags is None: + dof_axis_tag = get_dof_axis_tag_type(group, + use_tensor_product_fast_eval) + axis_tags = {i: (dof_axis_tag(),) for i in range(2)} + if ary_tags is None: + ary_tags = (NameHint("ref_inv_mass"),) + if use_tensor_product_fast_eval: + ary_tags += (TensorProductMassOperatorInverseTag(),) # type: ignore + + basis = get_element_group_basis( + group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + nodes = get_element_group_nodes( + group, use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + quadrature = get_accurate_quadrature_rule( + group, + required_exactness=2*(group.order + 1), + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy( + la.inv( + mp.nodal_quad_bilinear_form( + quadrature=quadrature, + test_basis=basis, + trial_basis=basis, + input_nodes=nodes + ) + ) + ) + ) + ) + + +@keyed_memoize_on_first_arg( + lambda face_group, vol_group, dtype, use_tensor_product_fast_eval: ( + face_group.discretization_key(), + vol_group.discretization_key(), + dtype, + use_tensor_product_fast_eval + ) +) +def reference_face_mass_matrix( + actx: ArrayContext, + face_group: NodalElementGroupBase, + vol_group: InterpolatoryElementGroupBase, + dtype, + ary_tags: tuple[Tag] | None = None, + axis_tags: Mapping[int, Tag] | None = None, + use_tensor_product_fast_eval: bool = True, + ) -> ArrayOrContainer: + + use_tensor_product_fast_eval = False + + if axis_tags is None: + axis_tags = { + 0: (DiscretizationDOFAxisTag(),), + 1: (DiscretizationFaceAxisTag(),), + 2: (DiscretizationDOFAxisTag(),) + } # type: ignore + if ary_tags is None: + ary_tags = (NameHint("face_mass"),) + + face_mass = np.empty( + (vol_group.nunit_dofs, + vol_group.mesh_el_group.nfaces, + face_group.nunit_dofs), + dtype=dtype + ) + + faces = get_faces_for_volume_group( + vol_group, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + vol_nodes = get_element_group_nodes( + vol_group, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + face_nodes = get_element_group_nodes( + face_group, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + face_basis = get_basis_for_face_group( + face_group, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + vol_basis = get_element_group_basis( + vol_group, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + for iface, face in enumerate(faces): + face_quadrature = get_quadrature_for_face( + face_group, + face, + required_exactness=(vol_group.order + face_group.order), + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + if face_basis is not None: + face_mass[:, iface, :] = mp.nodal_quad_bilinear_form( + quadrature=face_quadrature, + trial_basis=face_basis, + test_basis=vol_basis, + input_nodes=face_nodes, + output_nodes=vol_nodes, + mapping_function=face.map_to_volume + ) + + else: + face_mass[:, iface, :] = mp.nodal_quad_operator( + quadrature=face_quadrature, + test_basis=vol_basis, + nodes=vol_nodes, + mapping_function=face.map_to_volume + ) + + return actx.tag( + ary_tags, + tag_axes( + actx, + axis_tags, + actx.from_numpy(face_mass) + ) + ) + +# vim: foldmethod=marker diff --git a/grudge/models/euler.py b/grudge/models/euler.py index ef55ad63e..9382abbd5 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -358,7 +358,7 @@ def interp_to_quad(u): interface_fluxes = interface_fluxes + bc_fluxes return op.inverse_mass( - dcoll, + dcoll, dq, volume_fluxes - op.face_mass(dcoll, df, interface_fluxes) ) diff --git a/grudge/op.py b/grudge/op.py index 1937f35d4..125738cde 100644 --- a/grudge/op.py +++ b/grudge/op.py @@ -43,8 +43,6 @@ from __future__ import annotations -from meshmode.discretization import InterpolatoryElementGroupBase, NodalElementGroupBase - __copyright__ = """ Copyright (C) 2021 Andreas Kloeckner @@ -72,20 +70,36 @@ """ +from collections.abc import Callable, Iterable from functools import partial import numpy as np -import modepy as mp -from arraycontext import ArrayContext, ArrayOrContainer, map_array_container, tag_axes +from arraycontext import ( + ArrayContext, + ArrayOrContainer, + map_array_container, + tag_axes, +) +from meshmode.discretization import ( + Discretization, + ElementGroupBase, + InterpolatoryElementGroupBase, + NodalElementGroupBase, +) +from meshmode.discretization.poly_element import ( + TensorProductElementGroupBase, +) from meshmode.dof_array import DOFArray from meshmode.transform_metadata import ( + DiscretizationAmbientDimAxisTag, DiscretizationDOFAxisTag, DiscretizationElementAxisTag, - DiscretizationFaceAxisTag, - FirstAxisIsElementsTag, ) -from pytools import keyed_memoize_in +from modepy.tools import ( + reshape_array_for_tensor_product_space as fold, + unreshape_array_for_tensor_product_space as unfold, +) from pytools.obj_array import make_obj_array import grudge.dof_desc as dof_desc @@ -93,12 +107,21 @@ from grudge.dof_desc import ( DD_VOLUME_ALL, DISCR_TAG_BASE, + DISCR_TAG_QUAD, FACE_RESTR_ALL, DOFDesc, VolumeDomainTag, as_dofdesc, ) +from grudge.geometry import area_element, inverse_surface_metric_derivative_mat from grudge.interpolation import interp +from grudge.matrices import ( + reference_derivative_matrices, + reference_face_mass_matrix, + reference_inverse_mass_matrix, + reference_mass_matrix, + reference_stiffness_transpose_matrices, +) from grudge.projection import project from grudge.reductions import ( elementwise_integral, @@ -114,6 +137,7 @@ nodal_sum_loc, norm, ) +from grudge.tools import rec_map_subarrays from grudge.trace_pair import ( bdry_trace_pair, bv_trace_pair, @@ -125,6 +149,9 @@ project_tracepair, tracepair_with_discr_tag, ) +from grudge.transform.metadata import ( + TensorProductDOFAxisTag, +) __all__ = ( @@ -163,161 +190,280 @@ ) -# {{{ common derivative "kernels" +# {{{ tensor product operator application helper -def _single_axis_derivative_kernel( - actx, out_discr, in_discr, get_diff_mat, inv_jac_mat, xyz_axis, vec, - *, metric_in_matvec): - # This gets used from both the strong and the weak derivative. These differ - # in three ways: - # - which differentiation matrix gets used, - # - whether inv_jac_mat is pre-multiplied by a factor that includes the - # area element, and - # - whether the chain rule terms ("inv_jac_mat") sit outside (strong) - # or inside (weak) the matrix-vector product that carries out the - # derivative, cf. "metric_in_matvec". - return DOFArray( +def _single_axis_contraction( + actx: ArrayContext, + dim: int, + axis: int, + operator: ArrayOrContainer, + data: DOFArray, + tagged=None, + arg_names=None, + ) -> ArrayOrContainer: + """ + Generic routine to apply a 1D operator to a particular axis of *data*. + + The einsum specification is constructed based on the dimension of the + problem and can support up to 1 reduction axis and 22 non-reduction + (DOF) axes. The element axis is not counted since it is reserved. + """ + data = tag_axes( actx, - data=tuple( - # r for rst axis - actx.einsum("rej,rij,ej->ei" if metric_in_matvec else "rei,rij,ej->ei", - ijm_i[xyz_axis], - get_diff_mat( - actx, - out_element_group=out_grp, - in_element_group=in_grp), - vec_i, - arg_names=("inv_jac_t", "ref_stiffT_mat", "vec", ), - tagged=(FirstAxisIsElementsTag(),)) - - for out_grp, in_grp, vec_i, ijm_i in zip( - out_discr.groups, - in_discr.groups, - vec, - inv_jac_mat, strict=True))) - - -def _gradient_kernel(actx, out_discr, in_discr, get_diff_mat, inv_jac_mat, vec, - *, metric_in_matvec): - # See _single_axis_derivative_kernel for comments on the usage scenarios - # (both strong and weak derivative) and their differences. - per_group_grads = [ - # r for rst axis - # x for xyz axis - actx.einsum("xrej,rij,ej->xei" if metric_in_matvec else "xrei,rij,ej->xei", - ijm_i, - get_diff_mat( - actx, - out_element_group=out_grp, - in_element_group=in_grp - ), - vec_i, - arg_names=("inv_jac_t", "ref_stiffT_mat", "vec"), - tagged=(FirstAxisIsElementsTag(),)) - for out_grp, in_grp, vec_i, ijm_i in zip( - out_discr.groups, in_discr.groups, vec, - inv_jac_mat, strict=True)] + { + i: (DiscretizationElementAxisTag() if i == 0 else + TensorProductDOFAxisTag(i-1)) + for i in range(dim+1) + }, + data) + + operator_spec = "ij" + data_spec = f"e{"abcdfghklmn"[:axis]}j{"opqrstuvwxy"[:dim-axis-1]}" + out_spec = f"e{"abcdfghklmn"[:axis]}i{"opqrstuvwxy"[:dim-axis-1]}" + spec = operator_spec + "," + data_spec + "->" + out_spec + + return tag_axes( + actx, + { + i: (DiscretizationElementAxisTag() if i == 0 else + TensorProductDOFAxisTag(i-1)) + for i in range(dim+1) + }, + actx.einsum(spec, operator, data, arg_names=arg_names, + tagged=tagged)) + + +# }}} + + +# {{{ Derivative operator kernels + +def _single_axis_derivative_kernel( + actx: ArrayContext, + output_discr: Discretization, + input_discr: Discretization, + inv_jac_mat: Iterable, + xyz_axis: int, + vec: DOFArray, + compute_single_axis_derivative: Callable, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: + """ + Computes the *xyz_axis*th derivative of *vec* using + *compute_single_axis_derivative*. + + Strong and weak dispatch routines pass the corresponding + *compute_single_axis_derivative* routine so that the correct form is + computed. + """ + + group_data = [] + for output_group, input_group, vec_i, ijm_i in zip( + output_discr.groups, input_discr.groups, vec, inv_jac_mat, + strict=False): + + group_data.append(compute_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec_i, ijm_i, + use_tensor_product_fast_eval=use_tensor_product_fast_eval)) + + return DOFArray(actx, data=tuple(group_data)) + + +def _gradient_kernel( + actx: ArrayContext, + output_discr: Discretization, + input_discr: Discretization, + inv_jac_mat: Iterable, + vec: DOFArray, + compute_single_axis_derivative: Callable, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: + """ + Computes gradient of *vec* using *compute_single_axis_derivative* for each + entry in the gradient. + + Strong and weak dispatch routines provide the proper + *compute_single_axis_derivative* routine. + """ + + per_group_grads = [] + for output_group, input_group, vec_i, ijm_i in zip( + output_discr.groups, input_discr.groups, vec, inv_jac_mat, + strict=False): + + per_group_grads.append([compute_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec_i, ijm_i, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) + for xyz_axis in range(input_group.dim) + ]) return make_obj_array([ - DOFArray(actx, data=tuple([ # noqa: C409 - pgg_i[xyz_axis] for pgg_i in per_group_grads - ])) - for xyz_axis in range(out_discr.ambient_dim)]) - - -def _divergence_kernel(actx, out_discr, in_discr, get_diff_mat, inv_jac_mat, vec, - *, metric_in_matvec): - # See _single_axis_derivative_kernel for comments on the usage scenarios - # (both strong and weak derivative) and their differences. - per_group_divs = [ - # r for rst axis - # x for xyz axis - actx.einsum("xrej,rij,xej->ei" if metric_in_matvec else "xrei,rij,xej->ei", - ijm_i, - get_diff_mat( - actx, - out_element_group=out_grp, - in_element_group=in_grp - ), - vec_i, - arg_names=("inv_jac_t", "ref_stiffT_mat", "vec"), - tagged=(FirstAxisIsElementsTag(),)) - for out_grp, in_grp, vec_i, ijm_i in zip( - out_discr.groups, - in_discr.groups, - vec, - inv_jac_mat, strict=True)] + DOFArray(actx, data=tuple([ # noqa: C409 + pgg_i[xyz_axis] for pgg_i in per_group_grads + ])) + for xyz_axis in range(output_discr.ambient_dim) + ]) + + +def _divergence_kernel( + actx: ArrayContext, + output_discr: Discretization, + input_discr: Discretization, + inv_jac_mat: Iterable, + vec: DOFArray, + compute_single_axis_derivative: Callable, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: + """ + Computes divergence of *vec* by summing over each spatial derivative + computed by *compute_single_axis_derivative*. + + Strong and weak dispatch routines provide the proper + *compute_single_axis_derivative* routine. + """ + + per_group_divs = [] + for output_group, input_group, vec_i, ijm_i in zip( + output_discr.groups, input_discr.groups, vec, inv_jac_mat, + strict=False): + + per_group_divs.append(sum(compute_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec_ij, ijm_i, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) + for xyz_axis, vec_ij in enumerate(vec_i))) return DOFArray(actx, data=tuple(per_group_divs)) # }}} -# {{{ Derivative operators - -def _reference_derivative_matrices(actx: ArrayContext, - out_element_group: NodalElementGroupBase, - in_element_group: InterpolatoryElementGroupBase): - - @keyed_memoize_in( - actx, _reference_derivative_matrices, - lambda outgrp, ingrp: ( - outgrp.discretization_key(), - ingrp.discretization_key())) - def get_ref_derivative_mats( - out_grp: NodalElementGroupBase, - in_grp: InterpolatoryElementGroupBase): - return actx.freeze( - actx.tag_axis( - 1, DiscretizationDOFAxisTag(), - actx.from_numpy( - np.asarray( - mp.diff_matrices( - in_grp.basis_obj(), - out_grp.unit_nodes, - from_nodes=in_grp.unit_nodes, - ))))) - return get_ref_derivative_mats(out_element_group, in_element_group) - - -def _strong_scalar_grad(dcoll, dd_in, vec): - assert isinstance(dd_in.domain_tag, VolumeDomainTag) +# {{{ Strong derivative operators + +def _strong_tensor_product_single_axis_derivative( + actx: ArrayContext, + xyz_axis: int, + output_group: NodalElementGroupBase, + input_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer, + metrics: ArrayOrContainer + ) -> ArrayOrContainer: + + return sum(metrics[xyz_axis, rst_axis] * unfold( + output_group.space, + _single_axis_contraction( + actx, + input_group.dim, + rst_axis, + reference_derivative_matrices( + actx, output_group, input_group, + use_tensor_product_fast_eval=True)[0], + fold(input_group.space, vec), + arg_names=("diff_op_1d", "vec") + ) + ) + for rst_axis in range(input_group.dim) + ) - from grudge.geometry import inverse_surface_metric_derivative_mat - discr = dcoll.discr_from_dd(dd_in) - actx = vec.array_context +def _strong_simplicial_single_axis_derivative( + actx: ArrayContext, + xyz_axis: int, + output_group: NodalElementGroupBase, + input_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer, + metrics: ArrayOrContainer + ) -> ArrayOrContainer: + + return sum( + metrics[xyz_axis, rst_axis] * actx.einsum( + "ij,ej->ei", + reference_derivative_matrices( + actx, output_group, input_group, + use_tensor_product_fast_eval=False)[rst_axis], + vec, + arg_names=("diff_op", "vec")) + for rst_axis in range(input_group.dim) + ) - inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, - _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - return _gradient_kernel(actx, discr, discr, - _reference_derivative_matrices, inverse_jac_mat, vec, - metric_in_matvec=False) +def _strong_scalar_d_dx( + actx: ArrayContext, + xyz_axis: int, + output_group: NodalElementGroupBase, + input_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer, + metrics: ArrayOrContainer, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: + + if isinstance(input_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + return _strong_tensor_product_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec, metrics + ) + + return _strong_simplicial_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec, metrics + ) + + +def _strong_scalar_grad( + dcoll: DiscretizationCollection, + dd_in: DOFDesc, + vec: ArrayOrContainer, + *args, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: + + from grudge.geometry import inverse_surface_metric_derivative_mat + + assert isinstance(dd_in.domain_tag, VolumeDomainTag) -def _strong_scalar_div(dcoll, dd, vecs): - from arraycontext import get_container_context_recursively, serialize_container + discr = dcoll.discr_from_dd(dd_in) + actx = vec.array_context + metrics = inverse_surface_metric_derivative_mat(actx, dcoll, + dd=dd_in, _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + + return _gradient_kernel( + actx, discr, discr, metrics, vec, _strong_scalar_d_dx, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + +def _strong_scalar_div( + dcoll: DiscretizationCollection, + dd: DOFDesc, + vecs: ArrayOrContainer, + *args, + use_tensor_product_fast_eval: bool = True) -> ArrayOrContainer: + from arraycontext import ( + get_container_context_recursively, + serialize_container, + ) from grudge.geometry import inverse_surface_metric_derivative_mat assert isinstance(vecs, np.ndarray) assert vecs.shape == (dcoll.ambient_dim,) - discr = dcoll.discr_from_dd(dd) - actx = get_container_context_recursively(vecs) - vec = actx.np.stack([v for k, v in serialize_container(vecs)]) + assert actx is not None - inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd, + discr = dcoll.discr_from_dd(dd) + vec = actx.np.stack([v for _, v in serialize_container(vecs)]) + metrics = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd, _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - return _divergence_kernel(actx, discr, discr, - _reference_derivative_matrices, inverse_jac_mat, vec, - metric_in_matvec=False) + return _divergence_kernel( + actx, discr, discr, metrics, vec, _strong_scalar_d_dx, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) def local_grad( - dcoll: DiscretizationCollection, *args, nested=False) -> ArrayOrContainer: + dcoll: DiscretizationCollection, + *args, + nested=False, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: r"""Return the element-local gradient of a function :math:`f` represented by *vec*: @@ -348,13 +494,18 @@ def local_grad( from grudge.tools import rec_map_subarrays return rec_map_subarrays( - partial(_strong_scalar_grad, dcoll, dd_in), + partial(_strong_scalar_grad, dcoll, dd_in, + use_tensor_product_fast_eval=use_tensor_product_fast_eval), (), (dcoll.ambient_dim,), vec, scalar_cls=DOFArray, return_nested=nested,) def local_d_dx( - dcoll: DiscretizationCollection, xyz_axis, *args) -> ArrayOrContainer: + dcoll: DiscretizationCollection, + xyz_axis, + *args, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: r"""Return the element-local derivative along axis *xyz_axis* of a function :math:`f` represented by *vec*: @@ -382,22 +533,28 @@ def local_d_dx( raise TypeError("invalid number of arguments") if not isinstance(vec, DOFArray): - return map_array_container(partial(local_d_dx, dcoll, xyz_axis, dd), vec) + return map_array_container( + partial(local_d_dx, dcoll, xyz_axis, dd, + use_tensor_product_fast_eval=use_tensor_product_fast_eval), + vec) + + from grudge.geometry import inverse_surface_metric_derivative_mat discr = dcoll.discr_from_dd(dd) actx = vec.array_context - - from grudge.geometry import inverse_surface_metric_derivative_mat - inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd, + metrics = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd, _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) return _single_axis_derivative_kernel( - actx, discr, discr, - _reference_derivative_matrices, inverse_jac_mat, xyz_axis, vec, - metric_in_matvec=False) + actx, discr, discr, metrics, xyz_axis, vec, _strong_scalar_d_dx, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) -def local_div(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: +def local_div( + dcoll: DiscretizationCollection, + *args, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: r"""Return the element-local divergence of the vector function :math:`\mathbf{f}` represented by *vecs*: @@ -427,7 +584,9 @@ def local_div(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: from grudge.tools import rec_map_subarrays return rec_map_subarrays( - lambda vec: _strong_scalar_div(dcoll, dd, vec), + lambda vec: _strong_scalar_div( + dcoll, dd, vec, + use_tensor_product_fast_eval=use_tensor_product_fast_eval), (dcoll.ambient_dim,), (), vecs, scalar_cls=DOFArray) @@ -436,127 +595,199 @@ def local_div(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: # {{{ Weak derivative operators -def _reference_stiffness_transpose_matrices( - actx: ArrayContext, out_element_group, in_element_group): - @keyed_memoize_in( - actx, _reference_stiffness_transpose_matrices, - lambda out_grp, in_grp: (out_grp.discretization_key(), - in_grp.discretization_key())) - def get_ref_stiffness_transpose_mat(out_grp, in_grp): - if in_grp == out_grp: - mmat = mp.mass_matrix(out_grp.basis_obj(), out_grp.unit_nodes) - diff_matrices = mp.diff_matrices(out_grp.basis_obj(), out_grp.unit_nodes) - return actx.freeze( - actx.tag_axis(1, DiscretizationDOFAxisTag(), - actx.from_numpy( - np.asarray( - [dmat.T @ mmat.T for dmat in diff_matrices])))) - - from modepy import multi_vandermonde, vandermonde - - basis = out_grp.basis_obj() - vand = vandermonde(basis.functions, out_grp.unit_nodes) - grad_vand = multi_vandermonde(basis.gradients, in_grp.unit_nodes) - vand_inv_t = np.linalg.inv(vand).T - - if not isinstance(grad_vand, tuple): - # NOTE: special case for 1d - grad_vand = (grad_vand,) - - weights = in_grp.quadrature_rule().weights - return actx.freeze( - actx.from_numpy( - np.einsum( - "c,bz,acz->abc", - weights, - vand_inv_t, - grad_vand - ).copy() # contigify the array +def _weak_tensor_product_single_axis_derivative( + actx: ArrayContext, + xyz_axis: int, + output_group: NodalElementGroupBase, + input_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer, + metrics: ArrayOrContainer + ) -> ArrayOrContainer: + + vec_with_metrics = vec * metrics[xyz_axis] + weak_derivative = 0.0 + + stiff_t = reference_stiffness_transpose_matrices( + actx, input_group, output_group, + use_tensor_product_fast_eval=True + ) + + mass = reference_mass_matrix( + actx, input_group, output_group, + use_tensor_product_fast_eval=True + ) + + for rst_axis in range(input_group.dim): + apply_mass_axes = set(range(input_group.dim)) - {rst_axis} + weak_ref_derivative = fold(input_group.space, + vec_with_metrics[rst_axis]) + + for ax in apply_mass_axes: + weak_ref_derivative = _single_axis_contraction( + actx, + input_group.dim, + ax, + mass, + weak_ref_derivative, + arg_names=("mass_1d", "vec") ) + + weak_derivative += _single_axis_contraction( + actx, + input_group.dim, + rst_axis, + stiff_t, + weak_ref_derivative, + arg_names=("stiff_t_1d", "vec") + ) + + return tag_axes( + actx, + { + 0: DiscretizationElementAxisTag(), + 1: DiscretizationDOFAxisTag() + }, + unfold(output_group.space, weak_derivative) + ) + + +def _weak_simplicial_single_axis_derivative( + actx: ArrayContext, + xyz_axis: int, + output_group: NodalElementGroupBase, + input_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer, + metrics: ArrayOrContainer + ) -> ArrayOrContainer: + + vec_with_metrics = vec * metrics[xyz_axis] + + return sum( + actx.einsum( + "ij,ej->ei", + reference_stiffness_transpose_matrices( + actx, input_group, output_group, + use_tensor_product_fast_eval=False)[rst_axis], + vec_with_metrics[rst_axis], + arg_names=(f"stiffness_t_{rst_axis}", "vec_scaled")) + for rst_axis in range(input_group.dim) + ) + + +def _weak_scalar_d_dx( + actx: ArrayContext, + xyz_axis: int, + output_group: NodalElementGroupBase, + input_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer, + metrics: ArrayOrContainer, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: + + if input_group.dim == 1: + use_tensor_product_fast_eval = False + + if isinstance(input_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + return _weak_tensor_product_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec, metrics ) - return get_ref_stiffness_transpose_mat(out_element_group, - in_element_group) + return _weak_simplicial_single_axis_derivative( + actx, xyz_axis, output_group, input_group, vec, metrics + ) -def _weak_scalar_grad(dcoll, dd_in, vec): - from grudge.geometry import inverse_surface_metric_derivative_mat - dd_in = as_dofdesc(dd_in) - in_discr = dcoll.discr_from_dd(dd_in) - out_discr = dcoll.discr_from_dd(dd_in.with_discr_tag(DISCR_TAG_BASE)) +def _weak_scalar_grad(dcoll, dd_in, vec, *args, + use_tensor_product_fast_eval=True): + + # {{{ setup and grab discretizations actx = vec.array_context - inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, - times_area_element=True, - _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - return _gradient_kernel(actx, out_discr, in_discr, - _reference_stiffness_transpose_matrices, inverse_jac_mat, vec, - metric_in_matvec=True) + dd_in = as_dofdesc(dd_in) + input_discr = dcoll.discr_from_dd(dd_in) + output_discr = dcoll.discr_from_dd(dd_in.with_discr_tag(DISCR_TAG_BASE)) + # }}} -def _weak_scalar_div(dcoll, dd_in, vecs): - from arraycontext import get_container_context_recursively, serialize_container + # {{{ get metrics and apply scaling terms - from grudge.geometry import inverse_surface_metric_derivative_mat + metrics = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, + times_area_element=False, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - assert isinstance(vecs, np.ndarray) - assert vecs.shape == (dcoll.ambient_dim,) + # guaranteed "scalar" -> ensure axes are properly tagged + # WARNING: reshapes from TP fast operator eval can break propagation upward. + # hence, it is unwise to remove this. + vec_scaled = tag_axes( + actx, + { + 0: DiscretizationElementAxisTag(), + 1: DiscretizationDOFAxisTag() + }, + vec * area_element(actx, dcoll, dd=dd_in, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + ) + + # }}} + + return _gradient_kernel( + actx, output_discr, input_discr, metrics, vec_scaled, _weak_scalar_d_dx, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) - in_discr = dcoll.discr_from_dd(dd_in) - out_discr = dcoll.discr_from_dd(dd_in.with_discr_tag(DISCR_TAG_BASE)) + +def _weak_scalar_div(dcoll, dd_in, vecs, *args, + use_tensor_product_fast_eval=True): + from arraycontext import ( + get_container_context_recursively, + serialize_container, + ) + + # {{{ setup and grab discretizations actx = get_container_context_recursively(vecs) - vec = actx.np.stack([v for k, v in serialize_container(vecs)]) + assert actx is not None - inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, - times_area_element=True, - _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + vec = actx.np.stack([v for _, v in serialize_container(vecs)]) - return _divergence_kernel(actx, out_discr, in_discr, - _reference_stiffness_transpose_matrices, inverse_jac_mat, vec, - metric_in_matvec=True) + dd_in = as_dofdesc(dd_in) + input_discr = dcoll.discr_from_dd(dd_in) + output_discr = dcoll.discr_from_dd(dd_in.with_discr_tag(DISCR_TAG_BASE)) + # }}} -def weak_local_grad( - dcoll: DiscretizationCollection, *args, nested=False) -> ArrayOrContainer: - r"""Return the element-local weak gradient of the volume function - represented by *vec*. + # {{{ get metrics and apply scaling terms - May be called with ``(vec)`` or ``(dd_in, vec)``. + metrics = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, + times_area_element=False, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - Specifically, the function returns an object array where the :math:`i`-th - component is the weak derivative with respect to the :math:`i`-th coordinate - of a scalar function :math:`f`. See :func:`weak_local_d_dx` for further - information. For non-scalar :math:`f`, the function will return a nested object - array containing the component-wise weak derivatives. + # guaranteed "scalar" -> ensure axes are properly tagged + # WARNING: reshapes from TP fast operator eval can break propagation upward. + # hence, it is unwise to remove this. + vec_scaled = tag_axes( + actx, + { + 0: DiscretizationAmbientDimAxisTag(), + 1: DiscretizationElementAxisTag(), + 2: DiscretizationDOFAxisTag() + }, + vec * area_element(actx, dcoll, dd=dd_in, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + ) - :arg dd_in: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. - Defaults to the base volume discretization if not provided. - :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an - :class:`~arraycontext.ArrayContainer` of them. - :arg nested: return nested object arrays instead of a single multidimensional - array if *vec* is non-scalar - :returns: an object array (possibly nested) of - :class:`~meshmode.dof_array.DOFArray`\ s or - :class:`~arraycontext.ArrayContainer` of object arrays. - """ - if len(args) == 1: - vecs, = args - dd_in = DD_VOLUME_ALL - elif len(args) == 2: - dd_in, vecs = args - else: - raise TypeError("invalid number of arguments") + # }}} - from grudge.tools import rec_map_subarrays - return rec_map_subarrays( - partial(_weak_scalar_grad, dcoll, dd_in), - (), (dcoll.ambient_dim,), - vecs, scalar_cls=DOFArray, return_nested=nested) + return _divergence_kernel( + actx, output_discr, input_discr, metrics, vec_scaled, _weak_scalar_d_dx, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) -def weak_local_d_dx(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: +def weak_local_d_dx(dcoll: DiscretizationCollection, *args, + use_tensor_product_fast_eval=True) -> ArrayOrContainer: r"""Return the element-local weak derivative along axis *xyz_axis* of the volume function represented by *vec*. @@ -600,24 +831,69 @@ def weak_local_d_dx(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: vec ) - from grudge.geometry import inverse_surface_metric_derivative_mat - dd_in = as_dofdesc(dd_in) + in_discr = dcoll.discr_from_dd(dd_in) out_discr = dcoll.discr_from_dd(dd_in.with_discr_tag(DISCR_TAG_BASE)) actx = vec.array_context - inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, - times_area_element=True, - _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + metrics = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + vec_scaled = vec * area_element(actx, dcoll, dd=dd_in, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) return _single_axis_derivative_kernel( - actx, out_discr, in_discr, _reference_stiffness_transpose_matrices, - inverse_jac_mat, xyz_axis, vec, - metric_in_matvec=True) + actx, out_discr, in_discr, metrics, xyz_axis, vec_scaled, + _weak_scalar_d_dx, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + +def weak_local_grad( + dcoll: DiscretizationCollection, + *args, + nested=False, + use_tensor_product_fast_eval=True + ) -> ArrayOrContainer: + r"""Return the element-local weak gradient of the volume function + represented by *vec*. + + May be called with ``(vec)`` or ``(dd_in, vec)``. + + Specifically, the function returns an object array where the :math:`i`-th + component is the weak derivative with respect to the :math:`i`-th coordinate + of a scalar function :math:`f`. See :func:`weak_local_d_dx` for further + information. For non-scalar :math:`f`, the function will return a nested object + array containing the component-wise weak derivatives. + :arg dd_in: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. + Defaults to the base volume discretization if not provided. + :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an + :class:`~arraycontext.ArrayContainer` of them. + :arg nested: return nested object arrays instead of a single multidimensional + array if *vec* is non-scalar + :returns: an object array (possibly nested) of + :class:`~meshmode.dof_array.DOFArray`\ s or + :class:`~arraycontext.ArrayContainer` of object arrays. + """ + if len(args) == 1: + vecs, = args + dd_in = DD_VOLUME_ALL + elif len(args) == 2: + dd_in, vecs = args + else: + raise TypeError("invalid number of arguments") -def weak_local_div(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: + return rec_map_subarrays( + partial(_weak_scalar_grad, dcoll, dd_in, + use_tensor_product_fast_eval=use_tensor_product_fast_eval), + (), (dcoll.ambient_dim,), + vecs, scalar_cls=DOFArray, return_nested=nested, + ) + + +def weak_local_div(dcoll: DiscretizationCollection, + *args, + use_tensor_product_fast_eval=True) -> ArrayOrContainer: r"""Return the element-local weak divergence of the vector volume function represented by *vecs*. @@ -654,49 +930,64 @@ def weak_local_div(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: else: raise TypeError("invalid number of arguments") - from grudge.tools import rec_map_subarrays return rec_map_subarrays( - lambda vec: _weak_scalar_div(dcoll, dd_in, vec), - (dcoll.ambient_dim,), (), - vecs, scalar_cls=DOFArray) + lambda vec: _weak_scalar_div( + dcoll, dd_in, vec, + use_tensor_product_fast_eval=use_tensor_product_fast_eval), + (dcoll.ambient_dim,), (), vecs, scalar_cls=DOFArray + ) # }}} # {{{ Mass operator -def reference_mass_matrix(actx: ArrayContext, out_element_group, in_element_group): - @keyed_memoize_in( - actx, reference_mass_matrix, - lambda out_grp, in_grp: (out_grp.discretization_key(), - in_grp.discretization_key())) - def get_ref_mass_mat(out_grp, in_grp): - if out_grp == in_grp: - return actx.freeze( - actx.from_numpy( - mp.mass_matrix(out_grp.basis_obj(), out_grp.unit_nodes) - ) - ) +def _apply_mass_tensor_product( + actx: ArrayContext, + input_group: NodalElementGroupBase, + output_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer + ) -> ArrayOrContainer: + for xyz_axis in range(input_group.dim): + vec = _single_axis_contraction( + actx, + input_group.dim, + xyz_axis, + reference_mass_matrix( + actx, + output_group=output_group, + input_group=input_group, + use_tensor_product_fast_eval=True), + vec, + arg_names=("ref_mass_1d", "dofs")) - from modepy import vandermonde - basis = out_grp.basis_obj() - vand = vandermonde(basis.functions, out_grp.unit_nodes) - o_vand = vandermonde(basis.functions, in_grp.unit_nodes) - vand_inv_t = np.linalg.inv(vand).T + return vec - weights = in_grp.quadrature_rule().weights - return actx.freeze( - actx.tag_axis(0, DiscretizationDOFAxisTag(), - actx.from_numpy( - np.asarray( - np.einsum("j,ik,jk->ij", weights, vand_inv_t, o_vand), - order="C")))) - return get_ref_mass_mat(out_element_group, in_element_group) +def _apply_mass_simplicial( + actx: ArrayContext, + input_group: NodalElementGroupBase, + output_group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer + ) -> ArrayOrContainer: + return actx.einsum("ij,ej->ei", + reference_mass_matrix( + actx, + output_group=output_group, + input_group=input_group, + use_tensor_product_fast_eval=False), + vec, + arg_names=("mass_mat", "vec")) def _apply_mass_operator( - dcoll: DiscretizationCollection, dd_out, dd_in, vec): + dcoll: DiscretizationCollection, + dd_out: DOFDesc, + dd_in: DOFDesc, + vec: ArrayOrContainer, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: + if not isinstance(vec, DOFArray): return map_array_container( partial(_apply_mass_operator, dcoll, dd_out, dd_in), vec @@ -704,34 +995,43 @@ def _apply_mass_operator( from grudge.geometry import area_element - in_discr = dcoll.discr_from_dd(dd_in) - out_discr = dcoll.discr_from_dd(dd_out) + input_discr = dcoll.discr_from_dd(dd_in) + output_discr = dcoll.discr_from_dd(dd_out) actx = vec.array_context - area_elements = area_element(actx, dcoll, dd=dd_in, + vec = vec * area_element(actx, dcoll, dd=dd_in, _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - return DOFArray( - actx, - data=tuple( - actx.einsum("ij,ej,ej->ei", - reference_mass_matrix( + + group_data = [] + for input_group, output_group, vec_i in zip( + input_discr.groups, output_discr.groups, vec, strict=False): + if isinstance(input_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval and input_group.dim != 1: + group_data.append( + tag_axes( actx, - out_element_group=out_grp, - in_element_group=in_grp - ), - ae_i, - vec_i, - arg_names=("mass_mat", "jac", "vec"), - tagged=(FirstAxisIsElementsTag(),)) - - for in_grp, out_grp, ae_i, vec_i in zip( - in_discr.groups, out_discr.groups, area_elements, vec, - strict=True) - ) - ) + { + 0: DiscretizationElementAxisTag(), + 1: DiscretizationDOFAxisTag() + }, + unfold(output_group.space, + _apply_mass_tensor_product( + actx, + input_group, + output_group, + fold(input_group.space, vec_i)))) + ) + + else: + group_data.append(_apply_mass_simplicial( + actx, input_group, output_group, vec_i)) + + return DOFArray(actx, data=tuple(group_data)) -def mass(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: +def mass(dcoll: DiscretizationCollection, + *args, + use_tensor_product_fast_eval: bool = True) -> ArrayOrContainer: r"""Return the action of the DG mass matrix on a vector (or vectors) of :class:`~meshmode.dof_array.DOFArray`\ s, *vec*. In the case of *vec* being an :class:`~arraycontext.ArrayContainer`, @@ -767,303 +1067,382 @@ def mass(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: dd_out = dd_in.with_discr_tag(DISCR_TAG_BASE) - return _apply_mass_operator(dcoll, dd_out, dd_in, vec) + return _apply_mass_operator( + dcoll, dd_out, dd_in, vec, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) # }}} -# {{{ Mass inverse operator +# {{{ Face mass operator + +def _apply_face_mass_tensor_product( + actx: ArrayContext, + face_group: ElementGroupBase, + volume_group: NodalElementGroupBase, + vec: ArrayOrContainer, + ) -> ArrayOrContainer: -def reference_inverse_mass_matrix(actx: ArrayContext, element_group): - @keyed_memoize_in( - actx, reference_inverse_mass_matrix, - lambda grp: grp.discretization_key()) - def get_ref_inv_mass_mat(grp): - from modepy import inverse_mass_matrix - basis = grp.basis_obj() + vec = fold(face_group.space, vec) - return actx.freeze( - actx.tag_axis(0, DiscretizationDOFAxisTag(), - actx.from_numpy( - np.asarray( - inverse_mass_matrix(basis, grp.unit_nodes), - order="C")))) + face_mass_mat = reference_face_mass_matrix( + actx, + face_group=face_group, + vol_group=volume_group, + dtype=vec.dtype, + use_tensor_product_fast_eval=True + ) - return get_ref_inv_mass_mat(element_group) + return unfold(volume_group.space, vec) + + +def _apply_face_mass_simplicial( + actx: ArrayContext, + face_group: ElementGroupBase, + volume_group: NodalElementGroupBase, + vec: ArrayOrContainer + ) -> ArrayOrContainer: + + return actx.einsum( + "ifj,fej->ei", + reference_face_mass_matrix( + actx, + face_group=face_group, + vol_group=volume_group, + dtype=vec.dtype, + use_tensor_product_fast_eval=False + ), + vec.reshape( + volume_group.mesh_el_group.nfaces, + volume_group.nelements, + face_group.nunit_dofs), + arg_names=("ref_face_mass_mat", "vec") + ) -def _apply_inverse_mass_operator( - dcoll: DiscretizationCollection, dd_out, dd_in, vec): +def _apply_face_mass_operator( + dcoll: DiscretizationCollection, + dd_in: DOFDesc, + vec: ArrayOrContainer, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: if not isinstance(vec, DOFArray): return map_array_container( - partial(_apply_inverse_mass_operator, dcoll, dd_out, dd_in), vec + partial(_apply_face_mass_operator, dcoll, dd_in), vec ) from grudge.geometry import area_element - if dd_out != dd_in: - raise ValueError( - "Cannot compute inverse of a mass matrix mapping " - "between different element groups; inverse is not " - "guaranteed to be well-defined" - ) + dd_out = DOFDesc( + VolumeDomainTag(dd_in.domain_tag.volume_tag), + DISCR_TAG_BASE) + volm_discr = dcoll.discr_from_dd(dd_out) + face_discr = dcoll.discr_from_dd(dd_in) actx = vec.array_context - discr = dcoll.discr_from_dd(dd_in) - inv_area_elements = 1./area_element(actx, dcoll, dd=dd_in, + + assert len(face_discr.groups) == len(volm_discr.groups) + surf_area_elements = area_element(actx, dcoll, dd=dd_in, _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - group_data = [ - # Based on https://arxiv.org/pdf/1608.03836.pdf - # true_Minv ~ ref_Minv * ref_M * (1/jac_det) * ref_Minv - actx.einsum("ei,ij,ej->ei", - jac_inv, - reference_inverse_mass_matrix(actx, element_group=grp), - vec_i, - tagged=(FirstAxisIsElementsTag(),)) - for grp, jac_inv, vec_i in zip( - discr.groups, inv_area_elements, vec, strict=True)] + + group_data = [] + for vgroup, afgroup, vec_i, surf_ae_i in zip( + volm_discr.groups, face_discr.groups, vec, surf_area_elements, + strict=True): + + # fast evaluation only applicable to faces of 3D (or higher) elements + use_fast_eval = (use_tensor_product_fast_eval and volm_discr.dim > 3) + if isinstance(vgroup, TensorProductElementGroupBase) and use_fast_eval: + group_data.append( + _apply_face_mass_tensor_product( + actx, afgroup, vgroup, vec_i * surf_ae_i + ) + ) + + else: + group_data.append( + _apply_face_mass_simplicial( + actx, afgroup, vgroup, vec_i * surf_ae_i + ) + ) return DOFArray(actx, data=tuple(group_data)) -def inverse_mass(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: - r"""Return the action of the DG mass matrix inverse on a vector - (or vectors) of :class:`~meshmode.dof_array.DOFArray`\ s, *vec*. - In the case of *vec* being an :class:`~arraycontext.ArrayContainer`, - the inverse mass operator is applied component-wise. +def face_mass( + dcoll: DiscretizationCollection, + *args, + use_tensor_product_fast_eval: bool = True + ) -> ArrayOrContainer: + r"""Return the action of the DG face mass matrix on a vector (or vectors) + of :class:`~meshmode.dof_array.DOFArray`\ s, *vec*. In the case of + *vec* being an arbitrary :class:`~arraycontext.ArrayContainer`, + the face mass operator is applied component-wise. - For affine elements :math:`E`, the element-wise mass inverse - is computed directly as the inverse of the (physical) mass matrix: + May be called with ``(vec)`` or ``(dd_in, vec)``. - .. math:: + Specifically, this function applies the face mass matrix elementwise on a + vector of coefficients :math:`\mathbf{f}` as the sum of contributions for + each face :math:`f \subset \partial E`: - \left(\mathbf{M}_{J^e}\right)_{ij} = - \int_{\widehat{E}} \widehat{\phi}_i\cdot\widehat{\phi}_j J^e - \mathrm{d}\widehat{x}, + .. math:: - where :math:`\widehat{\phi}_i` are basis functions over the reference - element :math:`\widehat{E}`, and :math:`J^e` is the (constant) Jacobian - scaling factor (see :func:`grudge.geometry.area_element`). + \sum_{f=1}^{N_{\text{faces}} } \mathbf{M}_{f, E}\mathbf{f}|_f, - For non-affine :math:`E`, :math:`J^e` is not constant. In this case, a - weight-adjusted approximation is used instead following [Chan_2016]_: + where .. math:: - \mathbf{M}_{J^e}^{-1} \approx - \widehat{\mathbf{M}}^{-1}\mathbf{M}_{1/J^e}\widehat{\mathbf{M}}^{-1}, - - where :math:`\widehat{\mathbf{M}}` is the reference mass matrix on - :math:`\widehat{E}`. + \left(\mathbf{M}_{f, E}\right)_{ij} = + \int_{f \subset \partial E} \phi_i(s)\psi_j(s)\,\mathrm{d}s, - May be called with ``(vec)`` or ``(dd, vec)``. + where :math:`\phi_i` are (volume) polynomial basis functions on :math:`E` + evaluated on the face :math:`f`, and :math:`\psi_j` are basis functions for + a polynomial space defined on :math:`f`. + :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. + Defaults to the base ``"all_faces"`` discretization if not provided. :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an :class:`~arraycontext.ArrayContainer` of them. - :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. - Defaults to the base volume discretization if not provided. :returns: a :class:`~meshmode.dof_array.DOFArray` or an :class:`~arraycontext.ArrayContainer` like *vec*. """ + if len(args) == 1: vec, = args - dd = DD_VOLUME_ALL + dd_in = DD_VOLUME_ALL.trace(FACE_RESTR_ALL) elif len(args) == 2: - dd, vec = args + dd_in, vec = args else: raise TypeError("invalid number of arguments") - return _apply_inverse_mass_operator(dcoll, dd, dd, vec) + return _apply_face_mass_operator( + dcoll, dd_in, vec, + use_tensor_product_fast_eval=use_tensor_product_fast_eval) # }}} -# {{{ Face mass operator +# {{{ Mass inverse operator -def reference_face_mass_matrix( - actx: ArrayContext, face_element_group, vol_element_group, dtype): - @keyed_memoize_in( - actx, reference_mass_matrix, - lambda face_grp, vol_grp: (face_grp.discretization_key(), - vol_grp.discretization_key())) - def get_ref_face_mass_mat(face_grp, vol_grp): - nfaces = vol_grp.mesh_el_group.nfaces - assert face_grp.nelements == nfaces * vol_grp.nelements - - matrix = np.empty( - (vol_grp.nunit_dofs, - nfaces, - face_grp.nunit_dofs), - dtype=dtype - ) +def _apply_inverse_mass_tensor_product( + actx: ArrayContext, + group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer + ) -> ArrayOrContainer: - import modepy as mp - from meshmode.discretization import ElementGroupWithBasis - from meshmode.discretization.poly_element import QuadratureSimplexElementGroup - - n = vol_grp.order - m = face_grp.order - vol_basis = vol_grp.basis_obj() - faces = mp.faces_for_shape(vol_grp.shape) - - for iface, face in enumerate(faces): - # If the face group is defined on a higher-order - # quadrature grid, use the underlying quadrature rule - if isinstance(face_grp, QuadratureSimplexElementGroup): - face_quadrature = face_grp.quadrature_rule() - if face_quadrature.exact_to < m: - raise ValueError( - "The face quadrature rule is only exact for polynomials " - f"of total degree {face_quadrature.exact_to}. Please " - "ensure a quadrature rule is used that is at least " - f"exact for degree {m}." - ) - else: - # NOTE: This handles the general case where - # volume and surface quadrature rules may have different - # integration orders - face_quadrature = mp.quadrature_for_space( - mp.space_for_shape(face, 2*max(n, m)), - face - ) + vec = fold(group.space, vec) - # If the group has a nodal basis and is unisolvent, - # we use the basis on the face to compute the face mass matrix - if (isinstance(face_grp, ElementGroupWithBasis) - and face_grp.space.space_dim == face_grp.nunit_dofs): - - face_basis = face_grp.basis_obj() - - # Sanity check for face quadrature accuracy. Not integrating - # degree N + M polynomials here is asking for a bad time. - if face_quadrature.exact_to < m + n: - raise ValueError( - "The face quadrature rule is only exact for polynomials " - f"of total degree {face_quadrature.exact_to}. Please " - "ensure a quadrature rule is used that is at least " - f"exact for degree {n+m}." - ) - - matrix[:, iface, :] = mp.nodal_mass_matrix_for_face( - face, face_quadrature, - face_basis.functions, vol_basis.functions, - vol_grp.unit_nodes, - face_grp.unit_nodes, - ) - else: - # Otherwise, we use a routine that is purely quadrature-based - # (no need for explicit face basis functions) - matrix[:, iface, :] = mp.nodal_quad_mass_matrix_for_face( - face, - face_quadrature, - vol_basis.functions, - vol_grp.unit_nodes, - ) + for rst_axis in range(group.dim): + vec = _single_axis_contraction( + actx, group.dim, rst_axis, + reference_inverse_mass_matrix( + actx, group, + use_tensor_product_fast_eval=True), + vec, + arg_names=("base_mass_inv", "dofs")) - return actx.freeze( - tag_axes(actx, { - 0: DiscretizationDOFAxisTag(), - 2: DiscretizationDOFAxisTag() - }, - actx.from_numpy(matrix))) + return unfold(group.space, vec) - return get_ref_face_mass_mat(face_element_group, vol_element_group) +def _apply_inverse_mass_simplicial( + actx: ArrayContext, + group: InterpolatoryElementGroupBase, + vec: ArrayOrContainer + ) -> ArrayOrContainer: -def _apply_face_mass_operator(dcoll: DiscretizationCollection, dd_in, vec): + return actx.einsum( + "ij,ej->ei", + reference_inverse_mass_matrix( + actx, group, + use_tensor_product_fast_eval=False), + vec, + arg_names=("ref_inv_mass", "vec") + ) + + +# {{{ non-WADG + +def _apply_quad_inverse_mass( + actx: ArrayContext, + discr: Discretization, + vec: ArrayOrContainer, + inv_area_elements: ArrayOrContainer, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: + + group_data = [] + for input_group, output_group, vec_i, inv_ae_i in zip( + discr.groups, discr.groups, vec, inv_area_elements, strict=False): + if isinstance(input_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval and input_group.dim != 1: + group_data.append(_apply_inverse_mass_tensor_product( + actx, output_group, vec_i) * inv_ae_i) + else: + group_data.append(_apply_inverse_mass_simplicial( + actx, output_group, vec_i) * inv_ae_i) + + return DOFArray(actx, data=tuple(group_data)) + +# }}} + + +# {{{ WADG + +def _apply_wadg_inverse_mass( + actx: ArrayContext, + dcoll: DiscretizationCollection, + dd_out: DOFDesc, + dd_in: DOFDesc, + vec: ArrayOrContainer, + inv_area_elements: ArrayOrContainer, + use_tensor_product_fast_eval: bool = True + ) -> DOFArray: + + discr_base = dcoll.discr_from_dd(dd_out) + discr_quad = dcoll.discr_from_dd(dd_in) + + group_data = [] + for input_group, output_group, vec_i in zip( + discr_quad.groups, discr_base.groups, vec, strict=False): + + if isinstance(input_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + group_data.append(_apply_inverse_mass_tensor_product( + actx, output_group, vec_i)) + + else: + group_data.append(_apply_inverse_mass_simplicial( + actx, output_group, vec_i)) + + vec = inv_area_elements * project( + dcoll, dd_out, dd_in, DOFArray(actx, data=tuple(group_data))) + + group_data = [] + for input_group, output_group, vec_i in zip( + discr_quad.groups, discr_base.groups, vec, strict=False): + + if isinstance(input_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + group_data.append(_apply_inverse_mass_tensor_product( + actx, output_group, _apply_mass_tensor_product( + actx, input_group, output_group, vec_i))) + + else: + group_data.append(_apply_inverse_mass_simplicial( + actx, output_group, _apply_mass_simplicial( + actx, input_group, output_group, vec_i))) + + return DOFArray(actx, data=tuple(group_data)) + +# }}} + + +def _apply_inverse_mass( + dcoll: DiscretizationCollection, + dd_out: DOFDesc, + dd_in: DOFDesc, + vec: ArrayOrContainer, + uses_quadrature: bool = False, + use_tensor_product_fast_eval: bool = True): if not isinstance(vec, DOFArray): return map_array_container( - partial(_apply_face_mass_operator, dcoll, dd_in), vec + partial(_apply_inverse_mass, dcoll, dd_out, dd_in, + uses_quadrature, use_tensor_product_fast_eval), vec ) - from grudge.geometry import area_element - - dd_out = DOFDesc( - VolumeDomainTag(dd_in.domain_tag.volume_tag), - DISCR_TAG_BASE) + if dd_out != dd_in and not uses_quadrature: + raise ValueError( + "Cannot compute inverse of a mass matrix mapping " + "between different element groups unless using overintegration; " + "inverse is not guaranteed to be well-defined") - volm_discr = dcoll.discr_from_dd(dd_out) - face_discr = dcoll.discr_from_dd(dd_in) - dtype = vec.entry_dtype actx = vec.array_context + discr_base = dcoll.discr_from_dd(dd_out) + discr_quad = dcoll.discr_from_dd(dd_in) - assert len(face_discr.groups) == len(volm_discr.groups) - surf_area_elements = area_element(actx, dcoll, dd=dd_in, - _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) + inv_area_elements = 1./area_element(actx, dcoll, dd=dd_in, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - return DOFArray( - actx, - data=tuple( - actx.einsum("ifj,fej,fej->ei", - reference_face_mass_matrix( - actx, - face_element_group=afgrp, - vol_element_group=vgrp, - dtype=dtype), - actx.tag_axis(1, DiscretizationElementAxisTag(), - surf_ae_i.reshape( - vgrp.mesh_el_group.nfaces, - vgrp.nelements, - surf_ae_i.shape[-1])), - actx.tag_axis(0, DiscretizationFaceAxisTag(), - vec_i.reshape( - vgrp.mesh_el_group.nfaces, - vgrp.nelements, - afgrp.nunit_dofs)), - arg_names=("ref_face_mass_mat", "jac_surf", "vec"), - tagged=(FirstAxisIsElementsTag(),)) - - for vgrp, afgrp, vec_i, surf_ae_i in zip(volm_discr.groups, - face_discr.groups, - vec, - surf_area_elements, - strict=True))) - - -def face_mass(dcoll: DiscretizationCollection, *args) -> ArrayOrContainer: - r"""Return the action of the DG face mass matrix on a vector (or vectors) - of :class:`~meshmode.dof_array.DOFArray`\ s, *vec*. In the case of - *vec* being an arbitrary :class:`~arraycontext.ArrayContainer`, - the face mass operator is applied component-wise. + if discr_base == discr_quad: + return _apply_quad_inverse_mass( + actx, discr_base, vec, inv_area_elements, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) - May be called with ``(vec)`` or ``(dd_in, vec)``. + else: + if use_tensor_product_fast_eval: + from warnings import warn + warn("Fast operator evaluation + overintegration is not supported. " + "Defaulting to typical (full) operator evaluation.", + stacklevel=1) + + # see WADG: https://arxiv.org/pdf/1608.03836 + return _apply_wadg_inverse_mass( + actx, dcoll, dd_out, dd_in, vec, inv_area_elements, + use_tensor_product_fast_eval=False + ) - Specifically, this function applies the face mass matrix elementwise on a - vector of coefficients :math:`\mathbf{f}` as the sum of contributions for - each face :math:`f \subset \partial E`: + +def inverse_mass(dcoll: DiscretizationCollection, *args, + use_tensor_product_fast_eval=True) -> ArrayOrContainer: + r"""Return the action of the DG mass matrix inverse on a vector + (or vectors) of :class:`~meshmode.dof_array.DOFArray`\ s, *vec*. + In the case of *vec* being an :class:`~arraycontext.ArrayContainer`, + the inverse mass operator is applied component-wise. + + For affine elements :math:`E`, the element-wise mass inverse + is computed directly as the inverse of the (physical) mass matrix: .. math:: - \sum_{f=1}^{N_{\text{faces}} } \mathbf{M}_{f, E}\mathbf{f}|_f, + \left(\mathbf{M}_{J^e}\right)_{ij} = + \int_{\widehat{E}} \widehat{\phi}_i\cdot\widehat{\phi}_j J^e + \mathrm{d}\widehat{x}, - where + where :math:`\widehat{\phi}_i` are basis functions over the reference + element :math:`\widehat{E}`, and :math:`J^e` is the (constant) Jacobian + scaling factor (see :func:`grudge.geometry.area_element`). + + For non-affine :math:`E`, :math:`J^e` is not constant. In this case, a + weight-adjusted approximation is used instead following [Chan_2016]_: .. math:: - \left(\mathbf{M}_{f, E}\right)_{ij} = - \int_{f \subset \partial E} \phi_i(s)\psi_j(s)\,\mathrm{d}s, + \mathbf{M}_{J^e}^{-1} \approx + \widehat{\mathbf{M}}^{-1}\mathbf{M}_{1/J^e}\widehat{\mathbf{M}}^{-1}, - where :math:`\phi_i` are (volume) polynomial basis functions on :math:`E` - evaluated on the face :math:`f`, and :math:`\psi_j` are basis functions for - a polynomial space defined on :math:`f`. + where :math:`\widehat{\mathbf{M}}` is the reference mass matrix on :math:`\widehat{E}`. + + May be called with ``(vec)`` or ``(dd, vec)``. - :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. - Defaults to the base ``"all_faces"`` discretization if not provided. :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an :class:`~arraycontext.ArrayContainer` of them. + :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. + Defaults to the base volume discretization if not provided. :returns: a :class:`~meshmode.dof_array.DOFArray` or an :class:`~arraycontext.ArrayContainer` like *vec*. """ - if len(args) == 1: vec, = args - dd_in = DD_VOLUME_ALL.trace(FACE_RESTR_ALL) + dd = DD_VOLUME_ALL elif len(args) == 2: - dd_in, vec = args + dd, vec = args else: raise TypeError("invalid number of arguments") - return _apply_face_mass_operator(dcoll, dd_in, vec) + if dd.uses_quadrature(): + dd_in = dd.with_discr_tag(DISCR_TAG_QUAD) + dd_out = dd.with_discr_tag(DISCR_TAG_BASE) + else: + dd_in = dd + dd_out = dd + + return _apply_inverse_mass( + dcoll, dd_out, dd_in, vec, + uses_quadrature=dd.uses_quadrature(), + use_tensor_product_fast_eval=use_tensor_product_fast_eval) # }}} diff --git a/grudge/tools.py b/grudge/tools.py index 69405d62d..8dc89031f 100644 --- a/grudge/tools.py +++ b/grudge/tools.py @@ -36,6 +36,21 @@ import numpy as np from arraycontext import ArrayContext, ArrayOrContainer, ArrayOrContainerT +from meshmode.discretization.poly_element import ( + QuadratureSimplexElementGroup, + TensorProductElementGroupBase +) +from meshmode.discretization import ( + ElementGroupBase, + ElementGroupWithBasis, + NodalElementGroupBase, +) +from modepy import ( + Basis, + Face, + Quadrature, + TensorProductQuadrature +) from pytools import product @@ -79,6 +94,148 @@ def build_jacobian( # }}} +# {{{ discretization info extraction helpers + +def get_element_group_basis( + group: ElementGroupWithBasis, + use_tensor_product_fast_eval: bool = True + ) -> Basis: + + if isinstance(group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + + basis_obj = group.basis_obj().bases[0] + + if not sum(b == basis_obj for b in group.basis_obj().bases): + raise ValueError( + "Fast operator evaluation for tensor-product " + "discretizations requires that only a single " + "basis is used to construct the tensor-product") + + else: + basis_obj = group.basis_obj() + + return basis_obj + + +def get_element_group_nodes( + group: NodalElementGroupBase, + use_tensor_product_fast_eval: bool = True + ) -> np.ndarray: + + if isinstance(group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + return group.unit_nodes_1d + return group.unit_nodes + + +def get_accurate_quadrature_rule( + group: NodalElementGroupBase, + required_exactness: int | None = None, + use_tensor_product_fast_eval: bool = True, + ) -> Quadrature: + + import modepy as mp + + if not isinstance(group.quadrature_rule().exact_to, int): + return group.quadrature_rule() + + if required_exactness is None: + required_exactness = 2*group.order + + if group.quadrature_rule().exact_to < required_exactness: + quadrature_rule = mp.quadrature_for_space( + mp.space_for_shape(group.shape, required_exactness), + group.shape + ) + + else: + quadrature_rule = group.quadrature_rule() + + if isinstance(quadrature_rule, TensorProductQuadrature) and \ + use_tensor_product_fast_eval: + return quadrature_rule.quadratures[0] + + return quadrature_rule + + +def get_faces_for_volume_group( + volume_group: ElementGroupBase, + use_tensor_product_fast_eval: bool = True + ) -> tuple[Face, ...]: + + import modepy as mp + + if use_tensor_product_fast_eval and volume_group.dim > 2: + assert isinstance(volume_group, TensorProductElementGroupBase) + return mp.faces_for_shape(mp.Hypercube(volume_group.dim - 1)) + + return mp.faces_for_shape(volume_group.shape) + + +def get_quadrature_for_face( + face_group: NodalElementGroupBase, + face: Face, + required_exactness: int | None = None, + vol_group_order: int | None = None, + use_tensor_product_fast_eval: bool = True, + ) -> Quadrature: + + import modepy as mp + + if isinstance(face_group, QuadratureSimplexElementGroup): + if face_group.quadrature_rule().exact_to < face_group.order: + raise ValueError( + "The face quadrature rule is only exact for polynomials " + f"of total degree {face_group.quadrature_rule().exact_to}. " + "Please ensure a quadrature rule is used that is at least " + f"exact for degree {face_group.order}." + ) + + return face_group.quadrature_rule() + + if required_exactness is None: + if vol_group_order is None: + raise ValueError("Must supply one of `required_exactness` or " + "`vol_group_order` to construct a quadrature rule " + "accurate enough to evaluate the face mass matrix") + required_exactness = 2*max(face_group.order, vol_group_order) + + if not isinstance(face_group.quadrature_rule().exact_to, int): + return face_group.quadrature_rule() + + if face_group.quadrature_rule().exact_to < required_exactness: + quadrature_rule = mp.quadrature_for_space( + mp.space_for_shape(face, required_exactness), + face + ) + + else: + quadrature_rule = face_group.quadrature_rule() + + if isinstance(quadrature_rule, TensorProductQuadrature) and \ + use_tensor_product_fast_eval: + return quadrature_rule.quadratures[0] + + return quadrature_rule + + +def get_basis_for_face_group( + face_group: ElementGroupBase, + use_tensor_product_fast_eval: bool = True +) -> Basis | None: + + if isinstance(face_group, ElementGroupWithBasis): + return get_element_group_basis( + face_group, + use_tensor_product_fast_eval=use_tensor_product_fast_eval + ) + + return None + +# }}} + + # {{{ map_subarrays def map_subarrays( diff --git a/grudge/transform/mappers.py b/grudge/transform/mappers.py new file mode 100644 index 000000000..7b18c9aa6 --- /dev/null +++ b/grudge/transform/mappers.py @@ -0,0 +1,258 @@ +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2024 Addison Alvey-Blanco +Copyright (C) 2024 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. +""" + + +from meshmode.discretization import DiscretizationDOFAxisTag, DiscretizationElementAxisTag +import pytato as pt +from pytato import Array, Axis, ReductionDescriptor, Reshape +from pytato.analysis import is_einsum_similar_to_subscript +from pytato.array import Einsum, EinsumAxisDescriptor, EinsumElementwiseAxis, EinsumReductionAxis, immutabledict +from pytato.transform import ( + CombineMapper, + CopyMapper, + CopyMapperWithExtraArgs, +) + +from grudge.transform.metadata import ( + TensorProductDOFAxisTag, + TensorProductMassOperatorInverseTag, + TensorProductMassOperatorTag, + TensorProductOperatorAxisTag, + TensorProductStiffnessOperatorTag, +) + + +class MassCounter(CombineMapper): + def combine(self, *n_list): + return sum(n_list) + + def map_einsum(self, expr): + acc = 0 + for arg in expr.args: + if arg.tags_of_type(TensorProductMassOperatorTag): + acc += 1 + acc += self.rec(arg) + + return acc + + +class MassInverseCounter(CombineMapper): + def combine(self, *n_list): + return sum(n_list) + + def map_einsum(self, expr): + acc = 0 + for arg in expr.args: + if arg.tags_of_type(TensorProductMassOperatorInverseTag): + acc += 1 + acc += self.rec(arg) + + return acc + + +class InverseMassAttacher(CopyMapperWithExtraArgs): + def map_einsum( + self, + expr: Einsum, + inv_mass: Array, + access_descr: tuple[EinsumAxisDescriptor] + ) -> Array: + + new_args = [] + for iarg, arg in enumerate(expr.args): + if arg.tags_of_type(TensorProductMassOperatorTag): + if expr.access_descriptors[iarg] == access_descr: + a = inv_mass @ arg + a = a.copy(axes=tuple( + ax.tagged(TensorProductOperatorAxisTag()) + for ax in a.axes + )) + new_args.append(a) + continue + + elif arg.tags_of_type(TensorProductStiffnessOperatorTag): + if expr.access_descriptors[iarg] == access_descr: + a = inv_mass @ arg + a = a.copy(axes=tuple( + ax.tagged(TensorProductOperatorAxisTag(),) + for ax in a.axes + )) + new_args.append(a) + continue + + new_args.append(self.rec(arg, inv_mass, access_descr)) + + return expr.copy(args=tuple(new_args)) + + def map_reshape( + self, + expr: Reshape, + inv_mass: Array, + access_descr: tuple[EinsumAxisDescriptor] + ) -> Array: + + if isinstance(expr.array, Einsum): + if is_einsum_similar_to_subscript(expr.array, "ifj,fej->ei"): + _, nfaces, _ = expr.array.args[0].shape + + dim = int(nfaces / 2) + + output_ax = tuple( + descr for descr in access_descr + if isinstance(descr, EinsumElementwiseAxis) + ) + assert len(output_ax) == 1 + output_ax, = output_ax + + data_access_descrs = [] + for i in range(dim+1): + if i != output_ax.dim: + data_access_descrs.append(EinsumElementwiseAxis(i)) + else: + data_access_descrs.append(EinsumReductionAxis(0)) + access_descriptors = (access_descr, tuple(data_access_descrs)) + + redn_axis_to_redn_descr = immutabledict({ + EinsumReductionAxis(0): + ReductionDescriptor(tags=frozenset()) + }) + + axes = tuple( + Axis(tags=frozenset((DiscretizationElementAxisTag(),))) + if i == 0 + else Axis(tags=frozenset((TensorProductDOFAxisTag(i-1),))) + for i in range(dim+1) + ) + + return Einsum( + access_descriptors=access_descriptors, + args=(inv_mass, expr), + axes=axes, + redn_axis_to_redn_descr=redn_axis_to_redn_descr, + tags=frozenset() + ) + + return expr.copy(array=self.rec(expr.array, inv_mass, access_descr)) + + +class InverseMassDistributor(CopyMapper): + def map_einsum(self, expr: Einsum) -> Array: + for iarg, arg in enumerate(expr.args): + if not arg.tags_of_type(TensorProductMassOperatorInverseTag): + iarg_rec = iarg + break + + new_args = [] + for iarg, arg in enumerate(expr.args): + if arg.tags_of_type(TensorProductMassOperatorInverseTag): + return InverseMassAttacher()( + self.rec(expr.args[iarg_rec]), + arg, + expr.access_descriptors[iarg] + ) + + else: + new_args.append(self.rec(arg)) + + return expr.copy(args=tuple(new_args)) + + +def check_redundant_mass(expr: Einsum) -> bool: + found_mass = False + found_inverse_mass = False + + for arg in expr.args: + if arg.tags_of_type(TensorProductMassOperatorInverseTag): + found_inverse_mass = True + elif arg.tags_of_type(TensorProductMassOperatorTag): + found_mass = True + + return (found_inverse_mass and found_mass) + + +class RedundantMassRemover(CopyMapper): + def map_einsum(self, expr: Einsum) -> Array: + new_args = [] + for arg in expr.args: + if isinstance(arg, Einsum): + if check_redundant_mass(arg): + continue + new_args.append(self.rec(arg)) + + if len(new_args) == 1: + return self.rec(new_args[0]) + + return expr.copy(args=tuple(new_args)) + + +class FaceMassResultReshaper(CopyMapper): + def map_einsum(self, expr: Einsum) -> Array: + + new_args = [self.rec(arg) for arg in expr.args] + new_expr = expr.copy(args=tuple(new_args)) + + if is_einsum_similar_to_subscript(new_expr, "ifj,fej->ei"): + nfaces, nelts, _ = new_expr.args[1].shape + ndofs = expr.shape[1] + + from math import ceil + dim = ceil(nfaces / 2) + ndofs_1d = ceil(ndofs**(1/dim)) + + new_expr = new_expr.reshape(nelts, *(ndofs_1d,)*dim, order="F") + for i in range(dim+1): + new_expr = new_expr.with_tagged_axis( + i, ( + DiscretizationElementAxisTag() if i == 0 else + TensorProductDOFAxisTag(i-1) + ) + ) + + new_expr = new_expr.reshape(nelts, ndofs_1d**dim, order="F") + new_expr = new_expr.with_tagged_axis(0, + DiscretizationElementAxisTag()) + new_expr = new_expr.with_tagged_axis(1, DiscretizationDOFAxisTag()) + + return new_expr + + +def tensor_product_algebraic_transforms(dag): + # 0. preprocess face mass result (reshape to tp -> reshape from tp) + dag = FaceMassResultReshaper()(dag) + + # 1. distribute the inverse mass to: + # - einsums with stiffness + # - einsums with mass + # - face mass (insert applications between reshapes) + dag = InverseMassDistributor()(dag) + + # 2. remove einsums with mass and mass inverse + dag = RedundantMassRemover()(dag) + + # done + return dag diff --git a/grudge/transform/metadata.py b/grudge/transform/metadata.py new file mode 100644 index 000000000..3b9e4b6a6 --- /dev/null +++ b/grudge/transform/metadata.py @@ -0,0 +1,139 @@ +__copyright__ = "Copyright (C) 2024 Addison Alvey-Blanco" + +__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 pytato as pt +from pytato.transform import ArrayOrNames +from pytato.transform.metadata import ( + AxesTagsEquationCollector as BaseAxesTagsEquationCollector, + AxisIgnoredForPropagationTag, +) + +from arraycontext import ArrayContainer +from arraycontext.container.traversal import rec_map_array_container +from meshmode.discretization import ElementGroupBase +from meshmode.discretization.poly_element import TensorProductElementGroupBase +from meshmode.transform_metadata import ( + DiscretizationDOFAxisTag, + DiscretizationEntityAxisTag, +) +from pytools.tag import Tag, tag_dataclass + + +class OutputIsTensorProductDOFArrayOrdered(Tag): + """ + Signal an eager `arraycontext` to adjust strides to be compatible with + tensor product element DOFs when generating a `loopy` program. + """ + pass + + +@tag_dataclass +class TensorProductDOFAxisTag(DiscretizationEntityAxisTag): + """ + Signify an axis as containing the DOFs of a tensor product discretization. + `iaxis` is later interpreted to determine the relative update speed (i.e. + the stride) of each axis. + """ + iaxis: int + + +class TensorProductOperatorAxisTag(DiscretizationDOFAxisTag, + AxisIgnoredForPropagationTag): + """ + Signify an axis is part of a 1D operator applied to a tensor product + discretization. No tags will be propagated to or along axes containing this + tag. + """ + pass + + +class TensorProductMassOperatorTag(Tag): + """ + Tag an operator as being a reference mass operator. Used to realize an + algebraic simplification of redundant mass-times-mass-inverse operations + when using a tensor product discretization. + """ + pass + + +class TensorProductMassOperatorInverseTag(Tag): + """ + See `TensorProductMassOperatorTag`. + """ + pass + + +class TensorProductStiffnessOperatorTag(Tag): + """ + Similar to `TensorProductMassOperatorTag`. Used to implement an + associativity DAG transformation. + """ + pass + + +def get_dof_axis_tag_type( + element_group: ElementGroupBase, + use_tensor_product_fast_eval: bool = True + ) -> type[Tag]: + if isinstance(element_group, TensorProductElementGroupBase) and \ + use_tensor_product_fast_eval: + return TensorProductOperatorAxisTag + return DiscretizationDOFAxisTag + + +# {{{ solve for discretization metadata for arrays' axes + +class AxesTagsEquationCollector(BaseAxesTagsEquationCollector): + def map_reshape(self, expr: pt.Reshape) -> None: + super().map_reshape(expr) + # + # if (expr.size > 0 + # and (1 not in (expr.array.shape)) # leads to ambiguous newaxis + # and (set(expr.shape) <= (set(expr.array.shape) | {1}))): + # i_in_axis = 0 + # for i_out_axis, dim in enumerate(expr.shape): + # if dim != 1: + # assert dim == expr.array.shape[i_in_axis] + # self.record_equation( + # self.get_var_for_axis(expr.array, + # i_in_axis), + # self.get_var_for_axis(expr, + # i_out_axis) + # ) + # i_in_axis += 1 + # else: + # # print(f"Skipping: {expr.array.shape} -> {expr.shape}") + # # Wacky reshape => bail. + # pass + # pass + + +def unify_discretization_entity_tags(expr: ArrayContainer | ArrayOrNames + ) -> ArrayOrNames: + if not isinstance(expr, pt.Array | pt.DictOfNamedArrays): + return rec_map_array_container(unify_discretization_entity_tags, + expr) + + return pt.unify_axes_tags(expr, + tag_t=DiscretizationEntityAxisTag, + equations_collector_t=AxesTagsEquationCollector) + +# }}} diff --git a/grudge/version.py b/grudge/version.py index 4b139fed7..9ee4d2e75 100644 --- a/grudge/version.py +++ b/grudge/version.py @@ -4,7 +4,7 @@ def _parse_version(version: str) -> tuple[tuple[int, ...], str]: import re - m = re.match("^([0-9.]+)([a-z0-9]*?)$", VERSION_TEXT) + m = re.match(r"^([0-9.]+)([a-z0-9]*?)$", VERSION_TEXT) assert m is not None return tuple(int(nr) for nr in m.group(1).split(".")), m.group(2) diff --git a/test/gh-339-refined-1.msh b/test/gh-339-refined-1.msh new file mode 100644 index 000000000..95c34e825 --- /dev/null +++ b/test/gh-339-refined-1.msh @@ -0,0 +1,2062 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +941 +1 0 0 1 +2 0 0 0 +3 0 1 1 +4 0 1 0 +5 1 0 1 +6 1 0 0 +7 1 1 1 +8 1 1 0 +9 0 0 0.4999999999999999 +10 0 0 0.2499999999999999 +11 0 0 0.75 +12 0 0.4999999999999999 1 +13 0 0.2499999999999999 1 +14 0 0.75 1 +15 0 1 0.4999999999999999 +16 0 1 0.2499999999999999 +17 0 1 0.75 +18 0 0.4999999999999999 0 +19 0 0.2499999999999999 0 +20 0 0.75 0 +21 1 0 0.4999999999999999 +22 1 0 0.2499999999999999 +23 1 0 0.75 +24 1 0.4999999999999999 1 +25 1 0.2499999999999999 1 +26 1 0.75 1 +27 1 1 0.4999999999999999 +28 1 1 0.2499999999999999 +29 1 1 0.75 +30 1 0.4999999999999999 0 +31 1 0.2499999999999999 0 +32 1 0.75 0 +33 0.4999999999999999 0 0 +34 0.2499999999999999 0 0 +35 0.75 0 0 +36 0.4999999999999999 0 1 +37 0.2499999999999999 0 1 +38 0.75 0 1 +39 0.4999999999999999 1 0 +40 0.2499999999999999 1 0 +41 0.75 1 0 +42 0.4999999999999999 1 1 +43 0.2499999999999999 1 1 +44 0.75 1 1 +45 0 0.5 0.5 +46 0 0.25 0.75 +47 0 0.25 0.25 +48 0 0.75 0.75 +49 0 0.75 0.25 +50 0 0.1666666666666667 0.5 +51 0 0.5 0.8333333333333334 +52 0 0.5 0.1666666666666667 +53 0 0.8333333333333334 0.5 +54 0 0.08333333333333336 0.4999999999999999 +55 0 0.2083333333333334 0.375 +56 0 0.125 0.125 +57 0 0.1041666666666667 0.3124999999999999 +58 0 0.125 0.875 +59 0 0.2083333333333334 0.625 +60 0 0.1041666666666667 0.6875 +61 0 0.375 0.625 +62 0 0.375 0.375 +63 0 0.2916666666666667 0.5 +64 0 0.4999999999999999 0.9166666666666667 +65 0 0.375 0.7916666666666667 +66 0 0.3124999999999999 0.8958333333333334 +67 0 0.875 0.875 +68 0 0.625 0.7916666666666667 +69 0 0.6875 0.8958333333333334 +70 0 0.625 0.625 +71 0 0.5 0.7083333333333334 +72 0 0.4999999999999999 0.08333333333333336 +73 0 0.625 0.2083333333333334 +74 0 0.875 0.125 +75 0 0.6875 0.1041666666666667 +76 0 0.375 0.2083333333333334 +77 0 0.3124999999999999 0.1041666666666667 +78 0 0.625 0.375 +79 0 0.5 0.2916666666666667 +80 0 0.9166666666666667 0.4999999999999999 +81 0 0.7916666666666667 0.625 +82 0 0.8958333333333334 0.6875 +83 0 0.7916666666666667 0.375 +84 0 0.8958333333333334 0.3124999999999999 +85 0 0.7083333333333334 0.5 +86 1 0.5 0.5 +87 1 0.25 0.25 +88 1 0.25 0.75 +89 1 0.75 0.75 +90 1 0.75 0.25 +91 1 0.1666666666666667 0.5 +92 1 0.5 0.8333333333333334 +93 1 0.5 0.1666666666666667 +94 1 0.8333333333333334 0.5 +95 1 0.125 0.125 +96 1 0.2083333333333334 0.375 +97 1 0.08333333333333336 0.4999999999999999 +98 1 0.1041666666666667 0.3124999999999999 +99 1 0.375 0.375 +100 1 0.375 0.625 +101 1 0.2083333333333334 0.625 +102 1 0.2916666666666667 0.5 +103 1 0.125 0.875 +104 1 0.1041666666666667 0.6875 +105 1 0.375 0.7916666666666667 +106 1 0.4999999999999999 0.9166666666666667 +107 1 0.3124999999999999 0.8958333333333334 +108 1 0.625 0.625 +109 1 0.625 0.7916666666666667 +110 1 0.5 0.7083333333333334 +111 1 0.875 0.875 +112 1 0.6875 0.8958333333333334 +113 1 0.875 0.125 +114 1 0.625 0.2083333333333334 +115 1 0.4999999999999999 0.08333333333333336 +116 1 0.6875 0.1041666666666667 +117 1 0.625 0.375 +118 1 0.375 0.2083333333333334 +119 1 0.5 0.2916666666666667 +120 1 0.3124999999999999 0.1041666666666667 +121 1 0.7916666666666667 0.625 +122 1 0.9166666666666667 0.4999999999999999 +123 1 0.8958333333333334 0.6875 +124 1 0.7916666666666667 0.375 +125 1 0.7083333333333334 0.5 +126 1 0.8958333333333334 0.3124999999999999 +127 0.5 0 0.5 +128 0.25 0 0.25 +129 0.25 0 0.75 +130 0.75 0 0.75 +131 0.75 0 0.25 +132 0.1666666666666667 0 0.5 +133 0.5 0 0.8333333333333334 +134 0.5 0 0.1666666666666667 +135 0.8333333333333334 0 0.5 +136 0.08333333333333336 0 0.4999999999999999 +137 0.2083333333333334 0 0.625 +138 0.125 0 0.875 +139 0.1041666666666667 0 0.6875 +140 0.125 0 0.125 +141 0.2083333333333334 0 0.375 +142 0.1041666666666667 0 0.3124999999999999 +143 0.375 0 0.375 +144 0.375 0 0.625 +145 0.2916666666666667 0 0.5 +146 0.4999999999999999 0 0.9166666666666667 +147 0.625 0 0.7916666666666667 +148 0.875 0 0.875 +149 0.6875 0 0.8958333333333334 +150 0.375 0 0.7916666666666667 +151 0.3124999999999999 0 0.8958333333333334 +152 0.625 0 0.625 +153 0.5 0 0.7083333333333334 +154 0.4999999999999999 0 0.08333333333333336 +155 0.375 0 0.2083333333333334 +156 0.3124999999999999 0 0.1041666666666667 +157 0.875 0 0.125 +158 0.625 0 0.2083333333333334 +159 0.6875 0 0.1041666666666667 +160 0.625 0 0.375 +161 0.5 0 0.2916666666666667 +162 0.9166666666666667 0 0.4999999999999999 +163 0.7916666666666667 0 0.375 +164 0.8958333333333334 0 0.3124999999999999 +165 0.7916666666666667 0 0.625 +166 0.8958333333333334 0 0.6875 +167 0.7083333333333334 0 0.5 +168 0.5 1 0.5 +169 0.25 1 0.75 +170 0.25 1 0.25 +171 0.75 1 0.75 +172 0.75 1 0.25 +173 0.1666666666666667 1 0.5 +174 0.5 1 0.8333333333333334 +175 0.5 1 0.1666666666666667 +176 0.8333333333333334 1 0.5 +177 0.125 1 0.875 +178 0.2083333333333334 1 0.625 +179 0.08333333333333336 1 0.4999999999999999 +180 0.1041666666666667 1 0.6875 +181 0.375 1 0.625 +182 0.375 1 0.375 +183 0.2083333333333334 1 0.375 +184 0.2916666666666667 1 0.5 +185 0.125 1 0.125 +186 0.1041666666666667 1 0.3124999999999999 +187 0.875 1 0.875 +188 0.625 1 0.7916666666666667 +189 0.4999999999999999 1 0.9166666666666667 +190 0.6875 1 0.8958333333333334 +191 0.625 1 0.625 +192 0.375 1 0.7916666666666667 +193 0.5 1 0.7083333333333334 +194 0.3124999999999999 1 0.8958333333333334 +195 0.375 1 0.2083333333333334 +196 0.4999999999999999 1 0.08333333333333336 +197 0.3124999999999999 1 0.1041666666666667 +198 0.625 1 0.375 +199 0.625 1 0.2083333333333334 +200 0.5 1 0.2916666666666667 +201 0.875 1 0.125 +202 0.6875 1 0.1041666666666667 +203 0.7916666666666667 1 0.375 +204 0.9166666666666667 1 0.4999999999999999 +205 0.8958333333333334 1 0.3124999999999999 +206 0.7916666666666667 1 0.625 +207 0.7083333333333334 1 0.5 +208 0.8958333333333334 1 0.6875 +209 0.5 0.5 0 +210 0.25 0.75 0 +211 0.25 0.25 0 +212 0.75 0.25 0 +213 0.75 0.75 0 +214 0.1666666666666667 0.5 0 +215 0.5 0.1666666666666667 0 +216 0.5 0.8333333333333334 0 +217 0.8333333333333334 0.5 0 +218 0.08333333333333336 0.4999999999999999 0 +219 0.2083333333333334 0.375 0 +220 0.125 0.125 0 +221 0.1041666666666667 0.3124999999999999 0 +222 0.125 0.875 0 +223 0.2083333333333334 0.625 0 +224 0.1041666666666667 0.6875 0 +225 0.375 0.625 0 +226 0.375 0.375 0 +227 0.2916666666666667 0.5 0 +228 0.4999999999999999 0.08333333333333336 0 +229 0.625 0.2083333333333334 0 +230 0.875 0.125 0 +231 0.6875 0.1041666666666667 0 +232 0.375 0.2083333333333334 0 +233 0.3124999999999999 0.1041666666666667 0 +234 0.625 0.375 0 +235 0.5 0.2916666666666667 0 +236 0.4999999999999999 0.9166666666666667 0 +237 0.375 0.7916666666666667 0 +238 0.3124999999999999 0.8958333333333334 0 +239 0.875 0.875 0 +240 0.625 0.7916666666666667 0 +241 0.6875 0.8958333333333334 0 +242 0.625 0.625 0 +243 0.5 0.7083333333333334 0 +244 0.9166666666666667 0.4999999999999999 0 +245 0.7916666666666667 0.625 0 +246 0.8958333333333334 0.6875 0 +247 0.7916666666666667 0.375 0 +248 0.8958333333333334 0.3124999999999999 0 +249 0.7083333333333334 0.5 0 +250 0.5 0.5 1 +251 0.25 0.25 1 +252 0.25 0.75 1 +253 0.75 0.25 1 +254 0.75 0.75 1 +255 0.1666666666666667 0.5 1 +256 0.5 0.1666666666666667 1 +257 0.5 0.8333333333333334 1 +258 0.8333333333333334 0.5 1 +259 0.125 0.125 1 +260 0.2083333333333334 0.375 1 +261 0.08333333333333336 0.4999999999999999 1 +262 0.1041666666666667 0.3124999999999999 1 +263 0.375 0.375 1 +264 0.375 0.625 1 +265 0.2083333333333334 0.625 1 +266 0.2916666666666667 0.5 1 +267 0.125 0.875 1 +268 0.1041666666666667 0.6875 1 +269 0.875 0.125 1 +270 0.625 0.2083333333333334 1 +271 0.4999999999999999 0.08333333333333336 1 +272 0.6875 0.1041666666666667 1 +273 0.625 0.375 1 +274 0.375 0.2083333333333334 1 +275 0.5 0.2916666666666667 1 +276 0.3124999999999999 0.1041666666666667 1 +277 0.375 0.7916666666666667 1 +278 0.4999999999999999 0.9166666666666667 1 +279 0.3124999999999999 0.8958333333333334 1 +280 0.625 0.625 1 +281 0.625 0.7916666666666667 1 +282 0.5 0.7083333333333334 1 +283 0.875 0.875 1 +284 0.6875 0.8958333333333334 1 +285 0.7916666666666667 0.625 1 +286 0.9166666666666667 0.4999999999999999 1 +287 0.8958333333333334 0.6875 1 +288 0.7916666666666667 0.375 1 +289 0.7083333333333334 0.5 1 +290 0.8958333333333334 0.3124999999999999 1 +291 0.75 0.5 0.25 +292 0.5 0.5 0.5 +293 0.25 0.5 0.25 +294 0.5 0.25 0.25 +295 0.25 0.25 0.5 +296 0.75 0.25 0.5 +297 0.25 0.5 0.75 +298 0.25 0.75 0.5 +299 0.5 0.75 0.75 +300 0.75 0.5 0.75 +301 0.75 0.75 0.5 +302 0.5 0.75 0.25 +303 0.5 0.25 0.75 +304 0.5 0.5 0.3333333333333333 +305 0.6666666666666666 0.3333333333333333 0.3333333333333333 +306 0.3333333333333333 0.3333333333333333 0.3333333333333333 +307 0.5 0.3333333333333333 0.5 +308 0.5 0.375 0.375 +309 0.3333333333333333 0.6666666666666666 0.6666666666666666 +310 0.5 0.5 0.6666666666666666 +311 0.6666666666666666 0.6666666666666666 0.6666666666666666 +312 0.5 0.6666666666666666 0.5 +313 0.5 0.625 0.625 +314 0.6666666666666666 0.6666666666666666 0.3333333333333333 +315 0.3333333333333333 0.6666666666666666 0.3333333333333333 +316 0.5 0.625 0.375 +317 0.3333333333333333 0.3333333333333333 0.6666666666666666 +318 0.6666666666666666 0.3333333333333333 0.6666666666666666 +319 0.5 0.375 0.625 +320 0.8333333333333334 0.8333333333333334 0.6666666666666666 +321 0.8333333333333334 0.8333333333333334 0.3333333333333333 +322 0.875 0.875 0.5 +323 0.1666666666666667 0.8333333333333334 0.3333333333333333 +324 0.1666666666666667 0.8333333333333334 0.6666666666666666 +325 0.125 0.875 0.5 +326 0.6666666666666666 0.1666666666666667 0.8333333333333334 +327 0.3333333333333333 0.1666666666666667 0.8333333333333334 +328 0.5 0.125 0.875 +329 0.1666666666666667 0.6666666666666666 0.8333333333333334 +330 0.1666666666666667 0.3333333333333333 0.8333333333333334 +331 0.125 0.5 0.875 +332 0.3333333333333333 0.8333333333333334 0.8333333333333334 +333 0.6666666666666666 0.8333333333333334 0.8333333333333334 +334 0.5 0.875 0.875 +335 0.8333333333333334 0.6666666666666666 0.1666666666666667 +336 0.8333333333333334 0.3333333333333333 0.1666666666666667 +337 0.875 0.5 0.125 +338 0.8333333333333334 0.6666666666666666 0.8333333333333334 +339 0.8333333333333334 0.3333333333333333 0.8333333333333334 +340 0.875 0.5 0.875 +341 0.1666666666666667 0.3333333333333333 0.1666666666666667 +342 0.1666666666666667 0.6666666666666666 0.1666666666666667 +343 0.125 0.5 0.125 +344 0.1666666666666667 0.1666666666666667 0.6666666666666666 +345 0.1666666666666667 0.1666666666666667 0.3333333333333333 +346 0.125 0.125 0.5 +347 0.3333333333333333 0.8333333333333334 0.1666666666666667 +348 0.6666666666666666 0.8333333333333334 0.1666666666666667 +349 0.5 0.875 0.125 +350 0.6666666666666666 0.1666666666666667 0.1666666666666667 +351 0.3333333333333333 0.1666666666666667 0.1666666666666667 +352 0.5 0.125 0.125 +353 0.8333333333333334 0.1666666666666667 0.3333333333333333 +354 0.8333333333333334 0.1666666666666667 0.6666666666666666 +355 0.875 0.125 0.5 +356 0.25 0.25 0.75 +357 0.25 0.75 0.75 +358 0.75 0.75 0.75 +359 0.25 0.75 0.25 +360 0.75 0.75 0.25 +361 0.75 0.25 0.25 +362 0.75 0.25 0.75 +363 0.25 0.25 0.25 +364 0.625 0.5 0.125 +365 0.375 0.5 0.125 +366 0.5 0.375 0.125 +367 0.625 0.5 0.2916666666666666 +368 0.7083333333333333 0.4166666666666666 0.2916666666666666 +369 0.375 0.5 0.2916666666666666 +370 0.5 0.4375 0.3541666666666666 +371 0.2916666666666666 0.4166666666666666 0.2916666666666666 +372 0.5833333333333333 0.2916666666666666 0.2916666666666666 +373 0.4166666666666666 0.2916666666666666 0.2916666666666666 +374 0.5833333333333333 0.3541666666666666 0.3541666666666666 +375 0.4166666666666666 0.3541666666666666 0.3541666666666666 +376 0.5 0.5 0.2083333333333333 +377 0.6041666666666666 0.3958333333333333 0.2083333333333333 +378 0.3958333333333333 0.3958333333333333 0.2083333333333333 +379 0.6041666666666666 0.4270833333333333 0.3229166666666666 +380 0.3958333333333333 0.4270833333333333 0.3229166666666666 +381 0.4999999999999999 0.3229166666666666 0.3229166666666666 +382 0.5 0.4114583333333335 0.265625 +383 0.875 0.5 0.375 +384 0.75 0.5 0.5 +385 0.875 0.375 0.5 +386 0.5 0.5 0.4166666666666666 +387 0.5 0.4166666666666666 0.5 +388 0.7083333333333333 0.2916666666666666 0.4166666666666666 +389 0.625 0.2916666666666666 0.5 +390 0.5 0.3541666666666666 0.4375 +391 0.6875 0.5 0.3958333333333333 +392 0.7916666666666666 0.3958333333333333 0.3958333333333333 +393 0.6875 0.3958333333333333 0.5 +394 0.5 0.4270833333333333 0.4270833333333333 +395 0.6041666666666666 0.3229166666666666 0.4270833333333333 +396 0.6458333333333334 0.4114583333333335 0.4114583333333334 +397 0.125 0.5 0.375 +398 0.25 0.5 0.5 +399 0.125 0.375 0.5 +400 0.2916666666666666 0.2916666666666666 0.4166666666666666 +401 0.375 0.2916666666666666 0.5 +402 0.3125 0.5 0.3958333333333333 +403 0.2083333333333333 0.3958333333333333 0.3958333333333333 +404 0.3125 0.3958333333333333 0.5 +405 0.3958333333333333 0.3229166666666666 0.4270833333333333 +406 0.3541666666666667 0.4114583333333335 0.4114583333333333 +407 0.625 0.125 0.5 +408 0.5 0.125 0.375 +409 0.375 0.125 0.5 +410 0.6041666666666666 0.2083333333333333 0.3958333333333333 +411 0.5 0.2083333333333333 0.5 +412 0.3958333333333333 0.2083333333333333 0.3958333333333333 +413 0.5000000000000001 0.265625 0.4114583333333334 +414 0.375 0.5 0.875 +415 0.5 0.625 0.875 +416 0.625 0.5 0.875 +417 0.2916666666666666 0.5833333333333333 0.7083333333333333 +418 0.375 0.5 0.7083333333333333 +419 0.4166666666666666 0.7083333333333333 0.7083333333333333 +420 0.4166666666666666 0.6458333333333333 0.6458333333333333 +421 0.5833333333333333 0.7083333333333333 0.7083333333333333 +422 0.625 0.5 0.7083333333333333 +423 0.7083333333333333 0.5833333333333333 0.7083333333333333 +424 0.5 0.5625 0.6458333333333333 +425 0.5833333333333333 0.6458333333333333 0.6458333333333333 +426 0.3958333333333333 0.6041666666666666 0.7916666666666666 +427 0.5 0.5 0.7916666666666666 +428 0.6041666666666666 0.6041666666666666 0.7916666666666666 +429 0.3958333333333333 0.5729166666666666 0.6770833333333333 +430 0.4999999999999999 0.6770833333333333 0.6770833333333333 +431 0.6041666666666666 0.5729166666666666 0.6770833333333333 +432 0.5000000000000001 0.5885416666666665 0.7343750000000002 +433 0.125 0.5 0.625 +434 0.125 0.625 0.5 +435 0.2916666666666666 0.7083333333333333 0.5833333333333333 +436 0.375 0.7083333333333333 0.5 +437 0.5 0.5 0.5833333333333333 +438 0.5 0.5833333333333333 0.5 +439 0.5 0.6458333333333333 0.5625 +440 0.2083333333333333 0.6041666666666666 0.6041666666666666 +441 0.3125 0.5 0.6041666666666666 +442 0.3125 0.6041666666666666 0.5 +443 0.3958333333333333 0.6770833333333333 0.5729166666666666 +444 0.5 0.5729166666666666 0.5729166666666666 +445 0.3541666666666665 0.5885416666666665 0.5885416666666665 +446 0.5 0.875 0.625 +447 0.375 0.875 0.5 +448 0.625 0.875 0.5 +449 0.7083333333333333 0.7083333333333333 0.5833333333333333 +450 0.625 0.7083333333333333 0.5 +451 0.3958333333333333 0.7916666666666666 0.6041666666666666 +452 0.6041666666666666 0.7916666666666666 0.6041666666666666 +453 0.5 0.7916666666666666 0.5 +454 0.6041666666666666 0.6770833333333333 0.5729166666666666 +455 0.5 0.7343750000000002 0.5885416666666665 +456 0.875 0.5 0.625 +457 0.875 0.625 0.5 +458 0.6875 0.5 0.6041666666666666 +459 0.6875 0.6041666666666666 0.5 +460 0.7916666666666666 0.6041666666666666 0.6041666666666666 +461 0.6458333333333335 0.5885416666666667 0.5885416666666665 +462 0.5 0.6458333333333333 0.4375 +463 0.7083333333333333 0.7083333333333333 0.4166666666666666 +464 0.7083333333333333 0.5833333333333333 0.2916666666666666 +465 0.5 0.5625 0.3541666666666666 +466 0.5833333333333333 0.6458333333333333 0.3541666666666666 +467 0.7916666666666666 0.6041666666666666 0.3958333333333333 +468 0.5 0.5729166666666666 0.4270833333333333 +469 0.6041666666666666 0.6770833333333333 0.4270833333333333 +470 0.6041666666666666 0.5729166666666666 0.3229166666666666 +471 0.6458333333333335 0.5885416666666665 0.4114583333333335 +472 0.2916666666666666 0.7083333333333333 0.4166666666666666 +473 0.2916666666666666 0.5833333333333333 0.2916666666666666 +474 0.4166666666666666 0.6458333333333333 0.3541666666666666 +475 0.2083333333333333 0.6041666666666666 0.3958333333333333 +476 0.3958333333333333 0.6770833333333333 0.4270833333333333 +477 0.3958333333333333 0.5729166666666666 0.3229166666666666 +478 0.3541666666666667 0.5885416666666665 0.4114583333333335 +479 0.5 0.875 0.375 +480 0.5833333333333333 0.7083333333333333 0.2916666666666666 +481 0.4166666666666666 0.7083333333333333 0.2916666666666666 +482 0.6041666666666666 0.7916666666666666 0.3958333333333333 +483 0.3958333333333333 0.7916666666666666 0.3958333333333333 +484 0.4999999999999999 0.6770833333333333 0.3229166666666666 +485 0.5000000000000001 0.7343750000000002 0.4114583333333335 +486 0.5 0.625 0.125 +487 0.3958333333333333 0.6041666666666666 0.2083333333333333 +488 0.6041666666666666 0.6041666666666666 0.2083333333333333 +489 0.4999999999999999 0.5885416666666667 0.265625 +490 0.5 0.375 0.875 +491 0.2916666666666666 0.4166666666666666 0.7083333333333333 +492 0.5 0.4375 0.6458333333333333 +493 0.7083333333333333 0.4166666666666666 0.7083333333333333 +494 0.4166666666666666 0.2916666666666666 0.7083333333333333 +495 0.5833333333333333 0.2916666666666666 0.7083333333333333 +496 0.4166666666666666 0.3541666666666666 0.6458333333333333 +497 0.5833333333333333 0.3541666666666666 0.6458333333333333 +498 0.3958333333333333 0.3958333333333333 0.7916666666666666 +499 0.6041666666666666 0.3958333333333333 0.7916666666666666 +500 0.3958333333333333 0.4270833333333333 0.6770833333333333 +501 0.6041666666666666 0.4270833333333333 0.6770833333333333 +502 0.4999999999999999 0.3229166666666666 0.6770833333333333 +503 0.4999999999999999 0.4114583333333335 0.7343750000000002 +504 0.2916666666666666 0.2916666666666666 0.5833333333333333 +505 0.5 0.3541666666666666 0.5625 +506 0.2083333333333333 0.3958333333333333 0.6041666666666666 +507 0.5 0.4270833333333333 0.5729166666666666 +508 0.3958333333333333 0.3229166666666666 0.5729166666666666 +509 0.3541666666666667 0.4114583333333335 0.5885416666666665 +510 0.7083333333333333 0.2916666666666666 0.5833333333333333 +511 0.7916666666666666 0.3958333333333333 0.6041666666666666 +512 0.6041666666666666 0.3229166666666666 0.5729166666666666 +513 0.6458333333333334 0.4114583333333335 0.5885416666666665 +514 0.5 0.125 0.625 +515 0.3958333333333333 0.2083333333333333 0.6041666666666666 +516 0.6041666666666666 0.2083333333333333 0.6041666666666666 +517 0.5000000000000001 0.265625 0.5885416666666665 +518 0.7916666666666667 0.9166666666666667 0.7083333333333333 +519 0.7916666666666667 0.7916666666666667 0.5833333333333333 +520 0.8541666666666667 0.8541666666666667 0.5833333333333333 +521 0.7916666666666667 0.7916666666666667 0.4166666666666666 +522 0.7916666666666667 0.9166666666666667 0.2916666666666666 +523 0.8541666666666667 0.9375 0.5 +524 0.8541666666666667 0.8541666666666667 0.4166666666666666 +525 0.7083333333333334 0.8958333333333334 0.6041666666666666 +526 0.7083333333333334 0.8958333333333334 0.3958333333333333 +527 0.8229166666666667 0.9270833333333334 0.6041666666666666 +528 0.8229166666666667 0.8229166666666667 0.4999999999999999 +529 0.8229166666666667 0.9270833333333334 0.3958333333333333 +530 0.7656249999999998 0.9114583333333333 0.5000000000000001 +531 0.9166666666666667 0.7916666666666667 0.7083333333333333 +532 0.9375 0.8541666666666667 0.5 +533 0.8958333333333334 0.8958333333333334 0.7916666666666666 +534 0.9270833333333334 0.8229166666666667 0.6041666666666666 +535 0.9270833333333334 0.9270833333333334 0.5 +536 0.9114583333333333 0.9114583333333333 0.6458333333333331 +537 0.9166666666666667 0.7916666666666667 0.2916666666666666 +538 0.8958333333333334 0.7083333333333334 0.6041666666666666 +539 0.8958333333333334 0.7083333333333334 0.3958333333333333 +540 0.9270833333333334 0.8229166666666667 0.3958333333333333 +541 0.9114583333333334 0.7656249999999998 0.5000000000000001 +542 0.8958333333333334 0.8958333333333334 0.2083333333333333 +543 0.9114583333333333 0.911458333333333 0.3541666666666667 +544 0.2083333333333334 0.7916666666666667 0.4166666666666666 +545 0.2083333333333334 0.7916666666666667 0.5833333333333333 +546 0.08333333333333336 0.7916666666666667 0.2916666666666666 +547 0.1458333333333334 0.8541666666666667 0.4166666666666666 +548 0.08333333333333336 0.7916666666666667 0.7083333333333333 +549 0.1458333333333334 0.8541666666666667 0.5833333333333333 +550 0.0625 0.8541666666666667 0.5 +551 0.1041666666666667 0.7083333333333334 0.3958333333333333 +552 0.1041666666666667 0.7083333333333334 0.6041666666666666 +553 0.1770833333333334 0.8229166666666667 0.4999999999999999 +554 0.07291666666666669 0.8229166666666667 0.3958333333333333 +555 0.07291666666666669 0.8229166666666667 0.6041666666666666 +556 0.08854166666666669 0.7656249999999998 0.4999999999999999 +557 0.2083333333333334 0.9166666666666667 0.2916666666666666 +558 0.2083333333333334 0.9166666666666667 0.7083333333333333 +559 0.1458333333333334 0.9375 0.5 +560 0.2916666666666667 0.8958333333333334 0.3958333333333333 +561 0.2916666666666667 0.8958333333333334 0.6041666666666666 +562 0.1770833333333334 0.9270833333333334 0.3958333333333333 +563 0.1770833333333334 0.9270833333333334 0.6041666666666666 +564 0.2343749999999999 0.9114583333333333 0.5000000000000001 +565 0.1041666666666667 0.8958333333333334 0.2083333333333333 +566 0.07291666666666669 0.9270833333333334 0.4999999999999999 +567 0.08854166666666669 0.9114583333333334 0.3541666666666666 +568 0.1041666666666667 0.8958333333333334 0.7916666666666666 +569 0.08854166666666669 0.9114583333333333 0.6458333333333331 +570 0.5 0.1458333333333334 0.9375 +571 0.7083333333333333 0.2083333333333334 0.9166666666666667 +572 0.7083333333333333 0.08333333333333336 0.7916666666666667 +573 0.5 0.0625 0.8541666666666667 +574 0.5833333333333333 0.1458333333333334 0.8541666666666667 +575 0.7916666666666666 0.1041666666666667 0.8958333333333334 +576 0.5 0.07291666666666669 0.9270833333333334 +577 0.6041666666666666 0.1770833333333334 0.9270833333333334 +578 0.6041666666666666 0.07291666666666669 0.8229166666666667 +579 0.6458333333333335 0.08854166666666669 0.911458333333333 +580 0.2916666666666666 0.2083333333333334 0.9166666666666667 +581 0.2916666666666666 0.08333333333333336 0.7916666666666667 +582 0.4166666666666666 0.1458333333333334 0.8541666666666667 +583 0.2083333333333333 0.1041666666666667 0.8958333333333334 +584 0.3958333333333333 0.1770833333333334 0.9270833333333334 +585 0.3958333333333333 0.07291666666666669 0.8229166666666667 +586 0.3541666666666667 0.08854166666666669 0.911458333333333 +587 0.5833333333333333 0.2083333333333334 0.7916666666666667 +588 0.4166666666666666 0.2083333333333334 0.7916666666666667 +589 0.6041666666666666 0.2916666666666667 0.8958333333333334 +590 0.3958333333333333 0.2916666666666667 0.8958333333333334 +591 0.4999999999999999 0.1770833333333334 0.8229166666666667 +592 0.5000000000000001 0.234375 0.911458333333333 +593 0.3958333333333333 0.1041666666666667 0.7083333333333334 +594 0.6041666666666666 0.1041666666666667 0.7083333333333334 +595 0.4999999999999999 0.08854166666666669 0.7656249999999998 +596 0.2083333333333334 0.7083333333333333 0.9166666666666667 +597 0.2083333333333334 0.5833333333333333 0.7916666666666667 +598 0.1458333333333334 0.5833333333333333 0.8541666666666667 +599 0.2083333333333334 0.4166666666666666 0.7916666666666667 +600 0.2083333333333334 0.2916666666666666 0.9166666666666667 +601 0.1458333333333334 0.5 0.9375 +602 0.1458333333333334 0.4166666666666666 0.8541666666666667 +603 0.2916666666666667 0.6041666666666666 0.8958333333333334 +604 0.2916666666666667 0.3958333333333333 0.8958333333333334 +605 0.1770833333333334 0.6041666666666666 0.9270833333333334 +606 0.1770833333333334 0.4999999999999999 0.8229166666666667 +607 0.1770833333333334 0.3958333333333333 0.9270833333333334 +608 0.234375 0.5000000000000001 0.9114583333333333 +609 0.08333333333333336 0.7083333333333333 0.7916666666666667 +610 0.0625 0.5 0.8541666666666667 +611 0.1041666666666667 0.7916666666666666 0.8958333333333334 +612 0.07291666666666669 0.6041666666666666 0.8229166666666667 +613 0.07291666666666669 0.5 0.9270833333333334 +614 0.08854166666666669 0.6458333333333331 0.9114583333333333 +615 0.08333333333333336 0.2916666666666666 0.7916666666666667 +616 0.1041666666666667 0.6041666666666666 0.7083333333333334 +617 0.1041666666666667 0.3958333333333333 0.7083333333333334 +618 0.07291666666666669 0.3958333333333333 0.8229166666666667 +619 0.08854166666666669 0.5000000000000001 0.7656249999999998 +620 0.1041666666666667 0.2083333333333333 0.8958333333333334 +621 0.08854166666666669 0.3541666666666667 0.911458333333333 +622 0.2916666666666666 0.9166666666666667 0.7916666666666667 +623 0.4166666666666666 0.7916666666666667 0.7916666666666667 +624 0.4166666666666666 0.8541666666666667 0.8541666666666667 +625 0.5833333333333333 0.7916666666666667 0.7916666666666667 +626 0.7083333333333333 0.9166666666666667 0.7916666666666667 +627 0.5 0.9375 0.8541666666666667 +628 0.5833333333333333 0.8541666666666667 0.8541666666666667 +629 0.3958333333333333 0.8958333333333334 0.7083333333333334 +630 0.6041666666666666 0.8958333333333334 0.7083333333333334 +631 0.3958333333333333 0.9270833333333334 0.8229166666666667 +632 0.4999999999999999 0.8229166666666667 0.8229166666666667 +633 0.6041666666666666 0.9270833333333334 0.8229166666666667 +634 0.5000000000000001 0.9114583333333333 0.7656249999999998 +635 0.2916666666666666 0.7916666666666667 0.9166666666666667 +636 0.5 0.8541666666666667 0.9375 +637 0.2083333333333333 0.8958333333333334 0.8958333333333334 +638 0.3958333333333333 0.8229166666666667 0.9270833333333334 +639 0.5 0.9270833333333334 0.9270833333333334 +640 0.3541666666666666 0.9114583333333333 0.9114583333333333 +641 0.7083333333333333 0.7916666666666667 0.9166666666666667 +642 0.3958333333333333 0.7083333333333334 0.8958333333333334 +643 0.6041666666666666 0.7083333333333334 0.8958333333333334 +644 0.6041666666666666 0.8229166666666667 0.9270833333333334 +645 0.5 0.7656249999999998 0.9114583333333334 +646 0.7916666666666666 0.8958333333333334 0.8958333333333334 +647 0.6458333333333335 0.911458333333333 0.9114583333333333 +648 0.7916666666666667 0.7083333333333333 0.08333333333333336 +649 0.7916666666666667 0.5833333333333333 0.2083333333333334 +650 0.8541666666666667 0.5833333333333333 0.1458333333333334 +651 0.7916666666666667 0.4166666666666666 0.2083333333333334 +652 0.7916666666666667 0.2916666666666666 0.08333333333333336 +653 0.8541666666666667 0.5 0.0625 +654 0.8541666666666667 0.4166666666666666 0.1458333333333334 +655 0.7083333333333334 0.6041666666666666 0.1041666666666667 +656 0.7083333333333334 0.3958333333333333 0.1041666666666667 +657 0.8229166666666667 0.6041666666666666 0.07291666666666669 +658 0.8229166666666667 0.4999999999999999 0.1770833333333334 +659 0.8229166666666667 0.3958333333333333 0.07291666666666669 +660 0.7656249999999998 0.5000000000000001 0.08854166666666669 +661 0.9166666666666667 0.7083333333333333 0.2083333333333334 +662 0.9375 0.5 0.1458333333333334 +663 0.8958333333333334 0.7916666666666666 0.1041666666666667 +664 0.9270833333333334 0.6041666666666666 0.1770833333333334 +665 0.9270833333333334 0.5 0.07291666666666669 +666 0.9114583333333333 0.6458333333333331 0.08854166666666669 +667 0.9166666666666667 0.2916666666666666 0.2083333333333334 +668 0.8958333333333334 0.6041666666666666 0.2916666666666667 +669 0.8958333333333334 0.3958333333333333 0.2916666666666667 +670 0.9270833333333334 0.3958333333333333 0.1770833333333334 +671 0.9114583333333334 0.5000000000000001 0.234375 +672 0.8958333333333334 0.2083333333333333 0.1041666666666667 +673 0.9114583333333333 0.3541666666666667 0.08854166666666669 +674 0.8541666666666667 0.5 0.9375 +675 0.7916666666666667 0.7083333333333333 0.9166666666666667 +676 0.9166666666666667 0.7083333333333333 0.7916666666666667 +677 0.9375 0.5 0.8541666666666667 +678 0.8541666666666667 0.5833333333333333 0.8541666666666667 +679 0.8958333333333334 0.7916666666666666 0.8958333333333334 +680 0.9270833333333334 0.5 0.9270833333333334 +681 0.8229166666666667 0.6041666666666666 0.9270833333333334 +682 0.9270833333333334 0.6041666666666666 0.8229166666666667 +683 0.9114583333333333 0.6458333333333335 0.911458333333333 +684 0.7916666666666667 0.2916666666666666 0.9166666666666667 +685 0.9166666666666667 0.2916666666666666 0.7916666666666667 +686 0.8541666666666667 0.4166666666666666 0.8541666666666667 +687 0.8958333333333334 0.2083333333333333 0.8958333333333334 +688 0.8229166666666667 0.3958333333333333 0.9270833333333334 +689 0.9270833333333334 0.3958333333333333 0.8229166666666667 +690 0.9114583333333333 0.3541666666666667 0.911458333333333 +691 0.7916666666666667 0.5833333333333333 0.7916666666666667 +692 0.7916666666666667 0.4166666666666666 0.7916666666666667 +693 0.7083333333333334 0.6041666666666666 0.8958333333333334 +694 0.7083333333333334 0.3958333333333333 0.8958333333333334 +695 0.8229166666666667 0.4999999999999999 0.8229166666666667 +696 0.7656249999999998 0.5000000000000001 0.911458333333333 +697 0.8958333333333334 0.3958333333333333 0.7083333333333334 +698 0.8958333333333334 0.6041666666666666 0.7083333333333334 +699 0.911458333333333 0.4999999999999999 0.7656249999999998 +700 0.08333333333333336 0.2916666666666666 0.2083333333333334 +701 0.0625 0.5 0.1458333333333334 +702 0.08333333333333336 0.7083333333333333 0.2083333333333334 +703 0.2083333333333334 0.4166666666666666 0.2083333333333334 +704 0.2083333333333334 0.5833333333333333 0.2083333333333334 +705 0.1458333333333334 0.4166666666666666 0.1458333333333334 +706 0.1458333333333334 0.5833333333333333 0.1458333333333334 +707 0.1041666666666667 0.3958333333333333 0.2916666666666667 +708 0.1041666666666667 0.6041666666666666 0.2916666666666667 +709 0.07291666666666669 0.3958333333333333 0.1770833333333334 +710 0.07291666666666669 0.6041666666666666 0.1770833333333334 +711 0.1770833333333334 0.4999999999999999 0.1770833333333334 +712 0.08854166666666669 0.4999999999999999 0.234375 +713 0.2083333333333334 0.2916666666666666 0.08333333333333336 +714 0.1458333333333334 0.5 0.0625 +715 0.1041666666666667 0.2083333333333333 0.1041666666666667 +716 0.07291666666666669 0.5 0.07291666666666669 +717 0.1770833333333334 0.3958333333333333 0.07291666666666669 +718 0.08854166666666669 0.3541666666666666 0.08854166666666669 +719 0.2083333333333334 0.7083333333333333 0.08333333333333336 +720 0.1041666666666667 0.7916666666666666 0.1041666666666667 +721 0.1770833333333334 0.6041666666666666 0.07291666666666669 +722 0.08854166666666669 0.6458333333333334 0.08854166666666669 +723 0.2916666666666667 0.3958333333333333 0.1041666666666667 +724 0.2916666666666667 0.6041666666666666 0.1041666666666667 +725 0.234375 0.5000000000000001 0.08854166666666669 +726 0.0625 0.1458333333333334 0.5 +727 0.08333333333333336 0.2083333333333334 0.7083333333333333 +728 0.2083333333333334 0.08333333333333336 0.7083333333333333 +729 0.1458333333333334 0.0625 0.5 +730 0.1458333333333334 0.1458333333333334 0.5833333333333333 +731 0.1041666666666667 0.1041666666666667 0.7916666666666666 +732 0.07291666666666669 0.07291666666666669 0.5 +733 0.07291666666666669 0.1770833333333334 0.6041666666666666 +734 0.1770833333333334 0.07291666666666669 0.6041666666666666 +735 0.08854166666666669 0.08854166666666669 0.6458333333333335 +736 0.08333333333333336 0.2083333333333334 0.2916666666666666 +737 0.2083333333333334 0.08333333333333336 0.2916666666666666 +738 0.1458333333333334 0.1458333333333334 0.4166666666666666 +739 0.1041666666666667 0.1041666666666667 0.2083333333333333 +740 0.07291666666666669 0.1770833333333334 0.3958333333333333 +741 0.1770833333333334 0.07291666666666669 0.3958333333333333 +742 0.08854166666666669 0.08854166666666669 0.3541666666666667 +743 0.2083333333333334 0.2083333333333334 0.5833333333333333 +744 0.2083333333333334 0.2083333333333334 0.4166666666666666 +745 0.1041666666666667 0.2916666666666667 0.6041666666666666 +746 0.1041666666666667 0.2916666666666667 0.3958333333333333 +747 0.1770833333333334 0.1770833333333334 0.4999999999999999 +748 0.08854166666666669 0.234375 0.5000000000000001 +749 0.2916666666666667 0.1041666666666667 0.3958333333333333 +750 0.2916666666666667 0.1041666666666667 0.6041666666666666 +751 0.234375 0.08854166666666669 0.4999999999999999 +752 0.2916666666666666 0.9166666666666667 0.2083333333333334 +753 0.5 0.9375 0.1458333333333334 +754 0.7083333333333333 0.9166666666666667 0.2083333333333334 +755 0.4166666666666666 0.7916666666666667 0.2083333333333334 +756 0.5833333333333333 0.7916666666666667 0.2083333333333334 +757 0.4166666666666666 0.8541666666666667 0.1458333333333334 +758 0.5833333333333333 0.8541666666666667 0.1458333333333334 +759 0.3958333333333333 0.8958333333333334 0.2916666666666667 +760 0.6041666666666666 0.8958333333333334 0.2916666666666667 +761 0.3958333333333333 0.9270833333333334 0.1770833333333334 +762 0.6041666666666666 0.9270833333333334 0.1770833333333334 +763 0.4999999999999999 0.8229166666666667 0.1770833333333334 +764 0.4999999999999999 0.911458333333333 0.234375 +765 0.2916666666666666 0.7916666666666667 0.08333333333333336 +766 0.5 0.8541666666666667 0.0625 +767 0.2083333333333333 0.8958333333333334 0.1041666666666667 +768 0.5 0.9270833333333334 0.07291666666666669 +769 0.3958333333333333 0.8229166666666667 0.07291666666666669 +770 0.3541666666666666 0.911458333333333 0.08854166666666669 +771 0.7083333333333333 0.7916666666666667 0.08333333333333336 +772 0.7916666666666666 0.8958333333333334 0.1041666666666667 +773 0.6041666666666666 0.8229166666666667 0.07291666666666669 +774 0.6458333333333334 0.911458333333333 0.08854166666666669 +775 0.3958333333333333 0.7083333333333334 0.1041666666666667 +776 0.6041666666666666 0.7083333333333334 0.1041666666666667 +777 0.5000000000000001 0.7656249999999998 0.08854166666666669 +778 0.7083333333333333 0.2083333333333334 0.08333333333333336 +779 0.5 0.1458333333333334 0.0625 +780 0.7083333333333333 0.08333333333333336 0.2083333333333334 +781 0.5833333333333333 0.1458333333333334 0.1458333333333334 +782 0.5 0.0625 0.1458333333333334 +783 0.7916666666666666 0.1041666666666667 0.1041666666666667 +784 0.6041666666666666 0.1770833333333334 0.07291666666666669 +785 0.5 0.07291666666666669 0.07291666666666669 +786 0.6041666666666666 0.07291666666666669 0.1770833333333334 +787 0.6458333333333334 0.08854166666666669 0.08854166666666669 +788 0.2916666666666666 0.2083333333333334 0.08333333333333336 +789 0.5833333333333333 0.2083333333333334 0.2083333333333334 +790 0.4166666666666666 0.2083333333333334 0.2083333333333334 +791 0.4166666666666666 0.1458333333333334 0.1458333333333334 +792 0.6041666666666666 0.2916666666666667 0.1041666666666667 +793 0.3958333333333333 0.2916666666666667 0.1041666666666667 +794 0.3958333333333333 0.1770833333333334 0.07291666666666669 +795 0.4999999999999999 0.1770833333333334 0.1770833333333334 +796 0.5000000000000001 0.2343749999999999 0.08854166666666669 +797 0.2916666666666666 0.08333333333333336 0.2083333333333334 +798 0.2083333333333333 0.1041666666666667 0.1041666666666667 +799 0.3958333333333333 0.07291666666666669 0.1770833333333334 +800 0.3541666666666667 0.08854166666666669 0.08854166666666669 +801 0.6041666666666666 0.1041666666666667 0.2916666666666667 +802 0.3958333333333333 0.1041666666666667 0.2916666666666667 +803 0.5 0.08854166666666669 0.234375 +804 0.9375 0.1458333333333334 0.5 +805 0.9166666666666667 0.2083333333333334 0.2916666666666666 +806 0.7916666666666667 0.08333333333333336 0.2916666666666666 +807 0.8541666666666667 0.0625 0.5 +808 0.8541666666666667 0.1458333333333334 0.4166666666666666 +809 0.8958333333333334 0.1041666666666667 0.2083333333333333 +810 0.9270833333333334 0.07291666666666669 0.5 +811 0.9270833333333334 0.1770833333333334 0.3958333333333333 +812 0.8229166666666667 0.07291666666666669 0.3958333333333333 +813 0.911458333333333 0.08854166666666669 0.3541666666666667 +814 0.9166666666666667 0.2083333333333334 0.7083333333333333 +815 0.7916666666666667 0.08333333333333336 0.7083333333333333 +816 0.8541666666666667 0.1458333333333334 0.5833333333333333 +817 0.8958333333333334 0.1041666666666667 0.7916666666666666 +818 0.9270833333333334 0.1770833333333334 0.6041666666666666 +819 0.8229166666666667 0.07291666666666669 0.6041666666666666 +820 0.911458333333333 0.08854166666666669 0.6458333333333335 +821 0.7916666666666667 0.2083333333333334 0.4166666666666666 +822 0.7916666666666667 0.2083333333333334 0.5833333333333333 +823 0.8958333333333334 0.2916666666666667 0.3958333333333333 +824 0.8958333333333334 0.2916666666666667 0.6041666666666666 +825 0.8229166666666667 0.1770833333333334 0.4999999999999999 +826 0.911458333333333 0.234375 0.5 +827 0.7083333333333334 0.1041666666666667 0.6041666666666666 +828 0.7083333333333334 0.1041666666666667 0.3958333333333333 +829 0.7656249999999998 0.08854166666666669 0.5 +830 0.2083333333333334 0.2916666666666666 0.7916666666666667 +831 0.2916666666666666 0.2916666666666666 0.7083333333333333 +832 0.2083333333333334 0.2083333333333334 0.7083333333333333 +833 0.25 0.3541666666666666 0.75 +834 0.1458333333333334 0.25 0.75 +835 0.25 0.25 0.6458333333333333 +836 0.1770833333333333 0.3229166666666665 0.6770833333333337 +837 0.2916666666666666 0.2083333333333334 0.7916666666666667 +838 0.25 0.25 0.8541666666666667 +839 0.3541666666666666 0.25 0.75 +840 0.3229166666666669 0.3229166666666665 0.8229166666666665 +841 0.25 0.1458333333333334 0.75 +842 0.1770833333333334 0.1770833333333335 0.8229166666666666 +843 0.3229166666666665 0.1770833333333333 0.6770833333333335 +844 0.2083333333333334 0.7916666666666667 0.7083333333333333 +845 0.2083333333333334 0.7083333333333333 0.7916666666666667 +846 0.2916666666666666 0.7916666666666667 0.7916666666666667 +847 0.1458333333333334 0.75 0.75 +848 0.25 0.8541666666666667 0.75 +849 0.25 0.75 0.8541666666666667 +850 0.1770833333333334 0.8229166666666667 0.8229166666666663 +851 0.2916666666666666 0.7083333333333333 0.7083333333333333 +852 0.25 0.75 0.6458333333333333 +853 0.25 0.6458333333333333 0.75 +854 0.1770833333333334 0.6770833333333337 0.6770833333333335 +855 0.3541666666666666 0.75 0.75 +856 0.3229166666666669 0.8229166666666669 0.6770833333333334 +857 0.3229166666666669 0.677083333333333 0.8229166666666667 +858 0.7083333333333333 0.7916666666666667 0.7916666666666667 +859 0.7083333333333333 0.7083333333333333 0.7083333333333333 +860 0.7916666666666667 0.7083333333333333 0.7916666666666667 +861 0.6458333333333333 0.75 0.75 +862 0.75 0.75 0.8541666666666667 +863 0.75 0.6458333333333333 0.75 +864 0.677083333333333 0.6770833333333337 0.8229166666666667 +865 0.7916666666666667 0.7916666666666667 0.7083333333333333 +866 0.75 0.8541666666666667 0.75 +867 0.75 0.75 0.6458333333333333 +868 0.677083333333333 0.8229166666666665 0.6770833333333337 +869 0.8541666666666667 0.75 0.75 +870 0.8229166666666663 0.8229166666666666 0.8229166666666669 +871 0.8229166666666663 0.6770833333333335 0.677083333333333 +872 0.2083333333333334 0.7916666666666667 0.2916666666666666 +873 0.2916666666666666 0.7916666666666667 0.2083333333333334 +874 0.2083333333333334 0.7083333333333333 0.2083333333333334 +875 0.25 0.8541666666666667 0.25 +876 0.1458333333333334 0.75 0.25 +877 0.25 0.75 0.1458333333333334 +878 0.1770833333333333 0.822916666666667 0.1770833333333333 +879 0.2916666666666666 0.7083333333333333 0.2916666666666666 +880 0.25 0.75 0.3541666666666666 +881 0.3541666666666666 0.75 0.25 +882 0.3229166666666669 0.822916666666667 0.3229166666666667 +883 0.25 0.6458333333333333 0.25 +884 0.1770833333333334 0.6770833333333337 0.3229166666666667 +885 0.3229166666666665 0.6770833333333337 0.1770833333333333 +886 0.7916666666666667 0.7916666666666667 0.2916666666666666 +887 0.7916666666666667 0.7083333333333333 0.2083333333333334 +888 0.7083333333333333 0.7916666666666667 0.2083333333333334 +889 0.8541666666666667 0.75 0.25 +890 0.75 0.8541666666666667 0.25 +891 0.75 0.75 0.1458333333333334 +892 0.822916666666667 0.8229166666666667 0.1770833333333333 +893 0.7083333333333333 0.7083333333333333 0.2916666666666666 +894 0.75 0.75 0.3541666666666666 +895 0.75 0.6458333333333333 0.25 +896 0.822916666666667 0.6770833333333337 0.3229166666666667 +897 0.6458333333333333 0.75 0.25 +898 0.6770833333333337 0.8229166666666669 0.3229166666666667 +899 0.6770833333333337 0.677083333333333 0.1770833333333333 +900 0.7916666666666667 0.2916666666666666 0.2083333333333334 +901 0.7083333333333333 0.2916666666666666 0.2916666666666666 +902 0.7916666666666667 0.2083333333333334 0.2916666666666666 +903 0.75 0.3541666666666666 0.25 +904 0.8541666666666667 0.25 0.25 +905 0.75 0.25 0.3541666666666666 +906 0.8229166666666667 0.3229166666666665 0.3229166666666669 +907 0.7083333333333333 0.2083333333333334 0.2083333333333334 +908 0.75 0.25 0.1458333333333334 +909 0.6458333333333333 0.25 0.25 +910 0.6770833333333337 0.3229166666666665 0.1770833333333333 +911 0.75 0.1458333333333334 0.25 +912 0.8229166666666669 0.1770833333333335 0.1770833333333333 +913 0.677083333333333 0.1770833333333333 0.3229166666666667 +914 0.7916666666666667 0.2916666666666666 0.7916666666666667 +915 0.7916666666666667 0.2083333333333334 0.7083333333333333 +916 0.7083333333333333 0.2916666666666666 0.7083333333333333 +917 0.8541666666666667 0.25 0.75 +918 0.75 0.3541666666666666 0.75 +919 0.75 0.25 0.6458333333333333 +920 0.822916666666667 0.3229166666666667 0.6770833333333337 +921 0.7083333333333333 0.2083333333333334 0.7916666666666667 +922 0.75 0.25 0.8541666666666667 +923 0.75 0.1458333333333334 0.75 +924 0.822916666666667 0.1770833333333334 0.8229166666666665 +925 0.6458333333333333 0.25 0.75 +926 0.6770833333333337 0.3229166666666666 0.8229166666666666 +927 0.6770833333333337 0.1770833333333334 0.6770833333333335 +928 0.2083333333333334 0.2916666666666666 0.2083333333333334 +929 0.2083333333333334 0.2083333333333334 0.2916666666666666 +930 0.2916666666666666 0.2916666666666666 0.2916666666666666 +931 0.1458333333333334 0.25 0.25 +932 0.25 0.3541666666666666 0.25 +933 0.25 0.25 0.3541666666666666 +934 0.1770833333333334 0.3229166666666667 0.3229166666666669 +935 0.2916666666666666 0.2083333333333334 0.2083333333333334 +936 0.25 0.25 0.1458333333333334 +937 0.25 0.1458333333333334 0.25 +938 0.1770833333333334 0.1770833333333334 0.1770833333333333 +939 0.3541666666666666 0.25 0.25 +940 0.3229166666666669 0.3229166666666666 0.1770833333333333 +941 0.3229166666666669 0.1770833333333334 0.3229166666666667 +$EndNodes +$Elements +1112 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 15 2 0 7 7 +8 15 2 0 8 8 +9 1 2 0 1 2 10 +10 1 2 0 1 10 9 +11 1 2 0 1 9 11 +12 1 2 0 1 11 1 +13 1 2 0 2 1 13 +14 1 2 0 2 13 12 +15 1 2 0 2 12 14 +16 1 2 0 2 14 3 +17 1 2 0 3 4 16 +18 1 2 0 3 16 15 +19 1 2 0 3 15 17 +20 1 2 0 3 17 3 +21 1 2 0 4 2 19 +22 1 2 0 4 19 18 +23 1 2 0 4 18 20 +24 1 2 0 4 20 4 +25 1 2 0 5 6 22 +26 1 2 0 5 22 21 +27 1 2 0 5 21 23 +28 1 2 0 5 23 5 +29 1 2 0 6 5 25 +30 1 2 0 6 25 24 +31 1 2 0 6 24 26 +32 1 2 0 6 26 7 +33 1 2 0 7 8 28 +34 1 2 0 7 28 27 +35 1 2 0 7 27 29 +36 1 2 0 7 29 7 +37 1 2 0 8 6 31 +38 1 2 0 8 31 30 +39 1 2 0 8 30 32 +40 1 2 0 8 32 8 +41 1 2 0 9 2 34 +42 1 2 0 9 34 33 +43 1 2 0 9 33 35 +44 1 2 0 9 35 6 +45 1 2 0 10 1 37 +46 1 2 0 10 37 36 +47 1 2 0 10 36 38 +48 1 2 0 10 38 5 +49 1 2 0 11 4 40 +50 1 2 0 11 40 39 +51 1 2 0 11 39 41 +52 1 2 0 11 41 8 +53 1 2 0 12 3 43 +54 1 2 0 12 43 42 +55 1 2 0 12 42 44 +56 1 2 0 12 44 7 +57 3 2 0 1 2 10 57 56 +58 3 2 0 1 10 9 54 57 +59 3 2 0 1 57 54 50 55 +60 3 2 0 1 56 57 55 47 +61 3 2 0 1 9 11 60 54 +62 3 2 0 1 11 1 58 60 +63 3 2 0 1 60 58 46 59 +64 3 2 0 1 54 60 59 50 +65 3 2 0 1 47 55 63 62 +66 3 2 0 1 55 50 59 63 +67 3 2 0 1 63 59 46 61 +68 3 2 0 1 62 63 61 45 +69 3 2 0 1 1 13 66 58 +70 3 2 0 1 13 12 64 66 +71 3 2 0 1 66 64 51 65 +72 3 2 0 1 58 66 65 46 +73 3 2 0 1 12 14 69 64 +74 3 2 0 1 14 3 67 69 +75 3 2 0 1 69 67 48 68 +76 3 2 0 1 64 69 68 51 +77 3 2 0 1 46 65 71 61 +78 3 2 0 1 65 51 68 71 +79 3 2 0 1 71 68 48 70 +80 3 2 0 1 61 71 70 45 +81 3 2 0 1 4 20 75 74 +82 3 2 0 1 20 18 72 75 +83 3 2 0 1 75 72 52 73 +84 3 2 0 1 74 75 73 49 +85 3 2 0 1 18 19 77 72 +86 3 2 0 1 19 2 56 77 +87 3 2 0 1 77 56 47 76 +88 3 2 0 1 72 77 76 52 +89 3 2 0 1 49 73 79 78 +90 3 2 0 1 73 52 76 79 +91 3 2 0 1 79 76 47 62 +92 3 2 0 1 78 79 62 45 +93 3 2 0 1 3 17 82 67 +94 3 2 0 1 17 15 80 82 +95 3 2 0 1 82 80 53 81 +96 3 2 0 1 67 82 81 48 +97 3 2 0 1 15 16 84 80 +98 3 2 0 1 16 4 74 84 +99 3 2 0 1 84 74 49 83 +100 3 2 0 1 80 84 83 53 +101 3 2 0 1 48 81 85 70 +102 3 2 0 1 81 53 83 85 +103 3 2 0 1 85 83 49 78 +104 3 2 0 1 70 85 78 45 +105 3 2 0 2 6 95 98 22 +106 3 2 0 2 95 87 96 98 +107 3 2 0 2 98 96 91 97 +108 3 2 0 2 22 98 97 21 +109 3 2 0 2 87 99 102 96 +110 3 2 0 2 99 86 100 102 +111 3 2 0 2 102 100 88 101 +112 3 2 0 2 96 102 101 91 +113 3 2 0 2 21 97 104 23 +114 3 2 0 2 97 91 101 104 +115 3 2 0 2 104 101 88 103 +116 3 2 0 2 23 104 103 5 +117 3 2 0 2 5 103 107 25 +118 3 2 0 2 103 88 105 107 +119 3 2 0 2 107 105 92 106 +120 3 2 0 2 25 107 106 24 +121 3 2 0 2 88 100 110 105 +122 3 2 0 2 100 86 108 110 +123 3 2 0 2 110 108 89 109 +124 3 2 0 2 105 110 109 92 +125 3 2 0 2 24 106 112 26 +126 3 2 0 2 106 92 109 112 +127 3 2 0 2 112 109 89 111 +128 3 2 0 2 26 112 111 7 +129 3 2 0 2 8 113 116 32 +130 3 2 0 2 113 90 114 116 +131 3 2 0 2 116 114 93 115 +132 3 2 0 2 32 116 115 30 +133 3 2 0 2 90 117 119 114 +134 3 2 0 2 117 86 99 119 +135 3 2 0 2 119 99 87 118 +136 3 2 0 2 114 119 118 93 +137 3 2 0 2 30 115 120 31 +138 3 2 0 2 115 93 118 120 +139 3 2 0 2 120 118 87 95 +140 3 2 0 2 31 120 95 6 +141 3 2 0 2 7 111 123 29 +142 3 2 0 2 111 89 121 123 +143 3 2 0 2 123 121 94 122 +144 3 2 0 2 29 123 122 27 +145 3 2 0 2 89 108 125 121 +146 3 2 0 2 108 86 117 125 +147 3 2 0 2 125 117 90 124 +148 3 2 0 2 121 125 124 94 +149 3 2 0 2 27 122 126 28 +150 3 2 0 2 122 94 124 126 +151 3 2 0 2 126 124 90 113 +152 3 2 0 2 28 126 113 8 +153 3 2 0 3 1 11 139 138 +154 3 2 0 3 11 9 136 139 +155 3 2 0 3 139 136 132 137 +156 3 2 0 3 138 139 137 129 +157 3 2 0 3 9 10 142 136 +158 3 2 0 3 10 2 140 142 +159 3 2 0 3 142 140 128 141 +160 3 2 0 3 136 142 141 132 +161 3 2 0 3 129 137 145 144 +162 3 2 0 3 137 132 141 145 +163 3 2 0 3 145 141 128 143 +164 3 2 0 3 144 145 143 127 +165 3 2 0 3 5 38 149 148 +166 3 2 0 3 38 36 146 149 +167 3 2 0 3 149 146 133 147 +168 3 2 0 3 148 149 147 130 +169 3 2 0 3 36 37 151 146 +170 3 2 0 3 37 1 138 151 +171 3 2 0 3 151 138 129 150 +172 3 2 0 3 146 151 150 133 +173 3 2 0 3 130 147 153 152 +174 3 2 0 3 147 133 150 153 +175 3 2 0 3 153 150 129 144 +176 3 2 0 3 152 153 144 127 +177 3 2 0 3 2 34 156 140 +178 3 2 0 3 34 33 154 156 +179 3 2 0 3 156 154 134 155 +180 3 2 0 3 140 156 155 128 +181 3 2 0 3 33 35 159 154 +182 3 2 0 3 35 6 157 159 +183 3 2 0 3 159 157 131 158 +184 3 2 0 3 154 159 158 134 +185 3 2 0 3 128 155 161 143 +186 3 2 0 3 155 134 158 161 +187 3 2 0 3 161 158 131 160 +188 3 2 0 3 143 161 160 127 +189 3 2 0 3 6 22 164 157 +190 3 2 0 3 22 21 162 164 +191 3 2 0 3 164 162 135 163 +192 3 2 0 3 157 164 163 131 +193 3 2 0 3 21 23 166 162 +194 3 2 0 3 23 5 148 166 +195 3 2 0 3 166 148 130 165 +196 3 2 0 3 162 166 165 135 +197 3 2 0 3 131 163 167 160 +198 3 2 0 3 163 135 165 167 +199 3 2 0 3 167 165 130 152 +200 3 2 0 3 160 167 152 127 +201 3 2 0 4 3 177 180 17 +202 3 2 0 4 177 169 178 180 +203 3 2 0 4 180 178 173 179 +204 3 2 0 4 17 180 179 15 +205 3 2 0 4 169 181 184 178 +206 3 2 0 4 181 168 182 184 +207 3 2 0 4 184 182 170 183 +208 3 2 0 4 178 184 183 173 +209 3 2 0 4 15 179 186 16 +210 3 2 0 4 179 173 183 186 +211 3 2 0 4 186 183 170 185 +212 3 2 0 4 16 186 185 4 +213 3 2 0 4 7 187 190 44 +214 3 2 0 4 187 171 188 190 +215 3 2 0 4 190 188 174 189 +216 3 2 0 4 44 190 189 42 +217 3 2 0 4 171 191 193 188 +218 3 2 0 4 191 168 181 193 +219 3 2 0 4 193 181 169 192 +220 3 2 0 4 188 193 192 174 +221 3 2 0 4 42 189 194 43 +222 3 2 0 4 189 174 192 194 +223 3 2 0 4 194 192 169 177 +224 3 2 0 4 43 194 177 3 +225 3 2 0 4 4 185 197 40 +226 3 2 0 4 185 170 195 197 +227 3 2 0 4 197 195 175 196 +228 3 2 0 4 40 197 196 39 +229 3 2 0 4 170 182 200 195 +230 3 2 0 4 182 168 198 200 +231 3 2 0 4 200 198 172 199 +232 3 2 0 4 195 200 199 175 +233 3 2 0 4 39 196 202 41 +234 3 2 0 4 196 175 199 202 +235 3 2 0 4 202 199 172 201 +236 3 2 0 4 41 202 201 8 +237 3 2 0 4 8 201 205 28 +238 3 2 0 4 201 172 203 205 +239 3 2 0 4 205 203 176 204 +240 3 2 0 4 28 205 204 27 +241 3 2 0 4 172 198 207 203 +242 3 2 0 4 198 168 191 207 +243 3 2 0 4 207 191 171 206 +244 3 2 0 4 203 207 206 176 +245 3 2 0 4 27 204 208 29 +246 3 2 0 4 204 176 206 208 +247 3 2 0 4 208 206 171 187 +248 3 2 0 4 29 208 187 7 +249 3 2 0 5 2 19 221 220 +250 3 2 0 5 19 18 218 221 +251 3 2 0 5 221 218 214 219 +252 3 2 0 5 220 221 219 211 +253 3 2 0 5 18 20 224 218 +254 3 2 0 5 20 4 222 224 +255 3 2 0 5 224 222 210 223 +256 3 2 0 5 218 224 223 214 +257 3 2 0 5 211 219 227 226 +258 3 2 0 5 219 214 223 227 +259 3 2 0 5 227 223 210 225 +260 3 2 0 5 226 227 225 209 +261 3 2 0 5 6 35 231 230 +262 3 2 0 5 35 33 228 231 +263 3 2 0 5 231 228 215 229 +264 3 2 0 5 230 231 229 212 +265 3 2 0 5 33 34 233 228 +266 3 2 0 5 34 2 220 233 +267 3 2 0 5 233 220 211 232 +268 3 2 0 5 228 233 232 215 +269 3 2 0 5 212 229 235 234 +270 3 2 0 5 229 215 232 235 +271 3 2 0 5 235 232 211 226 +272 3 2 0 5 234 235 226 209 +273 3 2 0 5 4 40 238 222 +274 3 2 0 5 40 39 236 238 +275 3 2 0 5 238 236 216 237 +276 3 2 0 5 222 238 237 210 +277 3 2 0 5 39 41 241 236 +278 3 2 0 5 41 8 239 241 +279 3 2 0 5 241 239 213 240 +280 3 2 0 5 236 241 240 216 +281 3 2 0 5 210 237 243 225 +282 3 2 0 5 237 216 240 243 +283 3 2 0 5 243 240 213 242 +284 3 2 0 5 225 243 242 209 +285 3 2 0 5 8 32 246 239 +286 3 2 0 5 32 30 244 246 +287 3 2 0 5 246 244 217 245 +288 3 2 0 5 239 246 245 213 +289 3 2 0 5 30 31 248 244 +290 3 2 0 5 31 6 230 248 +291 3 2 0 5 248 230 212 247 +292 3 2 0 5 244 248 247 217 +293 3 2 0 5 213 245 249 242 +294 3 2 0 5 245 217 247 249 +295 3 2 0 5 249 247 212 234 +296 3 2 0 5 242 249 234 209 +297 3 2 0 6 1 259 262 13 +298 3 2 0 6 259 251 260 262 +299 3 2 0 6 262 260 255 261 +300 3 2 0 6 13 262 261 12 +301 3 2 0 6 251 263 266 260 +302 3 2 0 6 263 250 264 266 +303 3 2 0 6 266 264 252 265 +304 3 2 0 6 260 266 265 255 +305 3 2 0 6 12 261 268 14 +306 3 2 0 6 261 255 265 268 +307 3 2 0 6 268 265 252 267 +308 3 2 0 6 14 268 267 3 +309 3 2 0 6 5 269 272 38 +310 3 2 0 6 269 253 270 272 +311 3 2 0 6 272 270 256 271 +312 3 2 0 6 38 272 271 36 +313 3 2 0 6 253 273 275 270 +314 3 2 0 6 273 250 263 275 +315 3 2 0 6 275 263 251 274 +316 3 2 0 6 270 275 274 256 +317 3 2 0 6 36 271 276 37 +318 3 2 0 6 271 256 274 276 +319 3 2 0 6 276 274 251 259 +320 3 2 0 6 37 276 259 1 +321 3 2 0 6 3 267 279 43 +322 3 2 0 6 267 252 277 279 +323 3 2 0 6 279 277 257 278 +324 3 2 0 6 43 279 278 42 +325 3 2 0 6 252 264 282 277 +326 3 2 0 6 264 250 280 282 +327 3 2 0 6 282 280 254 281 +328 3 2 0 6 277 282 281 257 +329 3 2 0 6 42 278 284 44 +330 3 2 0 6 278 257 281 284 +331 3 2 0 6 284 281 254 283 +332 3 2 0 6 44 284 283 7 +333 3 2 0 6 7 283 287 26 +334 3 2 0 6 283 254 285 287 +335 3 2 0 6 287 285 258 286 +336 3 2 0 6 26 287 286 24 +337 3 2 0 6 254 280 289 285 +338 3 2 0 6 280 250 273 289 +339 3 2 0 6 289 273 253 288 +340 3 2 0 6 285 289 288 258 +341 3 2 0 6 24 286 290 25 +342 3 2 0 6 286 258 288 290 +343 3 2 0 6 290 288 253 269 +344 3 2 0 6 25 290 269 5 +345 5 2 0 1 209 364 376 365 366 377 382 378 +346 5 2 0 1 366 377 382 378 294 372 381 373 +347 5 2 0 1 364 291 367 376 377 368 379 382 +348 5 2 0 1 377 368 379 382 372 305 374 381 +349 5 2 0 1 365 376 369 293 378 382 380 371 +350 5 2 0 1 378 382 380 371 373 381 375 306 +351 5 2 0 1 376 367 304 369 382 379 370 380 +352 5 2 0 1 382 379 370 380 381 374 308 375 +353 5 2 0 1 291 383 391 367 368 392 396 379 +354 5 2 0 1 368 392 396 379 305 388 395 374 +355 5 2 0 1 383 86 384 391 392 385 393 396 +356 5 2 0 1 392 385 393 396 388 296 389 395 +357 5 2 0 1 367 391 386 304 379 396 394 370 +358 5 2 0 1 379 396 394 370 374 395 390 308 +359 5 2 0 1 391 384 292 386 396 393 387 394 +360 5 2 0 1 396 393 387 394 395 389 307 390 +361 5 2 0 1 293 369 402 397 371 380 406 403 +362 5 2 0 1 371 380 406 403 306 375 405 400 +363 5 2 0 1 369 304 386 402 380 370 394 406 +364 5 2 0 1 380 370 394 406 375 308 390 405 +365 5 2 0 1 397 402 398 45 403 406 404 399 +366 5 2 0 1 403 406 404 399 400 405 401 295 +367 5 2 0 1 402 386 292 398 406 394 387 404 +368 5 2 0 1 406 394 387 404 405 390 307 401 +369 5 2 0 1 127 407 410 408 409 411 413 412 +370 5 2 0 1 409 411 413 412 295 401 405 400 +371 5 2 0 1 407 296 388 410 411 389 395 413 +372 5 2 0 1 411 389 395 413 401 307 390 405 +373 5 2 0 1 408 410 372 294 412 413 381 373 +374 5 2 0 1 412 413 381 373 400 405 375 306 +375 5 2 0 1 410 388 305 372 413 395 374 381 +376 5 2 0 1 413 395 374 381 405 390 308 375 +377 5 2 0 1 250 414 426 415 416 427 432 428 +378 5 2 0 1 416 427 432 428 300 422 431 423 +379 5 2 0 1 414 297 417 426 427 418 429 432 +380 5 2 0 1 427 418 429 432 422 310 424 431 +381 5 2 0 1 415 426 419 299 428 432 430 421 +382 5 2 0 1 428 432 430 421 423 431 425 311 +383 5 2 0 1 426 417 309 419 432 429 420 430 +384 5 2 0 1 432 429 420 430 431 424 313 425 +385 5 2 0 1 297 433 440 417 418 441 445 429 +386 5 2 0 1 418 441 445 429 310 437 444 424 +387 5 2 0 1 433 45 434 440 441 398 442 445 +388 5 2 0 1 441 398 442 445 437 292 438 444 +389 5 2 0 1 417 440 435 309 429 445 443 420 +390 5 2 0 1 429 445 443 420 424 444 439 313 +391 5 2 0 1 440 434 298 435 445 442 436 443 +392 5 2 0 1 445 442 436 443 444 438 312 439 +393 5 2 0 1 299 419 451 446 421 430 455 452 +394 5 2 0 1 421 430 455 452 311 425 454 449 +395 5 2 0 1 419 309 435 451 430 420 443 455 +396 5 2 0 1 430 420 443 455 425 313 439 454 +397 5 2 0 1 446 451 447 168 452 455 453 448 +398 5 2 0 1 452 455 453 448 449 454 450 301 +399 5 2 0 1 451 435 298 447 455 443 436 453 +400 5 2 0 1 455 443 436 453 454 439 312 450 +401 5 2 0 1 86 384 458 456 457 459 461 460 +402 5 2 0 1 457 459 461 460 301 450 454 449 +403 5 2 0 1 384 292 437 458 459 438 444 461 +404 5 2 0 1 459 438 444 461 450 312 439 454 +405 5 2 0 1 456 458 422 300 460 461 431 423 +406 5 2 0 1 460 461 431 423 449 454 425 311 +407 5 2 0 1 458 437 310 422 461 444 424 431 +408 5 2 0 1 461 444 424 431 454 439 313 425 +409 5 2 0 1 86 384 459 457 383 391 471 467 +410 5 2 0 1 383 391 471 467 291 367 470 464 +411 5 2 0 1 384 292 438 459 391 386 468 471 +412 5 2 0 1 391 386 468 471 367 304 465 470 +413 5 2 0 1 457 459 450 301 467 471 469 463 +414 5 2 0 1 467 471 469 463 464 470 466 314 +415 5 2 0 1 459 438 312 450 471 468 462 469 +416 5 2 0 1 471 468 462 469 470 465 316 466 +417 5 2 0 1 292 398 442 438 386 402 478 468 +418 5 2 0 1 386 402 478 468 304 369 477 465 +419 5 2 0 1 398 45 434 442 402 397 475 478 +420 5 2 0 1 402 397 475 478 369 293 473 477 +421 5 2 0 1 438 442 436 312 468 478 476 462 +422 5 2 0 1 468 478 476 462 465 477 474 316 +423 5 2 0 1 442 434 298 436 478 475 472 476 +424 5 2 0 1 478 475 472 476 477 473 315 474 +425 5 2 0 1 301 450 453 448 463 469 485 482 +426 5 2 0 1 463 469 485 482 314 466 484 480 +427 5 2 0 1 450 312 436 453 469 462 476 485 +428 5 2 0 1 469 462 476 485 466 316 474 484 +429 5 2 0 1 448 453 447 168 482 485 483 479 +430 5 2 0 1 482 485 483 479 480 484 481 302 +431 5 2 0 1 453 436 298 447 485 476 472 483 +432 5 2 0 1 485 476 472 483 484 474 315 481 +433 5 2 0 1 209 365 376 364 486 487 489 488 +434 5 2 0 1 486 487 489 488 302 481 484 480 +435 5 2 0 1 365 293 369 376 487 473 477 489 +436 5 2 0 1 487 473 477 489 481 315 474 484 +437 5 2 0 1 364 376 367 291 488 489 470 464 +438 5 2 0 1 488 489 470 464 480 484 466 314 +439 5 2 0 1 376 369 304 367 489 477 465 470 +440 5 2 0 1 489 477 465 470 484 474 316 466 +441 5 2 0 1 250 414 427 416 490 498 503 499 +442 5 2 0 1 490 498 503 499 303 494 502 495 +443 5 2 0 1 414 297 418 427 498 491 500 503 +444 5 2 0 1 498 491 500 503 494 317 496 502 +445 5 2 0 1 416 427 422 300 499 503 501 493 +446 5 2 0 1 499 503 501 493 495 502 497 318 +447 5 2 0 1 427 418 310 422 503 500 492 501 +448 5 2 0 1 503 500 492 501 502 496 319 497 +449 5 2 0 1 297 433 441 418 491 506 509 500 +450 5 2 0 1 491 506 509 500 317 504 508 496 +451 5 2 0 1 433 45 398 441 506 399 404 509 +452 5 2 0 1 506 399 404 509 504 295 401 508 +453 5 2 0 1 418 441 437 310 500 509 507 492 +454 5 2 0 1 500 509 507 492 496 508 505 319 +455 5 2 0 1 441 398 292 437 509 404 387 507 +456 5 2 0 1 509 404 387 507 508 401 307 505 +457 5 2 0 1 300 422 458 456 493 501 513 511 +458 5 2 0 1 493 501 513 511 318 497 512 510 +459 5 2 0 1 422 310 437 458 501 492 507 513 +460 5 2 0 1 501 492 507 513 497 319 505 512 +461 5 2 0 1 456 458 384 86 511 513 393 385 +462 5 2 0 1 511 513 393 385 510 512 389 296 +463 5 2 0 1 458 437 292 384 513 507 387 393 +464 5 2 0 1 513 507 387 393 512 505 307 389 +465 5 2 0 1 127 409 515 514 407 411 517 516 +466 5 2 0 1 407 411 517 516 296 389 512 510 +467 5 2 0 1 409 295 504 515 411 401 508 517 +468 5 2 0 1 411 401 508 517 389 307 505 512 +469 5 2 0 1 514 515 494 303 516 517 502 495 +470 5 2 0 1 516 517 502 495 510 512 497 318 +471 5 2 0 1 515 504 317 494 517 508 496 502 +472 5 2 0 1 517 508 496 502 512 505 319 497 +473 5 2 0 1 168 191 525 448 198 207 530 526 +474 5 2 0 1 198 207 530 526 172 203 529 522 +475 5 2 0 1 191 171 518 525 207 206 527 530 +476 5 2 0 1 207 206 527 530 203 176 523 529 +477 5 2 0 1 448 525 519 301 526 530 528 521 +478 5 2 0 1 526 530 528 521 522 529 524 321 +479 5 2 0 1 525 518 320 519 530 527 520 528 +480 5 2 0 1 530 527 520 528 529 523 322 524 +481 5 2 0 1 171 187 533 518 206 208 536 527 +482 5 2 0 1 206 208 536 527 176 204 535 523 +483 5 2 0 1 187 7 111 533 208 29 123 536 +484 5 2 0 1 208 29 123 536 204 27 122 535 +485 5 2 0 1 518 533 531 320 527 536 534 520 +486 5 2 0 1 527 536 534 520 523 535 532 322 +487 5 2 0 1 533 111 89 531 536 123 121 534 +488 5 2 0 1 536 123 121 534 535 122 94 532 +489 5 2 0 1 301 519 538 457 521 528 541 539 +490 5 2 0 1 521 528 541 539 321 524 540 537 +491 5 2 0 1 519 320 531 538 528 520 534 541 +492 5 2 0 1 528 520 534 541 524 322 532 540 +493 5 2 0 1 457 538 108 86 539 541 125 117 +494 5 2 0 1 539 541 125 117 537 540 124 90 +495 5 2 0 1 538 531 89 108 541 534 121 125 +496 5 2 0 1 541 534 121 125 540 532 94 124 +497 5 2 0 1 8 28 205 201 113 126 543 542 +498 5 2 0 1 113 126 543 542 90 124 540 537 +499 5 2 0 1 28 27 204 205 126 122 535 543 +500 5 2 0 1 126 122 535 543 124 94 532 540 +501 5 2 0 1 201 205 203 172 542 543 529 522 +502 5 2 0 1 542 543 529 522 537 540 524 321 +503 5 2 0 1 205 204 176 203 543 535 523 529 +504 5 2 0 1 543 535 523 529 540 532 322 524 +505 5 2 0 1 45 434 551 78 70 552 556 85 +506 5 2 0 1 70 552 556 85 48 548 555 81 +507 5 2 0 1 434 298 544 551 552 545 553 556 +508 5 2 0 1 552 545 553 556 548 324 549 555 +509 5 2 0 1 78 551 546 49 85 556 554 83 +510 5 2 0 1 85 556 554 83 81 555 550 53 +511 5 2 0 1 551 544 323 546 556 553 547 554 +512 5 2 0 1 556 553 547 554 555 549 325 550 +513 5 2 0 1 298 447 560 544 545 561 564 553 +514 5 2 0 1 545 561 564 553 324 558 563 549 +515 5 2 0 1 447 168 182 560 561 181 184 564 +516 5 2 0 1 561 181 184 564 558 169 178 563 +517 5 2 0 1 544 560 557 323 553 564 562 547 +518 5 2 0 1 553 564 562 547 549 563 559 325 +519 5 2 0 1 560 182 170 557 564 184 183 562 +520 5 2 0 1 564 184 183 562 563 178 173 559 +521 5 2 0 1 49 546 565 74 83 554 567 84 +522 5 2 0 1 83 554 567 84 53 550 566 80 +523 5 2 0 1 546 323 557 565 554 547 562 567 +524 5 2 0 1 554 547 562 567 550 325 559 566 +525 5 2 0 1 74 565 185 4 84 567 186 16 +526 5 2 0 1 84 567 186 16 80 566 179 15 +527 5 2 0 1 565 557 170 185 567 562 183 186 +528 5 2 0 1 567 562 183 186 566 559 173 179 +529 5 2 0 1 3 177 568 67 17 180 569 82 +530 5 2 0 1 17 180 569 82 15 179 566 80 +531 5 2 0 1 177 169 558 568 180 178 563 569 +532 5 2 0 1 180 178 563 569 179 173 559 566 +533 5 2 0 1 67 568 548 48 82 569 555 81 +534 5 2 0 1 82 569 555 81 80 566 550 53 +535 5 2 0 1 568 558 324 548 569 563 549 555 +536 5 2 0 1 569 563 549 555 566 559 325 550 +537 5 2 0 1 5 38 272 269 148 149 579 575 +538 5 2 0 1 148 149 579 575 130 147 578 572 +539 5 2 0 1 38 36 271 272 149 146 576 579 +540 5 2 0 1 149 146 576 579 147 133 573 578 +541 5 2 0 1 269 272 270 253 575 579 577 571 +542 5 2 0 1 575 579 577 571 572 578 574 326 +543 5 2 0 1 272 271 256 270 579 576 570 577 +544 5 2 0 1 579 576 570 577 578 573 328 574 +545 5 2 0 1 36 37 276 271 146 151 586 576 +546 5 2 0 1 146 151 586 576 133 150 585 573 +547 5 2 0 1 37 1 259 276 151 138 583 586 +548 5 2 0 1 151 138 583 586 150 129 581 585 +549 5 2 0 1 271 276 274 256 576 586 584 570 +550 5 2 0 1 576 586 584 570 573 585 582 328 +551 5 2 0 1 276 259 251 274 586 583 580 584 +552 5 2 0 1 586 583 580 584 585 581 327 582 +553 5 2 0 1 253 270 275 273 571 577 592 589 +554 5 2 0 1 571 577 592 589 326 574 591 587 +555 5 2 0 1 270 256 274 275 577 570 584 592 +556 5 2 0 1 577 570 584 592 574 328 582 591 +557 5 2 0 1 273 275 263 250 589 592 590 490 +558 5 2 0 1 589 592 590 490 587 591 588 303 +559 5 2 0 1 275 274 251 263 592 584 580 590 +560 5 2 0 1 592 584 580 590 591 582 327 588 +561 5 2 0 1 127 144 153 152 514 593 595 594 +562 5 2 0 1 514 593 595 594 303 588 591 587 +563 5 2 0 1 144 129 150 153 593 581 585 595 +564 5 2 0 1 593 581 585 595 588 327 582 591 +565 5 2 0 1 152 153 147 130 594 595 578 572 +566 5 2 0 1 594 595 578 572 587 591 574 326 +567 5 2 0 1 153 150 133 147 595 585 573 578 +568 5 2 0 1 595 585 573 578 591 582 328 574 +569 5 2 0 1 250 264 603 414 263 266 608 604 +570 5 2 0 1 263 266 608 604 251 260 607 600 +571 5 2 0 1 264 252 596 603 266 265 605 608 +572 5 2 0 1 266 265 605 608 260 255 601 607 +573 5 2 0 1 414 603 597 297 604 608 606 599 +574 5 2 0 1 604 608 606 599 600 607 602 330 +575 5 2 0 1 603 596 329 597 608 605 598 606 +576 5 2 0 1 608 605 598 606 607 601 331 602 +577 5 2 0 1 252 267 611 596 265 268 614 605 +578 5 2 0 1 265 268 614 605 255 261 613 601 +579 5 2 0 1 267 3 67 611 268 14 69 614 +580 5 2 0 1 268 14 69 614 261 12 64 613 +581 5 2 0 1 596 611 609 329 605 614 612 598 +582 5 2 0 1 605 614 612 598 601 613 610 331 +583 5 2 0 1 611 67 48 609 614 69 68 612 +584 5 2 0 1 614 69 68 612 613 64 51 610 +585 5 2 0 1 297 597 616 433 599 606 619 617 +586 5 2 0 1 599 606 619 617 330 602 618 615 +587 5 2 0 1 597 329 609 616 606 598 612 619 +588 5 2 0 1 606 598 612 619 602 331 610 618 +589 5 2 0 1 433 616 70 45 617 619 71 61 +590 5 2 0 1 617 619 71 61 615 618 65 46 +591 5 2 0 1 616 609 48 70 619 612 68 71 +592 5 2 0 1 619 612 68 71 618 610 51 65 +593 5 2 0 1 1 13 262 259 58 66 621 620 +594 5 2 0 1 58 66 621 620 46 65 618 615 +595 5 2 0 1 13 12 261 262 66 64 613 621 +596 5 2 0 1 66 64 613 621 65 51 610 618 +597 5 2 0 1 259 262 260 251 620 621 607 600 +598 5 2 0 1 620 621 607 600 615 618 602 330 +599 5 2 0 1 262 261 255 260 621 613 601 607 +600 5 2 0 1 621 613 601 607 618 610 331 602 +601 5 2 0 1 168 181 629 446 191 193 634 630 +602 5 2 0 1 191 193 634 630 171 188 633 626 +603 5 2 0 1 181 169 622 629 193 192 631 634 +604 5 2 0 1 193 192 631 634 188 174 627 633 +605 5 2 0 1 446 629 623 299 630 634 632 625 +606 5 2 0 1 630 634 632 625 626 633 628 333 +607 5 2 0 1 629 622 332 623 634 631 624 632 +608 5 2 0 1 634 631 624 632 633 627 334 628 +609 5 2 0 1 169 177 637 622 192 194 640 631 +610 5 2 0 1 192 194 640 631 174 189 639 627 +611 5 2 0 1 177 3 267 637 194 43 279 640 +612 5 2 0 1 194 43 279 640 189 42 278 639 +613 5 2 0 1 622 637 635 332 631 640 638 624 +614 5 2 0 1 631 640 638 624 627 639 636 334 +615 5 2 0 1 637 267 252 635 640 279 277 638 +616 5 2 0 1 640 279 277 638 639 278 257 636 +617 5 2 0 1 299 623 642 415 625 632 645 643 +618 5 2 0 1 625 632 645 643 333 628 644 641 +619 5 2 0 1 623 332 635 642 632 624 638 645 +620 5 2 0 1 632 624 638 645 628 334 636 644 +621 5 2 0 1 415 642 264 250 643 645 282 280 +622 5 2 0 1 643 645 282 280 641 644 281 254 +623 5 2 0 1 642 635 252 264 645 638 277 282 +624 5 2 0 1 645 638 277 282 644 636 257 281 +625 5 2 0 1 7 44 190 187 283 284 647 646 +626 5 2 0 1 283 284 647 646 254 281 644 641 +627 5 2 0 1 44 42 189 190 284 278 639 647 +628 5 2 0 1 284 278 639 647 281 257 636 644 +629 5 2 0 1 187 190 188 171 646 647 633 626 +630 5 2 0 1 646 647 633 626 641 644 628 333 +631 5 2 0 1 190 189 174 188 647 639 627 633 +632 5 2 0 1 647 639 627 633 644 636 334 628 +633 5 2 0 1 209 242 655 364 234 249 660 656 +634 5 2 0 1 234 249 660 656 212 247 659 652 +635 5 2 0 1 242 213 648 655 249 245 657 660 +636 5 2 0 1 249 245 657 660 247 217 653 659 +637 5 2 0 1 364 655 649 291 656 660 658 651 +638 5 2 0 1 656 660 658 651 652 659 654 336 +639 5 2 0 1 655 648 335 649 660 657 650 658 +640 5 2 0 1 660 657 650 658 659 653 337 654 +641 5 2 0 1 213 239 663 648 245 246 666 657 +642 5 2 0 1 245 246 666 657 217 244 665 653 +643 5 2 0 1 239 8 113 663 246 32 116 666 +644 5 2 0 1 246 32 116 666 244 30 115 665 +645 5 2 0 1 648 663 661 335 657 666 664 650 +646 5 2 0 1 657 666 664 650 653 665 662 337 +647 5 2 0 1 663 113 90 661 666 116 114 664 +648 5 2 0 1 666 116 114 664 665 115 93 662 +649 5 2 0 1 291 649 668 383 651 658 671 669 +650 5 2 0 1 651 658 671 669 336 654 670 667 +651 5 2 0 1 649 335 661 668 658 650 664 671 +652 5 2 0 1 658 650 664 671 654 337 662 670 +653 5 2 0 1 383 668 117 86 669 671 119 99 +654 5 2 0 1 669 671 119 99 667 670 118 87 +655 5 2 0 1 668 661 90 117 671 664 114 119 +656 5 2 0 1 671 664 114 119 670 662 93 118 +657 5 2 0 1 6 31 248 230 95 120 673 672 +658 5 2 0 1 95 120 673 672 87 118 670 667 +659 5 2 0 1 31 30 244 248 120 115 665 673 +660 5 2 0 1 120 115 665 673 118 93 662 670 +661 5 2 0 1 230 248 247 212 672 673 659 652 +662 5 2 0 1 672 673 659 652 667 670 654 336 +663 5 2 0 1 248 244 217 247 673 665 653 659 +664 5 2 0 1 673 665 653 659 670 662 337 654 +665 5 2 0 1 7 26 287 283 111 112 683 679 +666 5 2 0 1 111 112 683 679 89 109 682 676 +667 5 2 0 1 26 24 286 287 112 106 680 683 +668 5 2 0 1 112 106 680 683 109 92 677 682 +669 5 2 0 1 283 287 285 254 679 683 681 675 +670 5 2 0 1 679 683 681 675 676 682 678 338 +671 5 2 0 1 287 286 258 285 683 680 674 681 +672 5 2 0 1 683 680 674 681 682 677 340 678 +673 5 2 0 1 24 25 290 286 106 107 690 680 +674 5 2 0 1 106 107 690 680 92 105 689 677 +675 5 2 0 1 25 5 269 290 107 103 687 690 +676 5 2 0 1 107 103 687 690 105 88 685 689 +677 5 2 0 1 286 290 288 258 680 690 688 674 +678 5 2 0 1 680 690 688 674 677 689 686 340 +679 5 2 0 1 290 269 253 288 690 687 684 688 +680 5 2 0 1 690 687 684 688 689 685 339 686 +681 5 2 0 1 254 285 289 280 675 681 696 693 +682 5 2 0 1 675 681 696 693 338 678 695 691 +683 5 2 0 1 285 258 288 289 681 674 688 696 +684 5 2 0 1 681 674 688 696 678 340 686 695 +685 5 2 0 1 280 289 273 250 693 696 694 416 +686 5 2 0 1 693 696 694 416 691 695 692 300 +687 5 2 0 1 289 288 253 273 696 688 684 694 +688 5 2 0 1 696 688 684 694 695 686 339 692 +689 5 2 0 1 86 100 110 108 456 697 699 698 +690 5 2 0 1 456 697 699 698 300 692 695 691 +691 5 2 0 1 100 88 105 110 697 685 689 699 +692 5 2 0 1 697 685 689 699 692 339 686 695 +693 5 2 0 1 108 110 109 89 698 699 682 676 +694 5 2 0 1 698 699 682 676 691 695 678 338 +695 5 2 0 1 110 105 92 109 699 689 677 682 +696 5 2 0 1 699 689 677 682 695 686 340 678 +697 5 2 0 1 45 62 79 78 397 707 712 708 +698 5 2 0 1 397 707 712 708 293 703 711 704 +699 5 2 0 1 62 47 76 79 707 700 709 712 +700 5 2 0 1 707 700 709 712 703 341 705 711 +701 5 2 0 1 78 79 73 49 708 712 710 702 +702 5 2 0 1 708 712 710 702 704 711 706 342 +703 5 2 0 1 79 76 52 73 712 709 701 710 +704 5 2 0 1 712 709 701 710 711 705 343 706 +705 5 2 0 1 47 56 77 76 700 715 718 709 +706 5 2 0 1 700 715 718 709 341 713 717 705 +707 5 2 0 1 56 2 19 77 715 220 221 718 +708 5 2 0 1 715 220 221 718 713 211 219 717 +709 5 2 0 1 76 77 72 52 709 718 716 701 +710 5 2 0 1 709 718 716 701 705 717 714 343 +711 5 2 0 1 77 19 18 72 718 221 218 716 +712 5 2 0 1 718 221 218 716 717 219 214 714 +713 5 2 0 1 49 73 75 74 702 710 722 720 +714 5 2 0 1 702 710 722 720 342 706 721 719 +715 5 2 0 1 73 52 72 75 710 701 716 722 +716 5 2 0 1 710 701 716 722 706 343 714 721 +717 5 2 0 1 74 75 20 4 720 722 224 222 +718 5 2 0 1 720 722 224 222 719 721 223 210 +719 5 2 0 1 75 72 18 20 722 716 218 224 +720 5 2 0 1 722 716 218 224 721 714 214 223 +721 5 2 0 1 209 226 723 365 225 227 725 724 +722 5 2 0 1 225 227 725 724 210 223 721 719 +723 5 2 0 1 226 211 713 723 227 219 717 725 +724 5 2 0 1 227 219 717 725 223 214 714 721 +725 5 2 0 1 365 723 703 293 724 725 711 704 +726 5 2 0 1 724 725 711 704 719 721 706 342 +727 5 2 0 1 723 713 341 703 725 717 705 711 +728 5 2 0 1 725 717 705 711 721 714 343 706 +729 5 2 0 1 1 11 60 58 138 139 735 731 +730 5 2 0 1 138 139 735 731 129 137 734 728 +731 5 2 0 1 11 9 54 60 139 136 732 735 +732 5 2 0 1 139 136 732 735 137 132 729 734 +733 5 2 0 1 58 60 59 46 731 735 733 727 +734 5 2 0 1 731 735 733 727 728 734 730 344 +735 5 2 0 1 60 54 50 59 735 732 726 733 +736 5 2 0 1 735 732 726 733 734 729 346 730 +737 5 2 0 1 9 10 57 54 136 142 742 732 +738 5 2 0 1 136 142 742 732 132 141 741 729 +739 5 2 0 1 10 2 56 57 142 140 739 742 +740 5 2 0 1 142 140 739 742 141 128 737 741 +741 5 2 0 1 54 57 55 50 732 742 740 726 +742 5 2 0 1 732 742 740 726 729 741 738 346 +743 5 2 0 1 57 56 47 55 742 739 736 740 +744 5 2 0 1 742 739 736 740 741 737 345 738 +745 5 2 0 1 46 59 63 61 727 733 748 745 +746 5 2 0 1 727 733 748 745 344 730 747 743 +747 5 2 0 1 59 50 55 63 733 726 740 748 +748 5 2 0 1 733 726 740 748 730 346 738 747 +749 5 2 0 1 61 63 62 45 745 748 746 399 +750 5 2 0 1 745 748 746 399 743 747 744 295 +751 5 2 0 1 63 55 47 62 748 740 736 746 +752 5 2 0 1 748 740 736 746 747 738 345 744 +753 5 2 0 1 127 143 145 144 409 749 751 750 +754 5 2 0 1 409 749 751 750 295 744 747 743 +755 5 2 0 1 143 128 141 145 749 737 741 751 +756 5 2 0 1 749 737 741 751 744 345 738 747 +757 5 2 0 1 144 145 137 129 750 751 734 728 +758 5 2 0 1 750 751 734 728 743 747 730 344 +759 5 2 0 1 145 141 132 137 751 741 729 734 +760 5 2 0 1 751 741 729 734 747 738 346 730 +761 5 2 0 1 168 182 200 198 479 759 764 760 +762 5 2 0 1 479 759 764 760 302 755 763 756 +763 5 2 0 1 182 170 195 200 759 752 761 764 +764 5 2 0 1 759 752 761 764 755 347 757 763 +765 5 2 0 1 198 200 199 172 760 764 762 754 +766 5 2 0 1 760 764 762 754 756 763 758 348 +767 5 2 0 1 200 195 175 199 764 761 753 762 +768 5 2 0 1 764 761 753 762 763 757 349 758 +769 5 2 0 1 170 185 197 195 752 767 770 761 +770 5 2 0 1 752 767 770 761 347 765 769 757 +771 5 2 0 1 185 4 40 197 767 222 238 770 +772 5 2 0 1 767 222 238 770 765 210 237 769 +773 5 2 0 1 195 197 196 175 761 770 768 753 +774 5 2 0 1 761 770 768 753 757 769 766 349 +775 5 2 0 1 197 40 39 196 770 238 236 768 +776 5 2 0 1 770 238 236 768 769 237 216 766 +777 5 2 0 1 172 199 202 201 754 762 774 772 +778 5 2 0 1 754 762 774 772 348 758 773 771 +779 5 2 0 1 199 175 196 202 762 753 768 774 +780 5 2 0 1 762 753 768 774 758 349 766 773 +781 5 2 0 1 201 202 41 8 772 774 241 239 +782 5 2 0 1 772 774 241 239 771 773 240 213 +783 5 2 0 1 202 196 39 41 774 768 236 241 +784 5 2 0 1 774 768 236 241 773 766 216 240 +785 5 2 0 1 209 225 775 486 242 243 777 776 +786 5 2 0 1 242 243 777 776 213 240 773 771 +787 5 2 0 1 225 210 765 775 243 237 769 777 +788 5 2 0 1 243 237 769 777 240 216 766 773 +789 5 2 0 1 486 775 755 302 776 777 763 756 +790 5 2 0 1 776 777 763 756 771 773 758 348 +791 5 2 0 1 775 765 347 755 777 769 757 763 +792 5 2 0 1 777 769 757 763 773 766 349 758 +793 5 2 0 1 6 230 231 35 157 783 787 159 +794 5 2 0 1 157 783 787 159 131 780 786 158 +795 5 2 0 1 230 212 229 231 783 778 784 787 +796 5 2 0 1 783 778 784 787 780 350 781 786 +797 5 2 0 1 35 231 228 33 159 787 785 154 +798 5 2 0 1 159 787 785 154 158 786 782 134 +799 5 2 0 1 231 229 215 228 787 784 779 785 +800 5 2 0 1 787 784 779 785 786 781 352 782 +801 5 2 0 1 212 234 235 229 778 792 796 784 +802 5 2 0 1 778 792 796 784 350 789 795 781 +803 5 2 0 1 234 209 226 235 792 366 793 796 +804 5 2 0 1 792 366 793 796 789 294 790 795 +805 5 2 0 1 229 235 232 215 784 796 794 779 +806 5 2 0 1 784 796 794 779 781 795 791 352 +807 5 2 0 1 235 226 211 232 796 793 788 794 +808 5 2 0 1 796 793 788 794 795 790 351 791 +809 5 2 0 1 33 228 233 34 154 785 800 156 +810 5 2 0 1 154 785 800 156 134 782 799 155 +811 5 2 0 1 228 215 232 233 785 779 794 800 +812 5 2 0 1 785 779 794 800 782 352 791 799 +813 5 2 0 1 34 233 220 2 156 800 798 140 +814 5 2 0 1 156 800 798 140 155 799 797 128 +815 5 2 0 1 233 232 211 220 800 794 788 798 +816 5 2 0 1 800 794 788 798 799 791 351 797 +817 5 2 0 1 127 408 801 160 143 802 803 161 +818 5 2 0 1 143 802 803 161 128 797 799 155 +819 5 2 0 1 408 294 789 801 802 790 795 803 +820 5 2 0 1 802 790 795 803 797 351 791 799 +821 5 2 0 1 160 801 780 131 161 803 786 158 +822 5 2 0 1 161 803 786 158 155 799 782 134 +823 5 2 0 1 801 789 350 780 803 795 781 786 +824 5 2 0 1 803 795 781 786 799 791 352 782 +825 5 2 0 1 6 22 98 95 157 164 813 809 +826 5 2 0 1 157 164 813 809 131 163 812 806 +827 5 2 0 1 22 21 97 98 164 162 810 813 +828 5 2 0 1 164 162 810 813 163 135 807 812 +829 5 2 0 1 95 98 96 87 809 813 811 805 +830 5 2 0 1 809 813 811 805 806 812 808 353 +831 5 2 0 1 98 97 91 96 813 810 804 811 +832 5 2 0 1 813 810 804 811 812 807 355 808 +833 5 2 0 1 21 23 104 97 162 166 820 810 +834 5 2 0 1 162 166 820 810 135 165 819 807 +835 5 2 0 1 23 5 103 104 166 148 817 820 +836 5 2 0 1 166 148 817 820 165 130 815 819 +837 5 2 0 1 97 104 101 91 810 820 818 804 +838 5 2 0 1 810 820 818 804 807 819 816 355 +839 5 2 0 1 104 103 88 101 820 817 814 818 +840 5 2 0 1 820 817 814 818 819 815 354 816 +841 5 2 0 1 87 96 102 99 805 811 826 823 +842 5 2 0 1 805 811 826 823 353 808 825 821 +843 5 2 0 1 96 91 101 102 811 804 818 826 +844 5 2 0 1 811 804 818 826 808 355 816 825 +845 5 2 0 1 99 102 100 86 823 826 824 385 +846 5 2 0 1 823 826 824 385 821 825 822 296 +847 5 2 0 1 102 101 88 100 826 818 814 824 +848 5 2 0 1 826 818 814 824 825 816 354 822 +849 5 2 0 1 127 152 167 160 407 827 829 828 +850 5 2 0 1 407 827 829 828 296 822 825 821 +851 5 2 0 1 152 130 165 167 827 815 819 829 +852 5 2 0 1 827 815 819 829 822 354 816 825 +853 5 2 0 1 160 167 163 131 828 829 812 806 +854 5 2 0 1 828 829 812 806 821 825 808 353 +855 5 2 0 1 167 165 135 163 829 819 807 812 +856 5 2 0 1 829 819 807 812 825 816 355 808 +857 5 2 0 1 45 433 617 61 399 506 836 745 +858 5 2 0 1 399 506 836 745 295 504 835 743 +859 5 2 0 1 433 297 599 617 506 491 833 836 +860 5 2 0 1 506 491 833 836 504 317 831 835 +861 5 2 0 1 61 617 615 46 745 836 834 727 +862 5 2 0 1 745 836 834 727 743 835 832 344 +863 5 2 0 1 617 599 330 615 836 833 830 834 +864 5 2 0 1 836 833 830 834 835 831 356 832 +865 5 2 0 1 297 414 604 599 491 498 840 833 +866 5 2 0 1 491 498 840 833 317 494 839 831 +867 5 2 0 1 414 250 263 604 498 490 590 840 +868 5 2 0 1 498 490 590 840 494 303 588 839 +869 5 2 0 1 599 604 600 330 833 840 838 830 +870 5 2 0 1 833 840 838 830 831 839 837 356 +871 5 2 0 1 604 263 251 600 840 590 580 838 +872 5 2 0 1 840 590 580 838 839 588 327 837 +873 5 2 0 1 46 615 620 58 727 834 842 731 +874 5 2 0 1 727 834 842 731 344 832 841 728 +875 5 2 0 1 615 330 600 620 834 830 838 842 +876 5 2 0 1 834 830 838 842 832 356 837 841 +877 5 2 0 1 58 620 259 1 731 842 583 138 +878 5 2 0 1 731 842 583 138 728 841 581 129 +879 5 2 0 1 620 600 251 259 842 838 580 583 +880 5 2 0 1 842 838 580 583 841 837 327 581 +881 5 2 0 1 127 514 515 409 144 593 843 750 +882 5 2 0 1 144 593 843 750 129 581 841 728 +883 5 2 0 1 514 303 494 515 593 588 839 843 +884 5 2 0 1 593 588 839 843 581 327 837 841 +885 5 2 0 1 409 515 504 295 750 843 835 743 +886 5 2 0 1 750 843 835 743 728 841 832 344 +887 5 2 0 1 515 494 317 504 843 839 831 835 +888 5 2 0 1 843 839 831 835 841 837 356 832 +889 5 2 0 1 3 67 568 177 267 611 850 637 +890 5 2 0 1 267 611 850 637 252 596 849 635 +891 5 2 0 1 67 48 548 568 611 609 847 850 +892 5 2 0 1 611 609 847 850 596 329 845 849 +893 5 2 0 1 177 568 558 169 637 850 848 622 +894 5 2 0 1 637 850 848 622 635 849 846 332 +895 5 2 0 1 568 548 324 558 850 847 844 848 +896 5 2 0 1 850 847 844 848 849 845 357 846 +897 5 2 0 1 48 70 552 548 609 616 854 847 +898 5 2 0 1 609 616 854 847 329 597 853 845 +899 5 2 0 1 70 45 434 552 616 433 440 854 +900 5 2 0 1 616 433 440 854 597 297 417 853 +901 5 2 0 1 548 552 545 324 847 854 852 844 +902 5 2 0 1 847 854 852 844 845 853 851 357 +903 5 2 0 1 552 434 298 545 854 440 435 852 +904 5 2 0 1 854 440 435 852 853 417 309 851 +905 5 2 0 1 169 558 561 181 622 848 856 629 +906 5 2 0 1 622 848 856 629 332 846 855 623 +907 5 2 0 1 558 324 545 561 848 844 852 856 +908 5 2 0 1 848 844 852 856 846 357 851 855 +909 5 2 0 1 181 561 447 168 629 856 451 446 +910 5 2 0 1 629 856 451 446 623 855 419 299 +911 5 2 0 1 561 545 298 447 856 852 435 451 +912 5 2 0 1 856 852 435 451 855 851 309 419 +913 5 2 0 1 250 414 603 264 415 426 857 642 +914 5 2 0 1 415 426 857 642 299 419 855 623 +915 5 2 0 1 414 297 597 603 426 417 853 857 +916 5 2 0 1 426 417 853 857 419 309 851 855 +917 5 2 0 1 264 603 596 252 642 857 849 635 +918 5 2 0 1 642 857 849 635 623 855 846 332 +919 5 2 0 1 603 597 329 596 857 853 845 849 +920 5 2 0 1 857 853 845 849 855 851 357 846 +921 5 2 0 1 250 415 643 280 416 428 864 693 +922 5 2 0 1 416 428 864 693 300 423 863 691 +923 5 2 0 1 415 299 625 643 428 421 861 864 +924 5 2 0 1 428 421 861 864 423 311 859 863 +925 5 2 0 1 280 643 641 254 693 864 862 675 +926 5 2 0 1 693 864 862 675 691 863 860 338 +927 5 2 0 1 643 625 333 641 864 861 858 862 +928 5 2 0 1 864 861 858 862 863 859 358 860 +929 5 2 0 1 299 446 630 625 421 452 868 861 +930 5 2 0 1 421 452 868 861 311 449 867 859 +931 5 2 0 1 446 168 191 630 452 448 525 868 +932 5 2 0 1 452 448 525 868 449 301 519 867 +933 5 2 0 1 625 630 626 333 861 868 866 858 +934 5 2 0 1 861 868 866 858 859 867 865 358 +935 5 2 0 1 630 191 171 626 868 525 518 866 +936 5 2 0 1 868 525 518 866 867 519 320 865 +937 5 2 0 1 254 641 646 283 675 862 870 679 +938 5 2 0 1 675 862 870 679 338 860 869 676 +939 5 2 0 1 641 333 626 646 862 858 866 870 +940 5 2 0 1 862 858 866 870 860 358 865 869 +941 5 2 0 1 283 646 187 7 679 870 533 111 +942 5 2 0 1 679 870 533 111 676 869 531 89 +943 5 2 0 1 646 626 171 187 870 866 518 533 +944 5 2 0 1 870 866 518 533 869 865 320 531 +945 5 2 0 1 86 457 460 456 108 538 871 698 +946 5 2 0 1 108 538 871 698 89 531 869 676 +947 5 2 0 1 457 301 449 460 538 519 867 871 +948 5 2 0 1 538 519 867 871 531 320 865 869 +949 5 2 0 1 456 460 423 300 698 871 863 691 +950 5 2 0 1 698 871 863 691 676 869 860 338 +951 5 2 0 1 460 449 311 423 871 867 859 863 +952 5 2 0 1 871 867 859 863 869 865 358 860 +953 5 2 0 1 4 185 565 74 222 767 878 720 +954 5 2 0 1 222 767 878 720 210 765 877 719 +955 5 2 0 1 185 170 557 565 767 752 875 878 +956 5 2 0 1 767 752 875 878 765 347 873 877 +957 5 2 0 1 74 565 546 49 720 878 876 702 +958 5 2 0 1 720 878 876 702 719 877 874 342 +959 5 2 0 1 565 557 323 546 878 875 872 876 +960 5 2 0 1 878 875 872 876 877 873 359 874 +961 5 2 0 1 170 182 560 557 752 759 882 875 +962 5 2 0 1 752 759 882 875 347 755 881 873 +963 5 2 0 1 182 168 447 560 759 479 483 882 +964 5 2 0 1 759 479 483 882 755 302 481 881 +965 5 2 0 1 557 560 544 323 875 882 880 872 +966 5 2 0 1 875 882 880 872 873 881 879 359 +967 5 2 0 1 560 447 298 544 882 483 472 880 +968 5 2 0 1 882 483 472 880 881 481 315 879 +969 5 2 0 1 49 546 551 78 702 876 884 708 +970 5 2 0 1 702 876 884 708 342 874 883 704 +971 5 2 0 1 546 323 544 551 876 872 880 884 +972 5 2 0 1 876 872 880 884 874 359 879 883 +973 5 2 0 1 78 551 434 45 708 884 475 397 +974 5 2 0 1 708 884 475 397 704 883 473 293 +975 5 2 0 1 551 544 298 434 884 880 472 475 +976 5 2 0 1 884 880 472 475 883 879 315 473 +977 5 2 0 1 209 486 775 225 365 487 885 724 +978 5 2 0 1 365 487 885 724 293 473 883 704 +979 5 2 0 1 486 302 755 775 487 481 881 885 +980 5 2 0 1 487 481 881 885 473 315 879 883 +981 5 2 0 1 225 775 765 210 724 885 877 719 +982 5 2 0 1 724 885 877 719 704 883 874 342 +983 5 2 0 1 775 755 347 765 885 881 873 877 +984 5 2 0 1 885 881 873 877 883 879 359 874 +985 5 2 0 1 8 113 542 201 239 663 892 772 +986 5 2 0 1 239 663 892 772 213 648 891 771 +987 5 2 0 1 113 90 537 542 663 661 889 892 +988 5 2 0 1 663 661 889 892 648 335 887 891 +989 5 2 0 1 201 542 522 172 772 892 890 754 +990 5 2 0 1 772 892 890 754 771 891 888 348 +991 5 2 0 1 542 537 321 522 892 889 886 890 +992 5 2 0 1 892 889 886 890 891 887 360 888 +993 5 2 0 1 90 117 539 537 661 668 896 889 +994 5 2 0 1 661 668 896 889 335 649 895 887 +995 5 2 0 1 117 86 457 539 668 383 467 896 +996 5 2 0 1 668 383 467 896 649 291 464 895 +997 5 2 0 1 537 539 521 321 889 896 894 886 +998 5 2 0 1 889 896 894 886 887 895 893 360 +999 5 2 0 1 539 457 301 521 896 467 463 894 +1000 5 2 0 1 896 467 463 894 895 464 314 893 +1001 5 2 0 1 172 522 526 198 754 890 898 760 +1002 5 2 0 1 754 890 898 760 348 888 897 756 +1003 5 2 0 1 522 321 521 526 890 886 894 898 +1004 5 2 0 1 890 886 894 898 888 360 893 897 +1005 5 2 0 1 198 526 448 168 760 898 482 479 +1006 5 2 0 1 760 898 482 479 756 897 480 302 +1007 5 2 0 1 526 521 301 448 898 894 463 482 +1008 5 2 0 1 898 894 463 482 897 893 314 480 +1009 5 2 0 1 209 364 655 242 486 488 899 776 +1010 5 2 0 1 486 488 899 776 302 480 897 756 +1011 5 2 0 1 364 291 649 655 488 464 895 899 +1012 5 2 0 1 488 464 895 899 480 314 893 897 +1013 5 2 0 1 242 655 648 213 776 899 891 771 +1014 5 2 0 1 776 899 891 771 756 897 888 348 +1015 5 2 0 1 655 649 335 648 899 895 887 891 +1016 5 2 0 1 899 895 887 891 897 893 360 888 +1017 5 2 0 1 86 383 669 99 385 392 906 823 +1018 5 2 0 1 385 392 906 823 296 388 905 821 +1019 5 2 0 1 383 291 651 669 392 368 903 906 +1020 5 2 0 1 392 368 903 906 388 305 901 905 +1021 5 2 0 1 99 669 667 87 823 906 904 805 +1022 5 2 0 1 823 906 904 805 821 905 902 353 +1023 5 2 0 1 669 651 336 667 906 903 900 904 +1024 5 2 0 1 906 903 900 904 905 901 361 902 +1025 5 2 0 1 291 364 656 651 368 377 910 903 +1026 5 2 0 1 368 377 910 903 305 372 909 901 +1027 5 2 0 1 364 209 234 656 377 366 792 910 +1028 5 2 0 1 377 366 792 910 372 294 789 909 +1029 5 2 0 1 651 656 652 336 903 910 908 900 +1030 5 2 0 1 903 910 908 900 901 909 907 361 +1031 5 2 0 1 656 234 212 652 910 792 778 908 +1032 5 2 0 1 910 792 778 908 909 789 350 907 +1033 5 2 0 1 87 667 672 95 805 904 912 809 +1034 5 2 0 1 805 904 912 809 353 902 911 806 +1035 5 2 0 1 667 336 652 672 904 900 908 912 +1036 5 2 0 1 904 900 908 912 902 361 907 911 +1037 5 2 0 1 95 672 230 6 809 912 783 157 +1038 5 2 0 1 809 912 783 157 806 911 780 131 +1039 5 2 0 1 672 652 212 230 912 908 778 783 +1040 5 2 0 1 912 908 778 783 911 907 350 780 +1041 5 2 0 1 127 408 410 407 160 801 913 828 +1042 5 2 0 1 160 801 913 828 131 780 911 806 +1043 5 2 0 1 408 294 372 410 801 789 909 913 +1044 5 2 0 1 801 789 909 913 780 350 907 911 +1045 5 2 0 1 407 410 388 296 828 913 905 821 +1046 5 2 0 1 828 913 905 821 806 911 902 353 +1047 5 2 0 1 410 372 305 388 913 909 901 905 +1048 5 2 0 1 913 909 901 905 911 907 361 902 +1049 5 2 0 1 86 100 697 456 385 824 920 511 +1050 5 2 0 1 385 824 920 511 296 822 919 510 +1051 5 2 0 1 100 88 685 697 824 814 917 920 +1052 5 2 0 1 824 814 917 920 822 354 915 919 +1053 5 2 0 1 456 697 692 300 511 920 918 493 +1054 5 2 0 1 511 920 918 493 510 919 916 318 +1055 5 2 0 1 697 685 339 692 920 917 914 918 +1056 5 2 0 1 920 917 914 918 919 915 362 916 +1057 5 2 0 1 88 103 687 685 814 817 924 917 +1058 5 2 0 1 814 817 924 917 354 815 923 915 +1059 5 2 0 1 103 5 269 687 817 148 575 924 +1060 5 2 0 1 817 148 575 924 815 130 572 923 +1061 5 2 0 1 685 687 684 339 917 924 922 914 +1062 5 2 0 1 917 924 922 914 915 923 921 362 +1063 5 2 0 1 687 269 253 684 924 575 571 922 +1064 5 2 0 1 924 575 571 922 923 572 326 921 +1065 5 2 0 1 300 692 694 416 493 918 926 499 +1066 5 2 0 1 493 918 926 499 318 916 925 495 +1067 5 2 0 1 692 339 684 694 918 914 922 926 +1068 5 2 0 1 918 914 922 926 916 362 921 925 +1069 5 2 0 1 416 694 273 250 499 926 589 490 +1070 5 2 0 1 499 926 589 490 495 925 587 303 +1071 5 2 0 1 694 684 253 273 926 922 571 589 +1072 5 2 0 1 926 922 571 589 925 921 326 587 +1073 5 2 0 1 127 152 827 407 514 594 927 516 +1074 5 2 0 1 514 594 927 516 303 587 925 495 +1075 5 2 0 1 152 130 815 827 594 572 923 927 +1076 5 2 0 1 594 572 923 927 587 326 921 925 +1077 5 2 0 1 407 827 822 296 516 927 919 510 +1078 5 2 0 1 516 927 919 510 495 925 916 318 +1079 5 2 0 1 827 815 354 822 927 923 915 919 +1080 5 2 0 1 927 923 915 919 925 921 362 916 +1081 5 2 0 1 45 62 707 397 399 746 934 403 +1082 5 2 0 1 399 746 934 403 295 744 933 400 +1083 5 2 0 1 62 47 700 707 746 736 931 934 +1084 5 2 0 1 746 736 931 934 744 345 929 933 +1085 5 2 0 1 397 707 703 293 403 934 932 371 +1086 5 2 0 1 403 934 932 371 400 933 930 306 +1087 5 2 0 1 707 700 341 703 934 931 928 932 +1088 5 2 0 1 934 931 928 932 933 929 363 930 +1089 5 2 0 1 47 56 715 700 736 739 938 931 +1090 5 2 0 1 736 739 938 931 345 737 937 929 +1091 5 2 0 1 56 2 220 715 739 140 798 938 +1092 5 2 0 1 739 140 798 938 737 128 797 937 +1093 5 2 0 1 700 715 713 341 931 938 936 928 +1094 5 2 0 1 931 938 936 928 929 937 935 363 +1095 5 2 0 1 715 220 211 713 938 798 788 936 +1096 5 2 0 1 938 798 788 936 937 797 351 935 +1097 5 2 0 1 293 703 723 365 371 932 940 378 +1098 5 2 0 1 371 932 940 378 306 930 939 373 +1099 5 2 0 1 703 341 713 723 932 928 936 940 +1100 5 2 0 1 932 928 936 940 930 363 935 939 +1101 5 2 0 1 365 723 226 209 378 940 793 366 +1102 5 2 0 1 378 940 793 366 373 939 790 294 +1103 5 2 0 1 723 713 211 226 940 936 788 793 +1104 5 2 0 1 940 936 788 793 939 935 351 790 +1105 5 2 0 1 127 143 749 409 408 802 941 412 +1106 5 2 0 1 408 802 941 412 294 790 939 373 +1107 5 2 0 1 143 128 737 749 802 797 937 941 +1108 5 2 0 1 802 797 937 941 790 351 935 939 +1109 5 2 0 1 409 749 744 295 412 941 933 400 +1110 5 2 0 1 412 941 933 400 373 939 930 306 +1111 5 2 0 1 749 737 345 744 941 937 929 933 +1112 5 2 0 1 941 937 929 933 939 935 363 930 +$EndElements diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index cf3ac2021..489401d2b 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -28,14 +28,12 @@ from grudge.array_context import ( PytestNumpyArrayContextFactory, - PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory, ) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory, - PytestPytatoPyOpenCLArrayContextFactory, + [PytestPytatoPyOpenCLArrayContextFactory, PytestNumpyArrayContextFactory]) import logging @@ -150,8 +148,8 @@ def rhs(x): assert np.allclose(np.diag(mat), 3*2*2 + 2) -@pytest.mark.parametrize("dim", [1, 2]) -@pytest.mark.parametrize("degree", [2, 4]) +@pytest.mark.parametrize("dim", [1]) +@pytest.mark.parametrize("degree", [2]) def test_wave_dt_estimate(actx_factory, dim, degree, visualize=False): actx = actx_factory() diff --git a/test/test_euler_model.py b/test/test_euler_model.py index 6e9e9c1d2..42ec3d48d 100644 --- a/test/test_euler_model.py +++ b/test/test_euler_model.py @@ -31,12 +31,12 @@ ) from grudge import op -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import PytestNumpyArrayContextFactory logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestNumpyArrayContextFactory]) @pytest.mark.parametrize("order", [1, 2, 3]) diff --git a/test/test_grudge.py b/test/test_grudge.py index b0849310b..3673a2991 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -42,13 +42,13 @@ from pytools.obj_array import flat_obj_array from grudge import dof_desc, geometry, op -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import PytestNumpyArrayContextFactory from grudge.discretization import make_discretization_collection logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestNumpyArrayContextFactory]) # {{{ mass operator trig integration @@ -231,6 +231,7 @@ def test_mass_surface_area(actx_factory, name): "warped_rect2", "warped_rect3", "gh-339-1", + "gh-339-1-overint", "gh-339-4", ]) def test_mass_operator_inverse(actx_factory, name): @@ -249,6 +250,7 @@ def test_mass_operator_inverse(actx_factory, name): builder = mesh_data.SpheroidMeshBuilder() elif name.startswith("warped_rect"): builder = mesh_data.WarpedRectMeshBuilder(dim=int(name[-1])) + overintegrate = True elif name == "gh-339-1": builder = mesh_data.GmshMeshBuilder3D("gh-339.msh") order = 1 @@ -273,14 +275,25 @@ def test_mass_operator_inverse(actx_factory, name): eoc = EOCRecorder() - for resolution in builder.resolutions: - mesh = builder.get_mesh(resolution) + if name.startswith("gh-339"): + resolutions = ["gh-339.msh"] + for i in [1]: + resolutions.append(f"gh-339-refined-{i}.msh") + else: + resolutions = builder.resolutions + + for resolution in resolutions: + if name.startswith("gh-339"): + builder = mesh_data.GmshMeshBuilder3D(resolution) + mesh = builder.get_mesh(None) + else: + mesh = builder.get_mesh(resolution) dcoll = make_discretization_collection( actx, mesh, discr_tag_to_group_factory={ dof_desc.DISCR_TAG_BASE: ( InterpolatoryEdgeClusteredGroupFactory(order)), dof_desc.DISCR_TAG_QUAD: ( - QuadratureGroupFactory(order)) + QuadratureGroupFactory(2*order + 1)) }) volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL) @@ -292,27 +305,41 @@ def test_mass_operator_inverse(actx_factory, name): def f(x): return actx.np.cos(4.0 * x[0]) - x_volm = actx.thaw(volume_discr.nodes()) - f_volm = f(x_volm) - if not overintegrate: - dd = dof_desc.DD_VOLUME_ALL - f_inv = op.inverse_mass( - dcoll, op.mass(dcoll, dd, f_volm) - ) - else: - dd_base_vol = dof_desc.as_dofdesc( - dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_BASE) - dd_quad_vol = dof_desc.as_dofdesc( - dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_QUAD) - f_inv = op.inverse_mass( - dcoll, op.mass(dcoll, dd_quad_vol, - op.project(dcoll, dd_base_vol, dd_quad_vol, f_volm))) - - inv_error = actx.to_numpy( + def f_inv(dcoll=dcoll, + volume_discr=volume_discr, + use_tensor_product_fast_eval=True): + x_volm = actx.thaw(volume_discr.nodes()) + f_volm = f(x_volm) + if not overintegrate: + dd = dof_desc.DD_VOLUME_ALL + f_inv = op.inverse_mass( + dcoll, op.mass(dcoll, dd, f_volm) + ) + else: + dd_base_vol = dof_desc.as_dofdesc( + dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_BASE) + dd_quad_vol = dof_desc.as_dofdesc( + dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_QUAD) + f_inv = op.inverse_mass( + dcoll, dd_quad_vol, op.mass(dcoll, dd_quad_vol, + op.project(dcoll, dd_base_vol, dd_quad_vol, f_volm), + use_tensor_product_fast_eval=use_tensor_product_fast_eval), + use_tensor_product_fast_eval=use_tensor_product_fast_eval) + + return actx.to_numpy( op.norm(dcoll, f_volm - f_inv, 2) / op.norm(dcoll, f_volm, 2)) # }}} + if isinstance(mesh.groups[0], TensorProductElementGroup): + inv_error_fast = f_inv(use_tensor_product_fast_eval=True) + inv_error_slow = f_inv(use_tensor_product_fast_eval=False) + + assert np.abs(inv_error_fast - inv_error_slow) < 1e-14 + inv_error = inv_error_fast + else: + inv_error = f_inv() + # compute max element size from grudge.dt_utils import h_max_from_volume @@ -322,9 +349,9 @@ def f(x): logger.info("inverse mass error\n%s", str(eoc)) - # NOTE: both cases give 1.0e-16-ish at the moment, but just to be on the - # safe side, choose a slightly larger tolerance - assert eoc.max_error() < 1.0e-14 + print("L^inf error:") + print(eoc) + assert eoc.max_error() <= 1.0e-14 or eoc.order_estimate() >= order - 0.5 # }}} diff --git a/test/test_metrics.py b/test/test_metrics.py index 1ee043b8a..a119d4e1d 100644 --- a/test/test_metrics.py +++ b/test/test_metrics.py @@ -34,7 +34,6 @@ from grudge.array_context import ( PytestNumpyArrayContextFactory, - PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory, ) from grudge.discretization import make_discretization_collection @@ -42,8 +41,7 @@ logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory, - PytestPytatoPyOpenCLArrayContextFactory, + [PytestPytatoPyOpenCLArrayContextFactory, PytestNumpyArrayContextFactory]) diff --git a/test/test_modal_connections.py b/test/test_modal_connections.py index d6bebdbd5..99261c30d 100644 --- a/test/test_modal_connections.py +++ b/test/test_modal_connections.py @@ -23,12 +23,12 @@ from arraycontext import pytest_generate_tests_for_array_contexts -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import PytestNumpyArrayContextFactory from grudge.discretization import make_discretization_collection pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestNumpyArrayContextFactory]) import pytest diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py index 342c52676..18e88ec93 100644 --- a/test/test_mpi_communication.py +++ b/test/test_mpi_communication.py @@ -35,7 +35,11 @@ from pytools.obj_array import flat_obj_array from grudge import dof_desc, op -from grudge.array_context import MPIPyOpenCLArrayContext, MPIPytatoArrayContext +from grudge.array_context import ( + MPINumpyArrayContext, + MPIPyOpenCLArrayContext, + MPIPytatoArrayContext, +) from grudge.discretization import make_discretization_collection from grudge.shortcuts import rk4_step @@ -49,7 +53,7 @@ class SimpleTag: # {{{ mpi test infrastructure -DISTRIBUTED_ACTXS = [MPIPyOpenCLArrayContext, MPIPytatoArrayContext] +DISTRIBUTED_ACTXS = [MPINumpyArrayContext, MPIPytatoArrayContext] def run_test_with_mpi(num_ranks, f, *args): @@ -87,6 +91,8 @@ def run_test_with_mpi_inner(): actx = actx_class(comm, queue, mpi_base_tag=15000) elif actx_class is MPIPyOpenCLArrayContext: actx = actx_class(comm, queue, force_device_scalars=True) + elif actx_class is MPINumpyArrayContext: + actx = actx_class(comm) else: raise ValueError("unknown actx_class") diff --git a/test/test_op.py b/test/test_op.py index 17f49a074..6126e4dfb 100644 --- a/test/test_op.py +++ b/test/test_op.py @@ -32,11 +32,18 @@ InterpolatoryEdgeClusteredGroupFactory, QuadratureGroupFactory, ) -from meshmode.mesh import BTAG_ALL +from meshmode.mesh import ( + BTAG_ALL, + SimplexElementGroup, + TensorProductElementGroup, +) from pytools.obj_array import make_obj_array from grudge import geometry, op -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import ( + PytestNumpyArrayContextFactory, + PytestPytatoPyOpenCLArrayContextFactory, +) from grudge.discretization import make_discretization_collection from grudge.dof_desc import ( DISCR_TAG_BASE, @@ -52,7 +59,8 @@ logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestPytatoPyOpenCLArrayContextFactory] +) # {{{ gradient @@ -60,14 +68,18 @@ @pytest.mark.parametrize("form", ["strong", "weak", "weak-overint"]) @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("order", [2, 3]) -@pytest.mark.parametrize("warp_mesh", [False, True]) +@pytest.mark.parametrize("warp_mesh", [True, False]) @pytest.mark.parametrize(("vectorize", "nested"), [ (False, False), (True, False), (True, True) ]) +@pytest.mark.parametrize("group_cls", [ + SimplexElementGroup, + TensorProductElementGroup +]) def test_gradient(actx_factory, form, dim, order, vectorize, nested, - warp_mesh, visualize=False): + warp_mesh, group_cls, visualize=False): actx = actx_factory() from pytools.convergence import EOCRecorder @@ -77,19 +89,27 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested, if warp_mesh: if dim == 1: pytest.skip("warped mesh in 1D not implemented") + + # NOTE: strong form fails for meshes with order > 1 + if group_cls == TensorProductElementGroup and form == "strong": + pytest.skip("strong form + mesh with order > 1 not implemented") + mesh = mgen.generate_warped_rect_mesh( - dim=dim, order=order, nelements_side=n) + dim=dim, order=order, nelements_side=n, + group_cls=group_cls) else: mesh = mgen.generate_regular_rect_mesh( a=(-1,)*dim, b=(1,)*dim, - nelements_per_axis=(n,)*dim) + nelements_per_axis=(n,)*dim, + group_cls=group_cls) dcoll = make_discretization_collection( - actx, mesh, - discr_tag_to_group_factory={ - DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order), - DISCR_TAG_QUAD: QuadratureGroupFactory(3 * order) - }) + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: + InterpolatoryEdgeClusteredGroupFactory(order), + DISCR_TAG_QUAD: QuadratureGroupFactory(3*order) + }) def f(x): result = 1 @@ -141,42 +161,47 @@ def get_flux(u_tpair, dcoll=dcoll): bdry_x = actx.thaw(dcoll.nodes(bdry_dd_base)) bdry_u = vectorize_if_requested(f(bdry_x)) - if form == "strong": - grad_u = ( - op.local_grad(dcoll, u, nested=nested) - # No flux terms because u doesn't have inter-el jumps + def compute_grad_u(): + if form == "strong": + grad_u = ( + op.local_grad(dcoll, u, nested=nested) + # No flux terms because u doesn't have inter-el jumps + ) + elif form.startswith("weak"): + assert form in ["weak", "weak-overint"] + if "overint" in form: + quad_discr_tag = DISCR_TAG_QUAD + else: + quad_discr_tag = DISCR_TAG_BASE + + allfaces_dd_base = as_dofdesc(FACE_RESTR_ALL, quad_discr_tag) + vol_dd_base = as_dofdesc(DTAG_VOLUME_ALL) + vol_dd_quad = vol_dd_base.with_discr_tag(quad_discr_tag) + bdry_dd_quad = bdry_dd_base.with_discr_tag(quad_discr_tag) + allfaces_dd_quad = allfaces_dd_base.with_discr_tag(quad_discr_tag) + + grad_u = op.inverse_mass(dcoll, vol_dd_quad, + -op.weak_local_grad(dcoll, vol_dd_quad, + op.project(dcoll, vol_dd_base, vol_dd_quad, u), + nested=nested) + + + op.face_mass(dcoll, + allfaces_dd_quad, + sum(get_flux( + op.project_tracepair(dcoll, allfaces_dd_quad, utpair)) + for utpair in op.interior_trace_pairs( + dcoll, u, volume_dd=vol_dd_base)) + + get_flux( + op.project_tracepair(dcoll, bdry_dd_quad, + bv_trace_pair(dcoll, bdry_dd_base, u, bdry_u))) + ) ) - elif form.startswith("weak"): - assert form in ["weak", "weak-overint"] - if "overint" in form: - quad_discr_tag = DISCR_TAG_QUAD else: - quad_discr_tag = DISCR_TAG_BASE + raise ValueError("Invalid form argument.") - allfaces_dd_base = as_dofdesc(FACE_RESTR_ALL, quad_discr_tag) - vol_dd_base = as_dofdesc(DTAG_VOLUME_ALL) - vol_dd_quad = vol_dd_base.with_discr_tag(quad_discr_tag) - bdry_dd_quad = bdry_dd_base.with_discr_tag(quad_discr_tag) - allfaces_dd_quad = allfaces_dd_base.with_discr_tag(quad_discr_tag) + return grad_u - grad_u = op.inverse_mass(dcoll, - -op.weak_local_grad(dcoll, vol_dd_quad, - op.project(dcoll, vol_dd_base, vol_dd_quad, u), - nested=nested) - + - op.face_mass(dcoll, - allfaces_dd_quad, - sum(get_flux( - op.project_tracepair(dcoll, allfaces_dd_quad, utpair)) - for utpair in op.interior_trace_pairs( - dcoll, u, volume_dd=vol_dd_base)) - + get_flux( - op.project_tracepair(dcoll, bdry_dd_quad, - bv_trace_pair(dcoll, bdry_dd_base, u, bdry_u))) - ) - ) - else: - raise ValueError("Invalid form argument.") + grad_u = actx.thaw(actx.freeze(actx.compile(compute_grad_u)())) if vectorize: expected_grad_u = make_obj_array( @@ -186,6 +211,8 @@ def get_flux(u_tpair, dcoll=dcoll): else: expected_grad_u = grad_f(x) + expected_grad_u = actx.thaw(actx.freeze(expected_grad_u)) + if visualize: # the code below does not handle the vectorized case assert not vectorize @@ -198,12 +225,11 @@ def get_flux(u_tpair, dcoll=dcoll): vis.write_vtk_file(filename, [ ("u", u), ("grad_u", grad_u), - ("expected_grad_u", expected_grad_u), - ], overwrite=True) + ("expected_grad_u", expected_grad_u), ], overwrite=True) rel_linf_err = actx.to_numpy( - op.norm(dcoll, grad_u - expected_grad_u, np.inf) - / op.norm(dcoll, expected_grad_u, np.inf)) + op.norm(dcoll, grad_u - expected_grad_u, 2) + / op.norm(dcoll, expected_grad_u, 2)) eoc_rec.add_data_point(1./n, rel_linf_err) print("L^inf error:") @@ -216,7 +242,7 @@ def get_flux(u_tpair, dcoll=dcoll): # {{{ divergence -@pytest.mark.parametrize("form", ["strong", "weak"]) +@pytest.mark.parametrize("form", ["strong", "weak", "weak-overint"]) @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("order", [2, 3]) @pytest.mark.parametrize(("vectorize", "nested"), [ @@ -224,8 +250,13 @@ def get_flux(u_tpair, dcoll=dcoll): (True, False), (True, True) ]) +@pytest.mark.parametrize("group_cls", [ + SimplexElementGroup, + TensorProductElementGroup +]) def test_divergence(actx_factory, form, dim, order, vectorize, nested, - visualize=False): + group_cls, visualize=False): + actx = actx_factory() from pytools.convergence import EOCRecorder @@ -234,9 +265,16 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested, for n in [4, 6, 8]: mesh = mgen.generate_regular_rect_mesh( a=(-1,)*dim, b=(1,)*dim, - nelements_per_axis=(n,)*dim) + nelements_per_axis=(n,)*dim, + group_cls=group_cls) - dcoll = make_discretization_collection(actx, mesh, order=order) + dcoll = make_discretization_collection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: + InterpolatoryEdgeClusteredGroupFactory(order), + DISCR_TAG_QUAD: QuadratureGroupFactory(3*order) + }) def f(x, dcoll=dcoll): result = make_obj_array([dcoll.zeros(actx) + (i+1) for i in range(dim)]) @@ -263,14 +301,14 @@ def div_f(x, dcoll=dcoll): result = result + deriv return result - x = actx.thaw(dcoll.nodes()) + def vectorize_if_requested(vec): + if vectorize: + return make_obj_array([(i+1)*vec for i in range(dim)]) + else: + return vec - if vectorize: - u = make_obj_array([(i+1)*f(x) for i in range(dim)]) - if not nested: - u = np.stack(u, axis=0) - else: - u = f(x) + x = actx.thaw(dcoll.nodes()) + u = vectorize_if_requested(f(x)) def get_flux(u_tpair, dcoll=dcoll): dd = u_tpair.dd @@ -281,24 +319,34 @@ def get_flux(u_tpair, dcoll=dcoll): flux = u_tpair.avg @ normal return op.project(dcoll, dd, dd_allfaces, flux) - dd_allfaces = as_dofdesc(FACE_RESTR_ALL) - if form == "strong": div_u = ( op.local_div(dcoll, u) # No flux terms because u doesn't have inter-el jumps ) - elif form == "weak": - div_u = op.inverse_mass(dcoll, - -op.weak_local_div(dcoll, u) + elif form.startswith("weak"): + assert form in ["weak", "weak-overint"] + if "overint" in form: + quad_discr_tag = DISCR_TAG_QUAD + else: + quad_discr_tag = DISCR_TAG_BASE + + allfaces_dd_base = as_dofdesc(FACE_RESTR_ALL, quad_discr_tag) + vol_dd_base = as_dofdesc(DTAG_VOLUME_ALL) + vol_dd_quad = vol_dd_base.with_discr_tag(quad_discr_tag) + allfaces_dd_quad = allfaces_dd_base.with_discr_tag(quad_discr_tag) + + div_u = op.inverse_mass(dcoll, vol_dd_quad, + -op.weak_local_div(dcoll, vol_dd_quad, + op.project(dcoll, vol_dd_base, vol_dd_quad, u)) + op.face_mass(dcoll, - dd_allfaces, + allfaces_dd_quad, # Note: no boundary flux terms here because u_ext == u_int == 0 - sum(get_flux(utpair) - for utpair in op.interior_trace_pairs(dcoll, u)) - ) - ) + sum(get_flux( + op.project_tracepair(dcoll, allfaces_dd_quad, utpair)) + for utpair in op.interior_trace_pairs( + dcoll, u, volume_dd=vol_dd_base)))) else: raise ValueError("Invalid form argument.") diff --git a/test/test_reductions.py b/test/test_reductions.py index 44850a54e..4d390a3f9 100644 --- a/test/test_reductions.py +++ b/test/test_reductions.py @@ -39,13 +39,13 @@ from pytools.obj_array import make_obj_array from grudge import op -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import PytestNumpyArrayContextFactory from grudge.discretization import make_discretization_collection logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestNumpyArrayContextFactory]) @pytest.mark.parametrize(("mesh_size", "with_initial"), [ diff --git a/test/test_tools.py b/test/test_tools.py index efd5747d5..b4db8b7ee 100644 --- a/test/test_tools.py +++ b/test/test_tools.py @@ -29,11 +29,11 @@ from arraycontext import pytest_generate_tests_for_array_contexts -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import PytestNumpyArrayContextFactory pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestNumpyArrayContextFactory]) import logging diff --git a/test/test_trace_pair.py b/test/test_trace_pair.py index 5cb060c7e..e851ce2af 100644 --- a/test/test_trace_pair.py +++ b/test/test_trace_pair.py @@ -29,14 +29,14 @@ from arraycontext import pytest_generate_tests_for_array_contexts from meshmode.dof_array import DOFArray -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import PytestNumpyArrayContextFactory from grudge.discretization import make_discretization_collection from grudge.trace_pair import TracePair logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestNumpyArrayContextFactory]) def test_trace_pair(actx_factory):