diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py index d4b977b129de7ef7019d4b237023c73095d91406..a96923dbb792ed9874595a3d84177a4f4d429635 100755 --- a/examples/dagrt-fusion.py +++ b/examples/dagrt-fusion.py @@ -651,6 +651,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper): # TODO: Not comprehensive. sym_op.InterpolationOperator, sym_op.RefFaceMassOperator, + sym_op.RefMassOperator, sym_op.RefInverseMassOperator, sym_op.OppositeInteriorFaceSwap)): val = self.map_profiled_essentially_elementwise_linear( diff --git a/grudge/execution.py b/grudge/execution.py index a95b15e75e43917688f2c3766651e3018ff9d526..b3faae20d4cc6862f81f3d2c2071ad1e53fa8a8b 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -63,7 +63,7 @@ class ExecutionMapper(mappers.Evaluator, discr = self.discrwb.discr_from_dd(expr.dd) result = discr.empty(self.queue, allocator=self.bound_op.allocator) - result.fill(1) + result.fill(1.0) return result def map_node_coordinate_component(self, expr): @@ -626,7 +626,9 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, )(sym_operator) dumper("before-global-to-reference", sym_operator) - sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator) + sym_operator = mappers.GlobalToReferenceMapper( + discrwb.ambient_dim, + dim=discrwb.dim)(sym_operator) dumper("before-distributed", sym_operator) diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index 51eb5202acaa1d1ffa6c4cbc8b293eef4a25e54e..7b36ee9f73a7fdf0143619a06617789f7a8a7862 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -923,7 +923,7 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper): % expr.function) def map_ones(self, expr): - return 1 + return 1.0 def map_node_coordinate_component(self, expr): return self.map_variable_ref_expr( diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 36b57d50aaf333fbd3303fd72f1425f5de085086..77f57a95619d3879433096b304b351f4670c8160 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -587,15 +587,19 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): reference elements, together with explicit multiplication by geometric factors. """ - def __init__(self, ambient_dim, dim=None): + def __init__(self, ambient_dim, dim=None, use_wadg=None): CSECachingMapperMixin.__init__(self) IdentityMapper.__init__(self) if dim is None: dim = ambient_dim + if use_wadg is None: + use_wadg = True + self.ambient_dim = ambient_dim self.dim = dim + self.use_wadg = use_wadg map_common_subexpression_uncached = \ IdentityMapper.map_common_subexpression @@ -605,40 +609,56 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): # if we encounter non-quadrature operators here, we know they # must be nodal. - if expr.op.dd_in.is_volume(): + dd_in = expr.op.dd_in + dd_out = expr.op.dd_out + + if dd_in.is_volume(): dim = self.dim else: dim = self.dim - 1 - jac_in = sym.area_element(self.ambient_dim, dim, dd=expr.op.dd_in) + jac_in = sym.area_element(self.ambient_dim, dim, dd=dd_in) jac_noquad = sym.area_element(self.ambient_dim, dim, - dd=expr.op.dd_in.with_qtag(sym.QTAG_NONE)) + dd=dd_in.with_qtag(sym.QTAG_NONE)) def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True): - jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in) - rec_field = self.rec(field) if with_jacobian: + jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in) rec_field = jac_tag * rec_field + + def imd(rst): + return sym.inverse_metric_derivative( + rst, expr.op.xyz_axis, + ambient_dim=self.ambient_dim, dim=self.dim, + dd=dd_in) + return sum( - sym.inverse_metric_derivative( - rst_axis, expr.op.xyz_axis, - ambient_dim=self.ambient_dim, dim=self.dim) - * ref_class(rst_axis, dd_in=dd_in)(rec_field) + ref_class(rst_axis, dd_in=dd_in)(imd(rst_axis) * rec_field) for rst_axis in range(self.dim)) if isinstance(expr.op, op.MassOperator): - return op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)( + return op.RefMassOperator(dd_in, dd_out)( jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.InverseMassOperator): - return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)( - 1/jac_in * self.rec(expr.field)) + if self.use_wadg: + # TODO: this is also required for flat curvilinear elements + # based on https://arxiv.org/pdf/1608.03836.pdf + return op.RefInverseMassOperator(dd_in, dd_out)( + op.RefMassOperator(dd_in, dd_out)( + 1.0/jac_in * op.RefInverseMassOperator(dd_in, dd_out)( + self.rec(expr.field)) + ) + ) + else: + return op.RefInverseMassOperator(dd_in, dd_out)( + 1/jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.FaceMassOperator): - jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1, - dd=expr.op.dd_in) - return op.RefFaceMassOperator(expr.op.dd_in, expr.op.dd_out)( + jac_in_surf = sym.area_element( + self.ambient_dim, self.dim - 1, dd=dd_in) + return op.RefFaceMassOperator(dd_in, dd_out)( jac_in_surf * self.rec(expr.field)) elif isinstance(expr.op, op.StiffnessOperator): @@ -650,12 +670,12 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): elif isinstance(expr.op, op.DiffOperator): return rewrite_derivative( op.RefDiffOperator, - expr.field, dd_in=expr.op.dd_in, with_jacobian=False) + expr.field, dd_in=dd_in, with_jacobian=False) elif isinstance(expr.op, op.StiffnessTOperator): return rewrite_derivative( op.RefStiffnessTOperator, - expr.field, dd_in=expr.op.dd_in) + expr.field, dd_in=dd_in) elif isinstance(expr.op, op.MInvSTOperator): return self.rec( diff --git a/test/ellipsoid.geo b/test/ellipsoid.geo new file mode 100644 index 0000000000000000000000000000000000000000..5dd9d0ecaa025f2762f404ac1ab27e3249af25e7 --- /dev/null +++ b/test/ellipsoid.geo @@ -0,0 +1,4 @@ +SetFactory("OpenCASCADE"); +Sphere(1) = {0, 0, 0, 0.5}; +Dilate {{0, 0, 0}, {0.5, 0.5, 1.0}} { Volume{1}; } +Mesh.Algorithm = 6; diff --git a/test/mesh_data.py b/test/mesh_data.py new file mode 100644 index 0000000000000000000000000000000000000000..d53b104011ec269791d32bddab2492ab3c428a52 --- /dev/null +++ b/test/mesh_data.py @@ -0,0 +1,108 @@ +import six +import numpy as np + + +class MeshBuilder: + order = 4 + + def __init__(self, **kwargs): + for k, v in six.iteritems(kwargs): + setattr(self, k, v) + + @property + def mesh_order(self): + return self.order + + @property + def ambient_dim(self): + raise NotImplementedError + + @property + def resolutions(self): + raise NotImplementedError + + def get_mesh(self, resolution, mesh_order): + raise NotImplementedError + + +class EllipseMeshBuilder(MeshBuilder): + ambient_dim = 2 + + radius = 3.1 + aspect_ratio = 2.0 + resolutions = [16, 32, 64, 128] + + @property + def curve_fn(self): + from meshmode.mesh.generation import ellipse + return lambda t: self.radius * ellipse(self.aspect_ratio, t) + + def get_mesh(self, resolution, mesh_order): + from meshmode.mesh.generation import make_curve_mesh + return make_curve_mesh( + self.curve_fn, + np.linspace(0.0, 1.0, resolution + 1), + mesh_order) + + +class SphereMeshBuilder(MeshBuilder): + ambient_dim = 3 + + resolutions = [0, 1, 2, 3] + radius = 1.0 + + def get_mesh(self, resolution, mesh_order): + from meshmode.mesh.generation import generate_icosphere + return generate_icosphere(self.radius, order=mesh_order, + uniform_refinement_rounds=resolution) + + +class SpheroidMeshBuilder(MeshBuilder): + ambient_dim = 3 + + order = 2 + resolutions = [1.0, 0.1, 0.05] + + @property + def radius(self): + return 0.25 + + @property + def aspect_ratio(self): + return 2.0 + + def get_mesh(self, resolution, mesh_order): + import os + filename = os.path.join(os.path.dirname(__file__), "ellipsoid.geo") + + from meshmode.mesh.io import generate_gmsh, FileSource + mesh = generate_gmsh( + FileSource(filename), 2, order=mesh_order, + other_options=[ + "-string", + "Mesh.CharacteristicLengthMax = %g;" % resolution]) + + from meshmode.mesh.processing import perform_flips + return perform_flips(mesh, np.ones(mesh.nelements)) + + +class BoxMeshBuilder(MeshBuilder): + ambient_dim = 2 + + mesh_order = 1 + resolutions = [8, 16, 32] + + a = (-0.5, -0.5, -0.5) + b = (+0.5, +0.5, +0.5) + + def get_mesh(self, resolution, mesh_order): + if not isinstance(resolution, (list, tuple)): + resolution = (resolution,) * self.ambient_dim + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=self.a, b=self.b, + n=resolution, + order=mesh_order) + + return mesh diff --git a/test/test_grudge.py b/test/test_grudge.py index 74512768e61f2aba372ad9d002a3d2f4b1778731..81c55f041b6bf1c5ca1c065f6071124fbc1e132b 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -30,21 +30,19 @@ import pyopencl as cl import pyopencl.array import pyopencl.clmath +from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields -import pytest # noqa - +import pytest from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) import logging - -logger = logging.getLogger(__name__) logging.basicConfig() +logger = logging.getLogger(__name__) logger.setLevel(logging.INFO) -from grudge import sym, bind, DGDiscretizationWithBoundaries - @pytest.mark.parametrize("dim", [2, 3]) def test_inverse_metric(ctx_factory, dim): @@ -157,6 +155,118 @@ def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): assert err_3 < 5.0e-10, err_3 +def _ellipse_surface_area(radius, aspect_ratio): + # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html + eccentricity = 1.0 - (1/aspect_ratio)**2 + + if abs(aspect_ratio - 2.0) < 1.0e-14: + # NOTE: hardcoded value so we don't need scipy for the test + ellip_e = 1.2110560275684594 + else: + from scipy.special import ellipe + ellip_e = ellipe(eccentricity) + + return 4.0 * radius * ellip_e + + +def _spheroid_surface_area(radius, aspect_ratio): + # https://en.wikipedia.org/wiki/Ellipsoid#Surface_area + a = 1.0 + c = aspect_ratio + + if aspect_ratio > 1: + e = np.sqrt(1.0 - (a/c)**2) + return 2.0 * np.pi * radius**2 * (1.0 + (c/a) / e * np.arcsin(e)) + else: + e = np.sqrt(1.0 - (c/a)**2) + return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e)) + + +@pytest.mark.parametrize("name", [ + "2-1-ellipse", "spheroid", "box2d", "box3d" + ]) +def test_mass_surface(ctx_factory, name): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # {{{ cases + + if name == "2-1-ellipse": + from mesh_data import EllipseMeshBuilder + builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio) + elif name == "spheroid": + from mesh_data import SpheroidMeshBuilder + builder = SpheroidMeshBuilder() + surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) + elif name == "box2d": + from mesh_data import BoxMeshBuilder + builder = BoxMeshBuilder(ambient_dim=2) + surface_area = 1.0 + elif name == "box3d": + from mesh_data import BoxMeshBuilder + builder = BoxMeshBuilder(ambient_dim=3) + surface_area = 1.0 + else: + raise ValueError("unknown geometry name: %s" % name) + + # }}} + + # {{{ convergence + + from pytools.convergence import EOCRecorder + eoc_area = EOCRecorder() + eoc_inv = EOCRecorder() + + for resolution in builder.resolutions: + mesh = builder.get_mesh(resolution, builder.mesh_order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + + x = discr.discr_from_dd("vol").nodes().get(queue) + logger.info("nnodes: %d", volume_discr.nnodes) + logger.info("bbox: %s", + [(np.min(x[n]), np.max(x[n])) for n in range(x.shape[0])]) + + # {{{ compute surface area + + dd = sym.DD_VOLUME + sym_op = sym.NodalSum(dd)(sym.MassOperator(dd, dd)(sym.Ones(dd))) + approx_surface_area = bind(discr, sym_op)(queue) + + logger.info("surface: got {:.5e} / expected {:.5e}".format( + approx_surface_area, surface_area)) + area_error = abs(approx_surface_area - surface_area) / abs(surface_area) + + # }}} + + # {{{ compute inverse mass + + sym_op = sym.InverseMassOperator(dd, dd)( + sym.MassOperator(dd, dd)(sym.Ones(dd))) + approx_inverse = bind(discr, sym_op)(queue) + + logger.info("inverse: got {:.5e} / expected {:.5e}".format( + cl.array.max(approx_inverse).get(queue), 1.0)) + inv_error = la.norm(approx_inverse.get(queue) - 1.0) + + # }}} + + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc_area.add_data_point(h_max, area_error) + eoc_inv.add_data_point(h_max, inv_error) + + # }}} + + logger.info("surface area error\n%s", str(eoc_area)) + logger.info("inverse mass error\n%s", str(eoc_inv)) + + assert eoc_area.max_error() < 1.0e-14 \ + or eoc_area.order_estimate() >= (builder.order - 1.0) + assert eoc_inv.max_error() < 5.0e-09 \ + or eoc_inv.order_estimate() >= (builder.order - 1.0) + + @pytest.mark.parametrize("dim", [1, 2, 3]) def test_tri_diff_mat(ctx_factory, dim, order=4): """Check differentiation matrix along the coordinate axes on a disk @@ -244,6 +354,7 @@ def test_2d_gauss_theorem(ctx_factory): ("disk", [0.1, 0.05]), ("rect2", [4, 8]), ("rect3", [4, 6]), + ("warped", [4, 6, 8]) ]) @pytest.mark.parametrize("op_type", ["strong", "weak"]) @pytest.mark.parametrize("flux_type", ["central"]) @@ -256,6 +367,9 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + if mesh_name == "warped" and op_type == "strong": + pytest.skip("expected failure for strong form on curvilinear grids") + from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() @@ -297,7 +411,12 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type dt_factor = 2 else: raise ValueError("dt_factor not known for %dd" % dim) + elif mesh_name == "warped": + dim = 2 + from meshmode.mesh.generation import generate_warped_rect_mesh + mesh = generate_warped_rect_mesh(dim, order, mesh_par) + dt_factor = 1 if dim == 2 else 2 else: raise ValueError("invalid mesh name: " + mesh_name) @@ -372,7 +491,11 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type abscissa_label="h", error_label="L2 Error")) - assert eoc_rec.order_estimate() > order + if mesh_name == "warped": + # NOTE: warped grids are hard.. be more forgiving + assert eoc_rec.order_estimate() > order - 0.25 + else: + assert eoc_rec.order_estimate() > order @pytest.mark.parametrize("order", [3, 4, 5]) @@ -626,7 +749,6 @@ if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker