From d167917f839ce951f94efd840b5d77b5ecd69580 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 21:26:39 -0500 Subject: [PATCH 1/8] 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 53bd4975..4251df40 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -352,7 +352,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type for event in dt_stepper.run(t_end=final_time): if isinstance(event, dt_stepper.StateComputed): step += 1 - print(event.t) + # print(event.t) last_t = event.t last_u = event.state_component -- GitLab From 82faafe584194cc1f0e9cae8bdd869d7c9ff6f01 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 29 Apr 2020 10:19:03 -0500 Subject: [PATCH 2/8] make mass operator work on surfaces --- grudge/execution.py | 6 ++- grudge/symbolic/compiler.py | 2 +- test/mesh_data.py | 91 +++++++++++++++++++++++++++++++++++++ test/test_grudge.py | 50 +++++++++++++++++--- 4 files changed, 139 insertions(+), 10 deletions(-) create mode 100644 test/mesh_data.py diff --git a/grudge/execution.py b/grudge/execution.py index a95b15e7..b3faae20 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 51eb5202..7b36ee9f 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/test/mesh_data.py b/test/mesh_data.py new file mode 100644 index 00000000..42c1eb74 --- /dev/null +++ b/test/mesh_data.py @@ -0,0 +1,91 @@ +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 MeshBuilderFactory: + def __init__(self, mesh_name, **kwargs): + self.mesh_name = "{}MeshBuilder".format(mesh_name) + self.mesh_kwargs = kwargs + + def __call__(self, **kwargs): + mesh_kwargs = self.mesh_kwargs.copy() + mesh_kwargs.update(kwargs) + + return globals()[self.mesh_name](**mesh_kwargs) + + +class EllipseMeshBuilder(MeshBuilder): + ambient_dim = 2 + + radius = 3.1 + aspect_ratio = 1.0 + resolutions = [8, 16, 32, 64] + + @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 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 4251df40..a7550d69 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -30,17 +30,18 @@ import pyopencl as cl import pyopencl.array import pyopencl.clmath +from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields +from mesh_data import MeshBuilderFactory -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__) - -from grudge import sym, bind, DGDiscretizationWithBoundaries +logging.basicConfig(level=logging.INFO) @pytest.mark.parametrize("dim", [2, 3]) @@ -154,6 +155,42 @@ def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): assert err_3 < 5.0e-10, err_3 +@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 @@ -621,7 +658,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 -- GitLab From 7666189118f2385fcedc17ca8e5246ed4558dd8c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 20:33:46 -0500 Subject: [PATCH 3/8] use geometries with non-constant jacobians --- test/ellipsoid.geo | 4 ++ test/mesh_data.py | 45 +++++++++++++++------- test/test_grudge.py | 92 +++++++++++++++++++++++++++++++++++++-------- 3 files changed, 112 insertions(+), 29 deletions(-) create mode 100644 test/ellipsoid.geo diff --git a/test/ellipsoid.geo b/test/ellipsoid.geo new file mode 100644 index 00000000..5dd9d0ec --- /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 index 42c1eb74..d53b1040 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -25,24 +25,12 @@ class MeshBuilder: raise NotImplementedError -class MeshBuilderFactory: - def __init__(self, mesh_name, **kwargs): - self.mesh_name = "{}MeshBuilder".format(mesh_name) - self.mesh_kwargs = kwargs - - def __call__(self, **kwargs): - mesh_kwargs = self.mesh_kwargs.copy() - mesh_kwargs.update(kwargs) - - return globals()[self.mesh_name](**mesh_kwargs) - - class EllipseMeshBuilder(MeshBuilder): ambient_dim = 2 radius = 3.1 - aspect_ratio = 1.0 - resolutions = [8, 16, 32, 64] + aspect_ratio = 2.0 + resolutions = [16, 32, 64, 128] @property def curve_fn(self): @@ -69,6 +57,35 @@ class SphereMeshBuilder(MeshBuilder): 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 diff --git a/test/test_grudge.py b/test/test_grudge.py index a7550d69..6da37954 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 @@ -155,40 +154,103 @@ def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): assert err_3 < 5.0e-10, err_3 -@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 d6eab5b21bc4a51e5ee4e714e0900401e8340588 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 May 2020 21:17:10 -0500 Subject: [PATCH 4/8] use elementwise reductions --- test/test_grudge.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 6da37954..374ddcf6 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -188,6 +188,8 @@ 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) @@ -204,6 +206,10 @@ def test_mass_surface(ctx_factory, name): else: raise ValueError("unknown geometry name: %s" % name) + # }}} + + # {{{ convergence + from pytools.convergence import EOCRecorder eoc_area = EOCRecorder() eoc_inv = EOCRecorder() @@ -218,7 +224,8 @@ def test_mass_surface(ctx_factory, name): logger.info("bbox: %s", [(np.min(x[n]), np.max(x[n])) for n in range(x.shape[0])]) - # compute surface area + # {{{ 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) @@ -227,7 +234,10 @@ def test_mass_surface(ctx_factory, name): approx_surface_area, surface_area)) area_error = abs(approx_surface_area - surface_area) / abs(surface_area) - # compute inverse mass + # }}} + + # {{{ compute inverse mass + sym_op = sym.InverseMassOperator(dd, dd)( sym.MassOperator(dd, dd)(sym.Ones(dd))) approx_inverse = bind(discr, sym_op)(queue) @@ -236,14 +246,14 @@ def test_mass_surface(ctx_factory, name): cl.array.max(approx_inverse).get(queue), 1.0)) inv_error = la.norm(approx_inverse.get(queue) - 1.0) - # compute max element size - if resolution > 1: - h_max = 1.0 / (resolution + 1) - else: - h_max = resolution + # }}} + + h_max = bind(discr, sym.h_max(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)) -- GitLab From 9a13f276d964088379f9790bbdf6225b457832c8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 4 May 2020 10:50:10 -0500 Subject: [PATCH 5/8] implemented weight-adjusted inverse --- grudge/symbolic/mappers/__init__.py | 35 ++++++++++++++++++++--------- test/test_grudge.py | 2 +- 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 6e6e259e..814fc47f 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -650,14 +650,17 @@ 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) @@ -673,17 +676,27 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): 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.ambient_dim != self.dim: + # 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): @@ -695,12 +708,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/test_grudge.py b/test/test_grudge.py index 374ddcf6..d132cb17 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -248,7 +248,7 @@ def test_mass_surface(ctx_factory, name): # }}} - h_max = bind(discr, sym.h_max(dd))(queue) + h_max = bind(discr, sym.h_max(discr.ambient_dim, dd=dd))(queue) eoc_area.add_data_point(h_max, area_error) eoc_inv.add_data_point(h_max, inv_error) -- GitLab From 45db11507cc3ef3674e7c6e220e41c416bb7b2d5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 5 May 2020 21:25:25 -0500 Subject: [PATCH 6/8] add test with curvilinear mesh --- grudge/symbolic/mappers/__init__.py | 23 +++++++++++++++-------- test/test_grudge.py | 14 ++++++++++++-- 2 files changed, 27 insertions(+), 10 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 814fc47f..bab27e72 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -632,15 +632,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 @@ -663,16 +667,19 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): 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): @@ -680,7 +687,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.InverseMassOperator): - if self.ambient_dim != self.dim: + 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)( diff --git a/test/test_grudge.py b/test/test_grudge.py index d132cb17..f82aca3d 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -351,6 +351,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"]) @@ -403,7 +404,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) @@ -461,7 +467,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type for event in dt_stepper.run(t_end=final_time): if isinstance(event, dt_stepper.StateComputed): step += 1 - # print(event.t) + print(event.t) last_t = event.t last_u = event.state_component @@ -480,7 +486,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]) -- GitLab From 2f245d671e6c969775fc1f5b9e744b4d516f8536 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 10 May 2020 11:15:16 -0500 Subject: [PATCH 7/8] tests: skip warped strong form test --- examples/dagrt-fusion.py | 1 + test/test_grudge.py | 3 +++ 2 files changed, 4 insertions(+) diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py index d4b977b1..a96923db 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/test/test_grudge.py b/test/test_grudge.py index f82aca3d..023687dc 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -363,6 +363,9 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type cl_ctx = cl.create_some_context() 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() -- GitLab From ce2158ca1c8950eb5690bdafde4d834dd76aac74 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 18 May 2020 10:07:26 -0500 Subject: [PATCH 8/8] renamed h_max to h_max_from_volume --- test/test_grudge.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 023687dc..ca9089b0 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -190,18 +190,21 @@ def test_mass_surface(ctx_factory, name): # {{{ cases - import mesh_data as data if name == "2-1-ellipse": - builder = data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + 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": - builder = data.SpheroidMeshBuilder() + from mesh_data import SpheroidMeshBuilder + builder = SpheroidMeshBuilder() surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) elif name == "box2d": - builder = data.BoxMeshBuilder(ambient_dim=2) + from mesh_data import BoxMeshBuilder + builder = BoxMeshBuilder(ambient_dim=2) surface_area = 1.0 elif name == "box3d": - builder = data.BoxMeshBuilder(ambient_dim=3) + from mesh_data import BoxMeshBuilder + builder = BoxMeshBuilder(ambient_dim=3) surface_area = 1.0 else: raise ValueError("unknown geometry name: %s" % name) @@ -248,7 +251,7 @@ def test_mass_surface(ctx_factory, name): # }}} - h_max = bind(discr, sym.h_max(discr.ambient_dim, dd=dd))(queue) + 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) -- GitLab