From 65eb51869728e9859c97581f95eb10f9d11b589d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 21:26:39 -0500 Subject: [PATCH 1/4] use h_max in advection convergence test --- test/test_grudge.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index ca9089b0..e24d1ce9 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -376,7 +376,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type if mesh_name == "segment": from meshmode.mesh.generation import generate_box_mesh mesh = generate_box_mesh( - [np.linspace(-1.0, 1.0, mesh_par)], + [np.linspa 2f245d67 - tests: skip warped strong form testce(-1.0, 1.0, mesh_par)], order=order) dim = 1 -- GitLab From 31d9929a0c9bb81d71d5a371f824ad3873f1dffd Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 29 Apr 2020 10:19:03 -0500 Subject: [PATCH 2/4] make mass operator work on surfaces --- test/test_grudge.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test/test_grudge.py b/test/test_grudge.py index e24d1ce9..2981668c 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -32,6 +32,7 @@ import pyopencl.clmath from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields +from mesh_data import MeshBuilderFactory import pytest from pyopencl.tools import ( # noqa @@ -266,6 +267,42 @@ def test_mass_surface(ctx_factory, name): or eoc_inv.order_estimate() >= (builder.order - 1.0) +@pytest.mark.parametrize(("mesh_builder_factory", "surface_area"), [ + (MeshBuilderFactory("Ellipse", radius=1.25), 2.0 * np.pi * 1.25), + (MeshBuilderFactory("Sphere", radius=1.25), 4.0 * np.pi * 1.25**2), + (MeshBuilderFactory("Box", ambient_dim=2), 1.0), + (MeshBuilderFactory("Box", ambient_dim=3), 1.0), + ]) +def test_mass_surface(ctx_factory, mesh_builder_factory, surface_area): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from pytools.convergence import EOCRecorder + eoc = EOCRecorder() + builder = mesh_builder_factory() + + for resolution in builder.resolutions: + mesh = builder.get_mesh(resolution, builder.mesh_order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) + + # 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("got {:.5e} / expected {:.5e}".format( + approx_surface_area, surface_area)) + area_error = abs(approx_surface_area - surface_area) / abs(surface_area) + + # compute max element size + h_max = 1.0 / (resolution + 1) + eoc.add_data_point(h_max, area_error) + + logger.info("\n%s", str(eoc)) + assert eoc.max_error() < 1.0e-14 \ + or eoc.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 -- GitLab From dada197456da81f6f54374a98993b6b5e1e8c851 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 20:33:46 -0500 Subject: [PATCH 3/4] use geometries with non-constant jacobians --- test/test_grudge.py | 92 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 77 insertions(+), 15 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 2981668c..f9631f2e 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -32,7 +32,6 @@ import pyopencl.clmath from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields -from mesh_data import MeshBuilderFactory import pytest from pyopencl.tools import ( # noqa @@ -267,40 +266,103 @@ def test_mass_surface(ctx_factory, name): or eoc_inv.order_estimate() >= (builder.order - 1.0) -@pytest.mark.parametrize(("mesh_builder_factory", "surface_area"), [ - (MeshBuilderFactory("Ellipse", radius=1.25), 2.0 * np.pi * 1.25), - (MeshBuilderFactory("Sphere", radius=1.25), 4.0 * np.pi * 1.25**2), - (MeshBuilderFactory("Box", ambient_dim=2), 1.0), - (MeshBuilderFactory("Box", ambient_dim=3), 1.0), +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, mesh_builder_factory, surface_area): +def test_mass_surface(ctx_factory, name): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) + import mesh_data as data + if name == "2-1-ellipse": + builder = data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio) + elif name == "spheroid": + builder = data.SpheroidMeshBuilder() + surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) + elif name == "box2d": + builder = data.BoxMeshBuilder(ambient_dim=2) + surface_area = 1.0 + elif name == "box3d": + builder = data.BoxMeshBuilder(ambient_dim=3) + surface_area = 1.0 + else: + raise ValueError("unknown geometry name: %s" % name) + from pytools.convergence import EOCRecorder - eoc = EOCRecorder() - builder = mesh_builder_factory() + 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("got {:.5e} / expected {:.5e}".format( + 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) + # compute max element size - h_max = 1.0 / (resolution + 1) - eoc.add_data_point(h_max, area_error) + if resolution > 1: + h_max = 1.0 / (resolution + 1) + else: + h_max = resolution + 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)) - logger.info("\n%s", str(eoc)) - assert eoc.max_error() < 1.0e-14 \ - or eoc.order_estimate() >= (builder.order - 1.0) + 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]) -- GitLab From 92f36bbb56f5144924bc5e71d36cd347d40d5d18 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 29 Apr 2020 14:56:36 -0500 Subject: [PATCH 4/4] make face operator and normals work on surfaces --- grudge/symbolic/mappers/__init__.py | 3 + grudge/symbolic/primitives.py | 60 +++++++----- test/test_grudge.py | 140 +++++++++++++++------------- 3 files changed, 114 insertions(+), 89 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index bab27e72..b6fcc0fe 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -870,6 +870,9 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): def map_ones(self, expr, enclosing_prec): return "Ones:" + self._format_dd(expr.dd) + def map_signed_face_ones(self, expr, enclosing_prec): + return "SignedOnes:" + self._format_dd(expr.dd) + # {{{ geometry data def map_node_coordinate_component(self, expr, enclosing_prec): diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index aa00e8e6..d7584539 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -526,7 +526,9 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None): dim = ambient_dim if dim == 0: - return MultiVector(np.array([_SignedFaceOnes(dd)])) + from pymbolic.geometric_algebra import get_euclidean_space + return MultiVector(_SignedFaceOnes(dd), + space=get_euclidean_space(ambient_dim)) from pytools import product return product( @@ -614,14 +616,8 @@ def area_element(ambient_dim, dim=None, dd=None): "area_element", cse_scope.DISCRETIZATION) -def mv_normal(dd, ambient_dim, dim=None): - """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`.""" - +def surface_normal(ambient_dim, dim=None, dd=None): dd = as_dofdesc(dd) - - if not dd.is_trace(): - raise ValueError("may only request normals on boundaries") - if dim is None: dim = ambient_dim - 1 @@ -632,23 +628,43 @@ def mv_normal(dd, ambient_dim, dim=None): / area_element(ambient_dim, dim, dd=dd) # Dorst Section 3.7.2 - mv = pder << pder.I.inv() + return cse(pder << pder.I.inv(), + "surface_normal", + cse_scope.DISCRETIZATION) - if dim == 0: - # NOTE: when the mesh is 0D, we do not have a clear notion of - # `exterior normal`, so we just take the tangent of the 1D element, - # interpolate it to the element faces and make it signed - tangent = parametrization_derivative( - ambient_dim, dim=1, dd=DD_VOLUME) - tangent = tangent / sqrt(tangent.norm_squared()) - from grudge.symbolic.operators import interp - interp = interp(DD_VOLUME, dd) - mv = MultiVector(np.array([ - mv.as_scalar() * interp(t) for t in tangent.as_vector() - ])) +def mv_normal(dd, ambient_dim, dim=None): + """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`.""" + dd = as_dofdesc(dd) + if not dd.is_trace(): + raise ValueError("may only request normals on boundaries") + + if dim is None: + dim = ambient_dim - 1 + + if dim == ambient_dim - 1: + return surface_normal(ambient_dim, dim=dim, dd=dd) - return cse(mv, "normal", cse_scope.DISCRETIZATION) + # NOTE: In the case of (d - 2)-dimensional curves, we don't really have + # enough information on the face to decide what an "exterior face normal" + # is (e.g the "normal" to a 1D curve in 3D space is actually a + # "normal plane") + # + # The trick done here is that we take the pseudoscalar in the + # volume, move it to the curve and then contract the face tangent with + # it to get the correct face normal vector. + assert dim == ambient_dim - 2 + + from grudge.symbolic.operators import interp + pder_volm = ( + pseudoscalar(ambient_dim, dim + 1, dd=DD_VOLUME) + .map(interp(DD_VOLUME, dd))) + pder_face = pseudoscalar(ambient_dim, dim, dd=dd) + mv = cse(pder_face << pder_volm, "face_normal", cse_scope.DISCRETIZATION) + + return cse(mv / sqrt(mv.norm_squared()), + "unit_face_normal", + cse_scope.DISCRETIZATION) def normal(dd, ambient_dim, dim=None, quadrature_tag=None): diff --git a/test/test_grudge.py b/test/test_grudge.py index f9631f2e..caf7f3f2 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -266,106 +266,112 @@ def test_mass_surface(ctx_factory, name): or eoc_inv.order_estimate() >= (builder.order - 1.0) -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 +@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) +def test_face_normal_surface(ctx_factory, name): + """Check that face normals are orthogonal to the surface normal""" - 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 + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + # {{{ cases + + import mesh_data as data + if name == "2-1-ellipse": + builder = data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + elif name == "spheroid": + builder = data.SpheroidMeshBuilder() else: - from scipy.special import ellipe - ellip_e = ellipe(eccentricity) + raise ValueError("unknown geometry name: %s" % name) - 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 + mesh = builder.get_mesh(builder.resolutions[-1], builder.mesh_order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) - 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)) + dv = sym.DD_VOLUME + df = sym.as_dofdesc(sym.FACE_RESTR_ALL) + surface_normal_sym = sym.surface_normal( + mesh.ambient_dim, dim=mesh.dim, dd=dv).map( + sym.interp(dv, df)).as_vector() + surface_normal = bind(discr, surface_normal_sym)(queue) + + face_normal_sym = sym.normal(df, mesh.ambient_dim, dim=mesh.dim - 1) + face_normal = bind(discr, face_normal_sym)(queue) + + n_dot_n = cl.array.sum( + sum(nv * nf for nv, nf in zip(surface_normal, face_normal)) + ).get(queue) + + t_dot_n = 0 + if mesh.ambient_dim == 3: + # NOTE: there's only one face tangent in 3d + face_tangent_sym = sym.forward_metric_derivative_vector( + mesh.ambient_dim, 0, dd=df) + face_tangent = bind(discr, face_tangent_sym)(queue) + + t_dot_n = cl.array.sum( + sum(tf * nf for tf, nf in zip(face_tangent, face_normal)) + ).get(queue) + + # }}} + + logger.info("dot: %.5e %.5e", n_dot_n, t_dot_n) + assert abs(n_dot_n) < 1.0e-12 + assert abs(t_dot_n) < 1.0e-12 + + +@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) +def test_face_mass_surface(ctx_factory, name): + """Check that integrating the exterior normals over all faces gives 0.""" -@pytest.mark.parametrize("name", [ - "2-1-ellipse", "spheroid", "box2d", "box3d" - ]) -def test_mass_surface(ctx_factory, name): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) + # {{{ cases + import mesh_data as data if name == "2-1-ellipse": builder = data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) - surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio) elif name == "spheroid": builder = data.SpheroidMeshBuilder() - surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) - elif name == "box2d": - builder = data.BoxMeshBuilder(ambient_dim=2) - surface_area = 1.0 - elif name == "box3d": - builder = data.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() + eoc = 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) + df = sym.as_dofdesc(sym.FACE_RESTR_ALL) + sym_op = sym.NodalSum(dd)(sym.FaceMassOperator(df, dd)( + sym.normal(df, mesh.ambient_dim, dim=mesh.dim - 1))) - 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) + value = la.norm(bind(discr, sym_op)(queue)) + logger.info("integral %.5e", value) # compute max element size - if resolution > 1: - h_max = 1.0 / (resolution + 1) - else: - h_max = resolution - eoc_area.add_data_point(h_max, area_error) - eoc_inv.add_data_point(h_max, inv_error) + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc.add_data_point(h_max, value) - 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) + logger.info("\n%s", str(eoc)) + assert eoc.max_error() < 1.0e-12 \ + or eoc.order_estimate() >= (builder.order - 1.0) -@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize("dim", [1, 2]) def test_tri_diff_mat(ctx_factory, dim, order=4): """Check differentiation matrix along the coordinate axes on a disk @@ -475,7 +481,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type if mesh_name == "segment": from meshmode.mesh.generation import generate_box_mesh mesh = generate_box_mesh( - [np.linspa 2f245d67 - tests: skip warped strong form testce(-1.0, 1.0, mesh_par)], + [np.linspace(-1.0, 1.0, mesh_par)], order=order) dim = 1 -- GitLab