From 6601f91ae6ed31ee1eae46d8c04dc0d0487fa74f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 4 Jun 2020 16:03:09 -0500 Subject: [PATCH 01/29] make the inverse mass operator work on surfaces --- examples/dagrt-fusion.py | 1 + grudge/execution.py | 6 +- grudge/symbolic/compiler.py | 2 +- grudge/symbolic/mappers/__init__.py | 45 ++++--- test/mesh_data.py | 130 +++++++++++++++++++++ test/test_grudge.py | 175 ++++++++++++++++++++++++++-- 6 files changed, 334 insertions(+), 25 deletions(-) create mode 100644 test/mesh_data.py 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/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/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 36b57d50..eef6cce7 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -597,6 +597,9 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): self.ambient_dim = ambient_dim self.dim = dim + # NOTE: only use WADG on surfaces at the moment + self.use_wadg = self.ambient_dim == (self.dim + 1) + map_common_subexpression_uncached = \ IdentityMapper.map_common_subexpression @@ -605,14 +608,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) @@ -628,21 +634,31 @@ 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.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): - return op.RefMassOperator()( + return op.RefMassOperator(dd_in=dd_in, dd_out=dd_out)( jac_noquad * self.rec( op.DiffOperator(expr.op.xyz_axis)(expr.field))) @@ -650,12 +666,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( @@ -805,6 +821,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): @@ -1174,7 +1193,7 @@ class ErrorChecker(CSECachingMapperMixin, IdentityMapper): def map_operator_binding(self, expr): if isinstance(expr.op, op.DiffOperatorBase): if (self.mesh is not None - and expr.op.xyz_axis >= self.mesh.dim): + and expr.op.xyz_axis >= self.mesh.ambient_dim): raise ValueError("optemplate tries to differentiate along a " "non-existent axis (e.g. Z in 2D)") diff --git a/test/mesh_data.py b/test/mesh_data.py new file mode 100644 index 00000000..c6ba21ee --- /dev/null +++ b/test/mesh_data.py @@ -0,0 +1,130 @@ +import six +import numpy as np + + +class MeshBuilder: + order = 4 + mesh_order = None + + def __init__(self, **kwargs): + for k, v in six.iteritems(kwargs): + setattr(self, k, v) + + if self.mesh_order is None: + self.mesh_order = 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 Curve2DMeshBuilder(MeshBuilder): + ambient_dim = 2 + resolutions = [16, 32, 64, 128] + + 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 EllipseMeshBuilder(Curve2DMeshBuilder): + radius = 3.1 + aspect_ratio = 2.0 + + @property + def curve_fn(self): + from meshmode.mesh.generation import ellipse + return lambda t: self.radius * ellipse(self.aspect_ratio, t) + + +class StarfishMeshBuilder(Curve2DMeshBuilder): + narms = 5 + amplitude = 0.25 + + @property + def curve_fn(self): + from meshmode.mesh.generation import NArmedStarfish + return NArmedStarfish(self.narms, self.amplitude) + + +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 + + mesh_order = 3 + resolutions = [1.0, 0.1, 0.05] + # resolutions = [1.0, 0.1, 0.05, 0.025, 0.015] + + @property + def radius(self): + return 0.25 + + @property + def diameter(self): + return 2.0 * self.radius + + @property + def aspect_ratio(self): + return 2.0 + + def get_mesh(self, resolution, mesh_order): + from meshmode.mesh.io import ScriptSource + source = ScriptSource(""" + SetFactory("OpenCASCADE"); + Sphere(1) = {{0, 0, 0, {r}}}; + Dilate {{ {{0, 0, 0}}, {{ {r}, {r}, {rr} }} }} {{ Volume{{1}}; }} + """.format(r=self.diameter, rr=self.aspect_ratio * self.diameter), + "geo" + ) + + from meshmode.mesh.io import generate_gmsh + mesh = generate_gmsh(source, 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 fc20811f..b41e003c 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -31,19 +31,17 @@ import pyopencl.array import pyopencl.clmath from pytools.obj_array import join_fields, make_obj_array +from grudge import sym, bind, DGDiscretizationWithBoundaries -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.setLevel(logging.INFO) - -from grudge import sym, bind, DGDiscretizationWithBoundaries +logging.basicConfig(level=logging.INFO) @pytest.mark.parametrize("dim", [2, 3]) @@ -157,6 +155,166 @@ 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 a < c: + 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_area(ctx_factory, name): + cl_ctx = cl.create_some_context() + 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 = 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("nelements: %d", volume_discr.mesh.nelements) + 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) + + # }}} + + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc.add_data_point(h_max, area_error) + + # }}} + + logger.info("surface area error\n%s", str(eoc)) + + assert eoc.max_error() < 1.0e-14 \ + or eoc.order_estimate() >= (builder.order - 1.0) + + +@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) +def test_surface_mass_operator_inverse(ctx_factory, name): + cl_ctx = cl.create_some_context() + 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) + else: + raise ValueError("unknown geometry name: %s" % name) + + # }}} + + # {{{ convergence + + from pytools.convergence import 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("nelements: %d", volume_discr.mesh.nelements) + logger.info("bbox: %s", + [(np.min(x[n]), np.max(x[n])) for n in range(x.shape[0])]) + + # {{{ compute inverse mass + + dd = sym.DD_VOLUME + sym_f = sym.cos(4.0 * sym.nodes(mesh.ambient_dim, dd)[0]) + sym_op = sym.InverseMassOperator(dd, dd)( + sym.MassOperator(dd, dd)(sym.var("f"))) + + f = bind(discr, sym_f)(queue) + f_inv = bind(discr, sym_op)(queue, f=f) + + logger.info("inverse: got {:.5e} / expected {:.5e}".format( + cl.array.max(f - f_inv).get(queue), 1.0)) + inv_error = la.norm(f.get(queue) - f_inv.get(queue)) \ + / la.norm(f.get(queue)) + + # }}} + + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc.add_data_point(h_max, inv_error) + + # }}} + + logger.info("inverse mass error\n%s", str(eoc)) + + assert eoc.max_error() < 5.0e-09 \ + 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 @@ -656,7 +814,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 a54e6ba89592a0cb34e886ce31cf295f61e57caa Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 4 Jun 2020 16:38:24 -0500 Subject: [PATCH 02/29] make face normals work on hypersurfaces --- grudge/symbolic/primitives.py | 63 ++++++++----- test/test_grudge.py | 161 ++++++++++++++++++++++++++++++++++ 2 files changed, 202 insertions(+), 22 deletions(-) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index d2f1a01f..06710462 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -538,7 +538,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( @@ -626,14 +628,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 @@ -644,23 +640,46 @@ 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 surface normal, move it to the + # face and then take a cross product with the face normal to get the + # correct exterior face normal vector. + assert dim == ambient_dim - 2 + + from grudge.symbolic.operators import interp + volm_normal = ( + surface_normal(ambient_dim, dim=dim + 1, dd=DD_VOLUME) + .map(interp(DD_VOLUME, dd))) + pder = pseudoscalar(ambient_dim, dim, dd=dd) + + mv = cse(-(volm_normal ^ pder) << volm_normal.I.inv(), + "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 b41e003c..3bf5ddf4 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -44,6 +44,8 @@ logger = logging.getLogger(__name__) logging.basicConfig(level=logging.INFO) +# {{{ inverse metric + @pytest.mark.parametrize("dim", [2, 3]) def test_inverse_metric(ctx_factory, dim): cl_ctx = ctx_factory() @@ -88,6 +90,10 @@ def test_inverse_metric(ctx_factory, dim): logger.info("error[%d, %d]: %.5e", i, j, err) assert err < 1.0e-12, (i, j, err) +# }}} + + +# {{{ mass operator trig integration @pytest.mark.parametrize("ambient_dim", [1, 2, 3]) @pytest.mark.parametrize("quad_tag", [sym.QTAG_NONE, "OVSMP"]) @@ -154,6 +160,10 @@ def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag): err_3 = abs(num_integral_3 - true_integral) assert err_3 < 5.0e-10, err_3 +# }}} + + +# {{{ mass operator surface area def _ellipse_surface_area(radius, aspect_ratio): # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html @@ -250,6 +260,10 @@ def test_mass_surface_area(ctx_factory, name): assert eoc.max_error() < 1.0e-14 \ or eoc.order_estimate() >= (builder.order - 1.0) +# }}} + + +# {{{ surface mass inverse @pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) def test_surface_mass_operator_inverse(ctx_factory, name): @@ -314,6 +328,119 @@ def test_surface_mass_operator_inverse(ctx_factory, name): assert eoc.max_error() < 5.0e-09 \ or eoc.order_estimate() >= (builder.order - 1.0) +# }}} + + +# {{{ surface face normal orthogonality + +def _avg_face_normal(dd, ambient_dim, dim=None): + dd = sym.as_dofdesc(dd) + assert dd.is_trace() + + if dim is None: + dim = ambient_dim - 1 + + face_normal_i = sym.normal(dd, ambient_dim, dim=dim) + face_normal_e = sym.OppositeInteriorFaceSwap()(face_normal_i) + face_normal = (face_normal_i - face_normal_e) / 2.0 + + return sym.cse( + face_normal / sym.sqrt(face_normal.dot(face_normal)), + "avg_face_normal", + sym.cse_scope.DISCRETIZATION) + + +@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) +def test_face_normal_surface(ctx_factory, mesh_name): + """Check that face normals are orthogonal to the surface normal""" + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # {{{ geometry + + if mesh_name == "2-1-ellipse": + from mesh_data import EllipseMeshBuilder + builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + elif mesh_name == "spheroid": + from mesh_data import SpheroidMeshBuilder + builder = SpheroidMeshBuilder() + else: + raise ValueError("unknown mesh name: %s" % mesh_name) + + mesh = builder.get_mesh(builder.resolutions[1], builder.mesh_order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) + + # }}} + + # {{{ symbolic + + dv = sym.DD_VOLUME + df = sym.as_dofdesc(sym.FACE_RESTR_INTERIOR) + + ambient_dim = mesh.ambient_dim + surf_dim = mesh.dim + face_dim = surf_dim - 1 + + sym_surf_normal = sym.interp(dv, df)( + sym.surface_normal(ambient_dim, dim=surf_dim, dd=dv).as_vector() + ) + sym_surf_normal = sym_surf_normal / sym.sqrt(sum(sym_surf_normal**2)) + + sym_face_normal_i = sym.normal(df, ambient_dim, dim=face_dim) + sym_face_normal_e = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_i) + + sym_face_normal_avg = _avg_face_normal(df, ambient_dim, dim=face_dim) + sym_face_normal_op = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_avg) + + if mesh.ambient_dim == 3: + # NOTE: there's only one face tangent in 3d + sym_face_tangent = (sym.pseudoscalar(ambient_dim, face_dim, dd=df) + / sym.area_element(ambient_dim, face_dim, dd=df)).as_vector() + + # }}} + + # {{{ checks + + rtol = 1.0e-14 + + surf_normal = bind(discr, sym_surf_normal)(queue) + + face_normal_i = bind(discr, sym_face_normal_i)(queue) + face_normal_e = bind(discr, sym_face_normal_e)(queue) + + face_normal_avg = bind(discr, sym_face_normal_avg)(queue) + face_normal_op = bind(discr, sym_face_normal_op)(queue) + + # check interpolated surface normal is orthogonal to face normal + error = la.norm(surf_normal.dot(face_normal_i).get(queue), np.inf) + logger.info("error[n_dot_i]: %.5e", error) + assert error < rtol + + # check angle between two neighboring elements + error = la.norm(face_normal_i.dot(face_normal_e).get(queue) + 1.0, np.inf) + logger.info("error[i_dot_e]: %.5e", error) + assert error > rtol + + # check uniqueness of normal on the two sides + error = la.norm(sum(face_normal_avg + face_normal_op).get(queue), np.inf) + logger.info("error[a_plus_o]: %.5e", error) + assert error < rtol + + # check orthogonality with face tangent + if ambient_dim == 3: + face_tangent = bind(discr, sym_face_tangent)(queue) + + error = la.norm(face_tangent.dot(face_normal_avg).get(queue), np.inf) + logger.info("error[t_dot_avg]: %.5e", error) + assert error < 5 * rtol + + # }}} + +# }}} + + +# {{{ diff operator @pytest.mark.parametrize("dim", [1, 2, 3]) def test_tri_diff_mat(ctx_factory, dim, order=4): @@ -354,6 +481,10 @@ def test_tri_diff_mat(ctx_factory, dim, order=4): logger.info("axis %d\n%s", axis, eoc_rec) assert eoc_rec.order_estimate() >= order +# }}} + + +# {{{ divergence theorem def test_2d_gauss_theorem(ctx_factory): """Verify Gauss's theorem explicitly on a mesh""" @@ -396,6 +527,10 @@ def test_2d_gauss_theorem(ctx_factory): assert abs(gauss_err) < 1e-13 +# }}} + + +# {{{ models: advection @pytest.mark.parametrize(("mesh_name", "mesh_pars"), [ ("segment", [8, 16, 32]), @@ -532,6 +667,10 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type assert eoc_rec.order_estimate() > order +# }}} + + +# {{{ models: maxwell @pytest.mark.parametrize("order", [3, 4, 5]) def test_convergence_maxwell(ctx_factory, order): @@ -602,6 +741,10 @@ def test_convergence_maxwell(ctx_factory, order): assert eoc_rec.order_estimate() > order +# }}} + + +# {{{ models: variable coefficient advection oversampling @pytest.mark.parametrize("order", [2, 3, 4]) def test_improvement_quadrature(ctx_factory, order): @@ -671,6 +814,10 @@ def test_improvement_quadrature(ctx_factory, order): assert (q_errs < errs).all() assert q_eoc > order +# }}} + + +# {{{ foreign points def test_foreign_points(ctx_factory): pytest.importorskip("sumpy") @@ -687,6 +834,10 @@ def test_foreign_points(ctx_factory): bind(pdiscr, sym.nodes(dim)**2)(queue) +# }}} + + +# {{{ operator collector determinism def test_op_collector_order_determinism(): class TestOperator(sym.Operator): @@ -712,6 +863,10 @@ def test_op_collector_order_determinism(): # The output order isn't significant, but it should always be the same. assert list(TestBoundOperatorCollector(TestOperator)(ob0 + ob1)) == [ob0, ob1] +# }}} + + +# {{{ bessel def test_bessel(ctx_factory): cl_ctx = ctx_factory() @@ -741,6 +896,10 @@ def test_bessel(ctx_factory): assert z < 1e-15 +# }}} + + +# {{{ function symbol def test_external_call(ctx_factory): cl_ctx = ctx_factory() @@ -805,6 +964,8 @@ def test_function_symbol_array(ctx_factory, array_type): norm = bind(discr, sym.norm(2, sym_x))(queue, x=x) assert isinstance(norm, float) +# }}} + # You can test individual routines by typing # $ python test_grudge.py 'test_routine()' -- GitLab From fc7a2a3507e3ab62fc832ba2f221a534910e6a25 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 4 Jun 2020 19:02:49 -0500 Subject: [PATCH 03/29] make stiffness operators work on hypersurfaces --- grudge/symbolic/mappers/__init__.py | 21 ++-- grudge/symbolic/primitives.py | 142 +++++++++++++++++++++-- test/mesh_data.py | 2 +- test/test_grudge.py | 172 +++++++++++++++++++++++++++- 4 files changed, 319 insertions(+), 18 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index eef6cce7..c8def898 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -621,17 +621,24 @@ 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) + def imd(rst): + return sym.inverse_surface_metric_derivative( + rst, expr.op.xyz_axis, + ambient_dim=self.ambient_dim, 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 - 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) - for rst_axis in range(self.dim)) + + return sum( + ref_class(rst_axis, dd_in=dd_in)(rec_field * imd(rst_axis)) + for rst_axis in range(self.dim)) + else: + return sum( + ref_class(rst_axis, dd_in=dd_in)(rec_field) * imd(rst_axis) + for rst_axis in range(self.dim)) if isinstance(expr.op, op.MassOperator): return op.RefMassOperator(dd_in, dd_out)( diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 06710462..fd963eba 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -494,7 +494,7 @@ def mv_nodes(ambient_dim, dd=None): return MultiVector(nodes(ambient_dim, dd)) -def forward_metric_derivative(xyz_axis, rst_axis, dd=None): +def forward_metric_nth_derivative(xyz_axis, ref_axes, dd=None): r""" Pointwise metric derivatives representing @@ -502,24 +502,44 @@ def forward_metric_derivative(xyz_axis, rst_axis, dd=None): \frac{\partial x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{rst\_axis}} } """ + + if isinstance(ref_axes, int): + ref_axes = ((ref_axes, 1),) + + if not isinstance(ref_axes, tuple): + raise ValueError("ref_axes must be a tuple") + + if tuple(sorted(ref_axes)) != ref_axes: + raise ValueError("ref_axes must be sorted") + + if len(dict(ref_axes)) != len(ref_axes): + raise ValueError("ref_axes must not contain an axis more than once") + if dd is None: dd = DD_VOLUME - inner_dd = dd.with_qtag(QTAG_NONE) - from grudge.symbolic.operators import RefDiffOperator - diff_op = RefDiffOperator(rst_axis, inner_dd) + from pytools import flatten + flat_ref_axes = flatten([rst_axis] * n for rst_axis, n in ref_axes) - result = diff_op(NodeCoordinateComponent(xyz_axis, inner_dd)) + from grudge.symbolic.operators import RefDiffOperator + result = NodeCoordinateComponent(xyz_axis, inner_dd) + for rst_axis in flat_ref_axes: + result = RefDiffOperator(rst_axis, inner_dd)(result) if dd.uses_quadrature(): from grudge.symbolic.operators import interp result = interp(inner_dd, dd)(result) - return cse( - result, - prefix="dx%d_dr%d" % (xyz_axis, rst_axis), - scope=cse_scope.DISCRETIZATION) + prefix = "dx%d_%s" % ( + xyz_axis, + "_".join("%sr%d" % ("d" * n, rst_axis) for rst_axis, n in ref_axes)) + + return cse(result, prefix, cse_scope.DISCRETIZATION) + + +def forward_metric_derivative(xyz_axis, rst_axis, dd=None): + return forward_metric_nth_derivative(xyz_axis, rst_axis, dd=dd) def forward_metric_derivative_vector(ambient_dim, rst_axis, dd=None): @@ -594,6 +614,7 @@ def forward_metric_derivative_mat(ambient_dim, dim=None, dd=None): result = np.zeros((ambient_dim, dim), dtype=np.object) for j in range(dim): result[:, j] = forward_metric_derivative_vector(ambient_dim, j, dd=dd) + return result @@ -610,6 +631,91 @@ def inverse_metric_derivative_mat(ambient_dim, dim=None, dd=None): return result +def first_fundamental_form(ambient_dim, dim=None, dd=None): + if dim is None: + dim = ambient_dim - 1 + + mder = forward_metric_derivative_mat(ambient_dim, dim=dim, dd=dd) + return cse(mder.T.dot(mder), "form1_mat", cse_scope.DISCRETIZATION) + + +def inverse_first_fundamental_form(ambient_dim, dim=None, dd=None): + if dim is None: + dim = ambient_dim - 1 + + if ambient_dim == dim: + inv_mder = inverse_metric_derivative_mat(ambient_dim, dim=dim, dd=dd) + inv_form1 = inv_mder.dot(inv_mder.T) + else: + form1 = first_fundamental_form(ambient_dim, dim, dd) + + if dim == 1: + inv_form1 = np.array([[1.0/form1[0, 0]]], dtype=np.object) + elif dim == 2: + (E, F), (_, G) = form1 # noqa: N806 + inv_form1 = 1.0 / (E * G - F * F) * np.array([ + [G, -F], [-F, E] + ], dtype=np.object) + else: + raise ValueError("%dD surfaces not supported" % dim) + + return cse(inv_form1, "inv_form1_mat", cse_scope.DISCRETIZATION) + + +def inverse_surface_metric_derivative(rst_axis, xyz_axis, + ambient_dim, dim=None, dd=None): + if dim is None: + dim = ambient_dim - 1 + + if ambient_dim == dim: + return inverse_metric_derivative(rst_axis, xyz_axis, + ambient_dim, dim=dim, dd=dd) + else: + inv_form1 = inverse_first_fundamental_form(ambient_dim, dim=dim, dd=dd) + imd = sum( + inv_form1[rst_axis, d] * forward_metric_derivative(xyz_axis, d, dd=dd) + for d in range(dim)) + + return cse(imd, + prefix="ds%d_dx%d" % (rst_axis, xyz_axis), + scope=cse_scope.DISCRETIZATION) + + +def second_fundamental_form(ambient_dim, dim=None, dd=None): + if dim is None: + dim = ambient_dim - 1 + + normal = surface_normal(ambient_dim, dim=dim, dd=dd).as_vector() + if dim == 1: + second_ref_axes = [((0, 2),)] + elif dim == 2: + second_ref_axes = [((0, 2),), ((0, 1), (1, 1)), ((1, 2),)] + else: + raise ValueError("%dD surfaces not supported" % dim) + + from pytools import flatten + form2 = np.empty((dim, dim), dtype=np.object) + for ref_axes in second_ref_axes: + i, j = flatten([rst_axis] * n for rst_axis, n in ref_axes) + + ruv = np.array([ + forward_metric_nth_derivative(xyz_axis, ref_axes, dd=dd) + for xyz_axis in range(ambient_dim)]) + form2[i, j] = form2[j, i] = normal.dot(ruv) + + return cse(form2, "form2_mat", cse_scope.DISCRETIZATION) + + +def shape_operator(ambient_dim, dim=None, dd=None): + if dim is None: + dim = ambient_dim - 1 + + inv_form1 = inverse_first_fundamental_form(ambient_dim, dim=dim, dd=dd) + form2 = second_fundamental_form(ambient_dim, dim=dim, dd=dd) + + return cse(-form2.dot(inv_form1), "shape_operator", cse_scope.DISCRETIZATION) + + def pseudoscalar(ambient_dim, dim=None, dd=None): if dim is None: dim = ambient_dim @@ -685,6 +791,24 @@ def mv_normal(dd, ambient_dim, dim=None): def normal(dd, ambient_dim, dim=None, quadrature_tag=None): return mv_normal(dd, ambient_dim, dim).as_vector() + +def summed_curvature(ambient_dim, dim=None, dd=None): + if dim is None: + dim = ambient_dim - 1 + + if ambient_dim == 1: + return 0.0 + + if ambient_dim == dim: + return 0.0 + + op = shape_operator(ambient_dim, dim=dim, dd=dd) + return cse(np.trace(op), "summed_curvature", cse_scope.DISCRETIZATION) + + +def mean_curvature(ambient_dim, dim=None, dd=None): + return 1.0 / (ambient_dim-1.0) * summed_curvature(ambient_dim, dim=dim, dd=dd) + # }}} diff --git a/test/mesh_data.py b/test/mesh_data.py index c6ba21ee..d1d04dac 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -74,7 +74,7 @@ class SpheroidMeshBuilder(MeshBuilder): mesh_order = 3 resolutions = [1.0, 0.1, 0.05] - # resolutions = [1.0, 0.1, 0.05, 0.025, 0.015] + # resolutions = [1.0, 0.1, 0.05, 0.03, 0.015] @property def radius(self): diff --git a/test/test_grudge.py b/test/test_grudge.py index 3bf5ddf4..e0d08dfc 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -527,6 +527,163 @@ def test_2d_gauss_theorem(ctx_factory): assert abs(gauss_err) < 1e-13 + +@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) +def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): + r"""Check the surface divergence theorem. + + .. math:: + + \int_Sigma \phi \nabla_i f_i = + \int_\Sigma \nabla_i \phi f_i + + \int_\Sigma \kappa \phi f_i n_i + + \int_{\partial \Sigma} \phi f_i m_i + + where :math:`n_i` is the surface normal and :class:`m_i` is the + face normal (which should be orthogonal to both the surface normal + and the face tangent). + """ + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # {{{ cases + + if mesh_name == "2-1-ellipse": + from mesh_data import EllipseMeshBuilder + builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + elif mesh_name == "spheroid": + from mesh_data import SpheroidMeshBuilder + builder = SpheroidMeshBuilder() + elif mesh_name == "circle": + from mesh_data import EllipseMeshBuilder + builder = EllipseMeshBuilder(radius=1.0, aspect_ratio=1.0) + elif mesh_name == "starfish": + from mesh_data import StarfishMeshBuilder + builder = StarfishMeshBuilder() + elif mesh_name == "sphere": + from mesh_data import SphereMeshBuilder + builder = SphereMeshBuilder(radius=1.0, mesh_order=16) + else: + raise ValueError("unknown mesh name: %s" % mesh_name) + + # }}} + + # {{{ convergene + + def f(x): + return join_fields( + sym.sin(3*x[1]) + sym.cos(3*x[0]) + 1.0, + sym.sin(2*x[0]) + sym.cos(x[1]), + 3.0 * sym.cos(x[0] / 2) + sym.cos(x[1]), + )[:ambient_dim] + + from pytools.convergence import EOCRecorder + eoc_global = EOCRecorder() + eoc_local = EOCRecorder() + + theta = np.pi / 3.33 + ambient_dim = builder.ambient_dim + if ambient_dim == 2: + mesh_rotation = np.array([ + [np.cos(theta), -np.sin(theta)], + [np.sin(theta), np.cos(theta)], + ]) + else: + mesh_rotation = np.array([ + [1.0, 0.0, 0.0], + [0.0, np.cos(theta), -np.sin(theta)], + [0.0, np.sin(theta), np.cos(theta)], + ]) + + mesh_offset = np.array([0.33, -0.21, 0.0])[:ambient_dim] + + for i, resolution in enumerate(builder.resolutions): + from meshmode.mesh.processing import affine_map + mesh = builder.get_mesh(resolution, builder.mesh_order) + # mesh = affine_map(mesh, A=mesh_rotation, b=mesh_offset) + + from meshmode.discretization.poly_element import \ + QuadratureSimplexGroupFactory + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order, + quad_tag_to_group_factory={ + "product": QuadratureSimplexGroupFactory(2 * builder.order) + }) + + volume = discr.discr_from_dd(sym.DD_VOLUME) + assert len(volume.groups) == 1 + logger.info("nnodes: %d", volume.nnodes) + logger.info("nelements: %d", volume.mesh.nelements) + + dd = sym.DD_VOLUME + dq = dd.with_qtag("product") + df = sym.as_dofdesc(sym.FACE_RESTR_ALL) + ambient_dim = discr.ambient_dim + dim = discr.dim + + # variables + sym_f = f(sym.nodes(ambient_dim, dd=dd)) + sym_f_quad = f(sym.nodes(ambient_dim, dd=dq)) + sym_kappa = sym.summed_curvature(ambient_dim, dim=dim, dd=dq) + sym_normal = sym.surface_normal(ambient_dim, dim=dim, dd=dq).as_vector() + + sym_face_normal = sym.normal(df, ambient_dim, dim=dim - 1) + sym_face_f = sym.interp(dd, df)(sym_f) + + # operators + sym_stiff = sum( + sym.StiffnessOperator(d)(f) for d, f in enumerate(sym_f) + ) + sym_stiff_t = sum( + sym.StiffnessTOperator(d)(f) for d, f in enumerate(sym_f) + ) + sym_k = sym.MassOperator(dq, dd)(sym_kappa * sym_f_quad.dot(sym_normal)) + sym_flux = sym.FaceMassOperator()(sym_face_f.dot(sym_face_normal)) + + # sum everything up and check the result + sym_op_global = sym.NodalSum(dd)( + sym_stiff - (sym_stiff_t + sym_k)) + sym_op_local = sym.ElementwiseSumOperator(dd)( + sym_stiff - (sym_stiff_t + sym_k + sym_flux)) + + op_global = bind(discr, sym_op_global)(queue) + op_local = bind(discr, sym_op_local)(queue).get(queue) + + err_global = la.norm(op_global) + err_local = la.norm( + la.norm(volume.groups[0].view(op_local), np.inf, axis=1), + np.inf) + logger.info("errors: global %.5e local %.5e", err_global, err_local) + + # compute max element size + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc_global.add_data_point(h_max, err_global) + eoc_local.add_data_point(h_max, err_local) + + if visualize: + r = cl.array.to_device(queue, op_local) + r = cl.clmath.log10(cl.clmath.fabs(r) + 1.0e-16) + + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(queue, discr, vis_order=builder.order) + + filename = "test_surface_divergence_theorem_error_{:04d}".format(i) + vis.write_vtk_file(filename, [ + ("r", r) + ], overwrite=True, legend=False) + + # }}} + + order = min(builder.order, builder.mesh_order) - 0.5 + logger.info("\n%s", str(eoc_global)) + logger.info("\n%s", str(eoc_local)) + + assert eoc_global.max_error() < 1.0e-12 \ + or eoc_global.order_estimate() >= order + + assert eoc_local.max_error() < 1.0e-12 \ + or eoc_local.order_estimate() >= order + # }}} @@ -537,6 +694,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"]) @@ -549,6 +707,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() @@ -590,7 +751,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) @@ -665,7 +831,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 # }}} -- GitLab From a3a23a9eecca8b701bd7b0bb63f47e5df5e6f278 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 4 Jun 2020 19:52:29 -0500 Subject: [PATCH 04/29] flake8 fixes --- grudge/symbolic/primitives.py | 2 +- test/mesh_data.py | 10 +++++----- test/test_grudge.py | 4 +--- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index fd963eba..ec75087e 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -673,7 +673,7 @@ def inverse_surface_metric_derivative(rst_axis, xyz_axis, else: inv_form1 = inverse_first_fundamental_form(ambient_dim, dim=dim, dd=dd) imd = sum( - inv_form1[rst_axis, d] * forward_metric_derivative(xyz_axis, d, dd=dd) + inv_form1[rst_axis, d]*forward_metric_derivative(xyz_axis, d, dd=dd) for d in range(dim)) return cse(imd, diff --git a/test/mesh_data.py b/test/mesh_data.py index d1d04dac..f508b520 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -91,11 +91,11 @@ class SpheroidMeshBuilder(MeshBuilder): def get_mesh(self, resolution, mesh_order): from meshmode.mesh.io import ScriptSource source = ScriptSource(""" - SetFactory("OpenCASCADE"); - Sphere(1) = {{0, 0, 0, {r}}}; - Dilate {{ {{0, 0, 0}}, {{ {r}, {r}, {rr} }} }} {{ Volume{{1}}; }} - """.format(r=self.diameter, rr=self.aspect_ratio * self.diameter), - "geo" + SetFactory("OpenCASCADE"); + Sphere(1) = {{0, 0, 0, {r}}}; + Dilate {{ {{0, 0, 0}}, {{ {r}, {r}, {rr} }} }} {{ Volume{{1}}; }} + """.format(r=self.diameter, rr=self.aspect_ratio * self.diameter), + "geo" ) from meshmode.mesh.io import generate_gmsh diff --git a/test/test_grudge.py b/test/test_grudge.py index e0d08dfc..366494ec 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -275,11 +275,9 @@ def test_surface_mass_operator_inverse(ctx_factory, name): 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) else: raise ValueError("unknown geometry name: %s" % name) @@ -601,7 +599,7 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): for i, resolution in enumerate(builder.resolutions): from meshmode.mesh.processing import affine_map mesh = builder.get_mesh(resolution, builder.mesh_order) - # mesh = affine_map(mesh, A=mesh_rotation, b=mesh_offset) + mesh = affine_map(mesh, A=mesh_rotation, b=mesh_offset) from meshmode.discretization.poly_element import \ QuadratureSimplexGroupFactory -- GitLab From 30d2f6119198bc52724a030c62c2a4898aff41fd Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 4 Jun 2020 20:08:18 -0500 Subject: [PATCH 05/29] add a dim to PointsDiscretization --- grudge/discretization.py | 5 +++++ test/mesh_data.py | 3 ++- test/test_grudge.py | 18 ++---------------- 3 files changed, 9 insertions(+), 17 deletions(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index d24fb2c0..7847c4ba 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -388,9 +388,14 @@ class PointsDiscretization(DiscretizationBase): self.mpi_communicator = None + @property def ambient_dim(self): return self._nodes.shape[0] + @property + def dim(self): + return self.ambient_dim + @property def mesh(self): return self diff --git a/test/mesh_data.py b/test/mesh_data.py index f508b520..c94e497c 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -102,7 +102,8 @@ class SpheroidMeshBuilder(MeshBuilder): mesh = generate_gmsh(source, 2, order=mesh_order, other_options=[ "-string", - "Mesh.CharacteristicLengthMax = %g;" % resolution]) + "Mesh.CharacteristicLengthMax = %g;" % resolution], + target_unit="MM") from meshmode.mesh.processing import perform_flips return perform_flips(mesh, np.ones(mesh.nelements)) diff --git a/test/test_grudge.py b/test/test_grudge.py index 366494ec..7d867470 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -455,7 +455,7 @@ def test_tri_diff_mat(ctx_factory, dim, order=4): from pytools.convergence import EOCRecorder axis_eoc_recs = [EOCRecorder() for axis in range(dim)] - for n in [10, 20]: + for n in [4, 8, 16]: mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, n=(n,)*dim, order=4) @@ -692,7 +692,6 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): ("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"]) @@ -705,9 +704,6 @@ 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() @@ -749,12 +745,6 @@ 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) @@ -829,11 +819,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type abscissa_label="h", error_label="L2 Error")) - 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 + assert eoc_rec.order_estimate() > order # }}} -- GitLab From 40d609d609864349f4b91ab6224dc2beb4159ed2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 8 Jun 2020 20:02:08 -0500 Subject: [PATCH 06/29] use correct h_max on surfaces --- test/mesh_data.py | 9 +++++---- test/test_grudge.py | 19 +++++++++++-------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/test/mesh_data.py b/test/mesh_data.py index c94e497c..9cc144aa 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -73,8 +73,8 @@ class SpheroidMeshBuilder(MeshBuilder): ambient_dim = 3 mesh_order = 3 - resolutions = [1.0, 0.1, 0.05] - # resolutions = [1.0, 0.1, 0.05, 0.03, 0.015] + resolutions = [1.0, 0.11, 0.05] + # resolutions = [1.0, 0.11, 0.05, 0.03, 0.015] @property def radius(self): @@ -101,8 +101,9 @@ class SpheroidMeshBuilder(MeshBuilder): from meshmode.mesh.io import generate_gmsh mesh = generate_gmsh(source, 2, order=mesh_order, other_options=[ - "-string", - "Mesh.CharacteristicLengthMax = %g;" % resolution], + "-optimize_ho", + "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution + ], target_unit="MM") from meshmode.mesh.processing import perform_flips diff --git a/test/test_grudge.py b/test/test_grudge.py index 7d867470..c8501d42 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -250,7 +250,8 @@ def test_mass_surface_area(ctx_factory, name): # }}} - h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + h_max = bind(discr, sym.h_max_from_volume( + discr.ambient_dim, dim=discr.dim, dd=dd))(queue) eoc.add_data_point(h_max, area_error) # }}} @@ -258,7 +259,7 @@ def test_mass_surface_area(ctx_factory, name): logger.info("surface area error\n%s", str(eoc)) assert eoc.max_error() < 1.0e-14 \ - or eoc.order_estimate() >= (builder.order - 1.0) + or eoc.order_estimate() > builder.order # }}} @@ -316,7 +317,8 @@ def test_surface_mass_operator_inverse(ctx_factory, name): # }}} - h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + h_max = bind(discr, sym.h_max_from_volume( + discr.ambient_dim, dim=discr.dim, dd=dd))(queue) eoc.add_data_point(h_max, inv_error) # }}} @@ -324,7 +326,7 @@ def test_surface_mass_operator_inverse(ctx_factory, name): logger.info("inverse mass error\n%s", str(eoc)) assert eoc.max_error() < 5.0e-09 \ - or eoc.order_estimate() >= (builder.order - 1.0) + or eoc.order_estimate() > builder.order # }}} @@ -477,7 +479,7 @@ def test_tri_diff_mat(ctx_factory, dim, order=4): for axis, eoc_rec in enumerate(axis_eoc_recs): logger.info("axis %d\n%s", axis, eoc_rec) - assert eoc_rec.order_estimate() >= order + assert eoc_rec.order_estimate() > order # }}} @@ -654,7 +656,8 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): logger.info("errors: global %.5e local %.5e", err_global, err_local) # compute max element size - h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + h_max = bind(discr, sym.h_max_from_volume( + discr.ambient_dim, dim=discr.dim, dd=dd))(queue) eoc_global.add_data_point(h_max, err_global) eoc_local.add_data_point(h_max, err_local) @@ -677,10 +680,10 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): logger.info("\n%s", str(eoc_local)) assert eoc_global.max_error() < 1.0e-12 \ - or eoc_global.order_estimate() >= order + or eoc_global.order_estimate() > order - 0.5 assert eoc_local.max_error() < 1.0e-12 \ - or eoc_local.order_estimate() >= order + or eoc_local.order_estimate() > order - 0.5 # }}} -- GitLab From 1a4d368e84b29195b1de4fd3f15de924b48c3eb2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 19 Jun 2020 18:50:21 -0500 Subject: [PATCH 07/29] add a h_min value as well --- grudge/symbolic/operators.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index c4afc0ec..77f735b9 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -742,6 +742,26 @@ def h_max_from_volume(ambient_dim, dim=None, dd=None): ) )**(1.0/dim) +def h_min_from_volume(ambient_dim, dim=None, dd=None): + """Defines a characteristic length based on the volume of the elements. + This length may not be representative if the elements have very high + aspect ratios. + """ + + import grudge.symbolic.primitives as prim + if dd is None: + dd = prim.DD_VOLUME + dd = prim.as_dofdesc(dd) + + if dim is None: + dim = ambient_dim + + return NodalMin(dd_in=dd)( + ElementwiseSumOperator(dd)( + MassOperator(dd_in=dd)(prim.Ones(dd)) + ) + )**(1.0/dim) + # }}} # vim: foldmethod=marker -- GitLab From e050cd087b1bc2283c9aa851520737e661ca98b8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 19 Jun 2020 18:51:15 -0500 Subject: [PATCH 08/29] add some more tests for the surface face normal --- test/mesh_data.py | 2 +- test/test_grudge.py | 52 ++++++++++++++++++++++++++------------------- 2 files changed, 31 insertions(+), 23 deletions(-) diff --git a/test/mesh_data.py b/test/mesh_data.py index 9cc144aa..2ea2cc47 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -72,7 +72,7 @@ class SphereMeshBuilder(MeshBuilder): class SpheroidMeshBuilder(MeshBuilder): ambient_dim = 3 - mesh_order = 3 + mesh_order = 4 resolutions = [1.0, 0.11, 0.05] # resolutions = [1.0, 0.11, 0.05, 0.03, 0.015] diff --git a/test/test_grudge.py b/test/test_grudge.py index c8501d42..e27fe03b 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -333,21 +333,12 @@ def test_surface_mass_operator_inverse(ctx_factory, name): # {{{ surface face normal orthogonality -def _avg_face_normal(dd, ambient_dim, dim=None): - dd = sym.as_dofdesc(dd) - assert dd.is_trace() +def _avg_face_normal(x, side=-1): + x_i = x + x_e = sym.OppositeInteriorFaceSwap()(x_i) + x_avg = (x_i + side * x_e) / 2.0 - if dim is None: - dim = ambient_dim - 1 - - face_normal_i = sym.normal(dd, ambient_dim, dim=dim) - face_normal_e = sym.OppositeInteriorFaceSwap()(face_normal_i) - face_normal = (face_normal_i - face_normal_e) / 2.0 - - return sym.cse( - face_normal / sym.sqrt(face_normal.dot(face_normal)), - "avg_face_normal", - sym.cse_scope.DISCRETIZATION) + return x_avg / sym.sqrt(x_avg.dot(x_avg)) @pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) @@ -368,9 +359,13 @@ def test_face_normal_surface(ctx_factory, mesh_name): else: raise ValueError("unknown mesh name: %s" % mesh_name) - mesh = builder.get_mesh(builder.resolutions[1], builder.mesh_order) + mesh = builder.get_mesh(builder.resolutions[0], builder.mesh_order) discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + logger.info("nnodes: %d", volume_discr.nnodes) + logger.info("nelements: %d", volume_discr.mesh.nelements) + # }}} # {{{ symbolic @@ -379,24 +374,26 @@ def test_face_normal_surface(ctx_factory, mesh_name): df = sym.as_dofdesc(sym.FACE_RESTR_INTERIOR) ambient_dim = mesh.ambient_dim - surf_dim = mesh.dim - face_dim = surf_dim - 1 + dim = mesh.dim sym_surf_normal = sym.interp(dv, df)( - sym.surface_normal(ambient_dim, dim=surf_dim, dd=dv).as_vector() + sym.surface_normal(ambient_dim, dim=dim, dd=dv).as_vector() ) sym_surf_normal = sym_surf_normal / sym.sqrt(sum(sym_surf_normal**2)) - sym_face_normal_i = sym.normal(df, ambient_dim, dim=face_dim) + sym_surf_normal_avg = _avg_face_normal(sym_surf_normal, side=1.0) + + sym_face_normal_i = sym.normal(df, ambient_dim, dim=dim - 1) sym_face_normal_e = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_i) - sym_face_normal_avg = _avg_face_normal(df, ambient_dim, dim=face_dim) + sym_face_normal_avg = _avg_face_normal(sym_face_normal_i) sym_face_normal_op = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_avg) if mesh.ambient_dim == 3: # NOTE: there's only one face tangent in 3d - sym_face_tangent = (sym.pseudoscalar(ambient_dim, face_dim, dd=df) - / sym.area_element(ambient_dim, face_dim, dd=df)).as_vector() + sym_face_tangent = ( + sym.pseudoscalar(ambient_dim, dim - 1, dd=df) + / sym.area_element(ambient_dim, dim - 1, dd=df)).as_vector() # }}} @@ -405,6 +402,7 @@ def test_face_normal_surface(ctx_factory, mesh_name): rtol = 1.0e-14 surf_normal = bind(discr, sym_surf_normal)(queue) + surf_normal_avg = bind(discr, sym_surf_normal_avg)(queue) face_normal_i = bind(discr, sym_face_normal_i)(queue) face_normal_e = bind(discr, sym_face_normal_e)(queue) @@ -417,6 +415,16 @@ def test_face_normal_surface(ctx_factory, mesh_name): logger.info("error[n_dot_i]: %.5e", error) assert error < rtol + # check averaged ones are also orthogonal + error = la.norm(surf_normal_avg.dot(face_normal_avg).get(queue), np.inf) + logger.info("error[a_dot_a]: %.5e", error) + assert error < rtol + + # check averaged face normal and interpolated face normal + error = la.norm(surf_normal.dot(face_normal_avg).get(queue), np.inf) + logger.info("error[n_dot_a]: %.5e", error) + assert error > rtol + # check angle between two neighboring elements error = la.norm(face_normal_i.dot(face_normal_e).get(queue) + 1.0, np.inf) logger.info("error[i_dot_e]: %.5e", error) -- GitLab From 5397437bddb516b9f5d8c78a1744e773f736d16e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Jun 2020 19:07:34 -0500 Subject: [PATCH 09/29] add a surface advection example --- examples/advection/surface.py | 211 ++++++++++++++++++++++++++++++++++ grudge/models/advection.py | 61 ++++++++++ 2 files changed, 272 insertions(+) create mode 100644 examples/advection/surface.py diff --git a/examples/advection/surface.py b/examples/advection/surface.py new file mode 100644 index 00000000..69ab3a90 --- /dev/null +++ b/examples/advection/surface.py @@ -0,0 +1,211 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2020 Alexandru Fikl" + +__license__ = """ +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +""" + +import os +import numpy as np +import numpy.linalg as la + +import pyopencl as cl +import pyopencl.array +import pyopencl.clmath + +from grudge import bind, sym +from pytools.obj_array import make_obj_array + +import logging +logger = logging.getLogger(__name__) + + +# {{{ plotting (keep in sync with `var-velocity.py`) + +class Plotter: + def __init__(self, queue, discr, order, visualize=True): + self.queue = queue + self.ambient_dim = discr.ambient_dim + self.dim = discr.dim + + self.visualize = visualize + if not self.visualize: + return + + if self.ambient_dim == 2: + import matplotlib.pyplot as pt + self.fig = pt.figure(figsize=(8, 8), dpi=300) + + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + x = volume_discr.nodes().with_queue(queue) + self.x = (2.0 * np.pi * (x[1] < 0) + cl.clmath.atan2(x[1], x[0])).get(queue) + elif self.ambient_dim == 3: + from grudge.shortcuts import make_visualizer + self.vis = make_visualizer(discr, vis_order=order) + else: + raise ValueError("unsupported dimension") + + def __call__(self, evt, basename, overwrite=True): + if not self.visualize: + return + + if self.ambient_dim == 2: + u = evt.state_component.get(self.queue) + + filename = "%s.png" % basename + if not overwrite and os.path.exists(filename): + from meshmode import FileExistsError + raise FileExistsError("output file '%s' already exists" % filename) + + ax = self.fig.gca() + ax.grid() + ax.plot(self.x, u, '-') + ax.plot(self.x, u, 'k.') + ax.set_xlabel(r"$\theta$") + ax.set_ylabel("$u$") + ax.set_title("t = {:.2f}".format(evt.t)) + + self.fig.savefig(filename) + self.fig.clf() + elif self.ambient_dim == 3: + self.vis.write_vtk_file("%s.vtu" % basename, [ + ("u", evt.state_component) + ], overwrite=overwrite) + else: + raise ValueError("unsupported dimension") + +# }}} + + +def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # {{{ parameters + + # sphere radius + radius = 1.0 + # sphere resolution + resolution = 64 if dim == 2 else 1 + + # cfl + dt_factor = 1.0 + # final time + final_time = 1.0 + + # velocity field + sym_x = sym.nodes(dim) + c = make_obj_array([ + -sym_x[1], sym_x[0], 0.0 + ])[:dim] + norm_c = sym.sqrt((c**2).sum()) + # flux + flux_type = "central" + + # }}} + + # {{{ discretization + + if dim == 2: + from meshmode.mesh.generation import make_curve_mesh, ellipse + mesh = make_curve_mesh( + lambda t: radius * ellipse(1.0, t), + np.linspace(0.0, 1.0, resolution + 1), + order) + elif dim == 3: + from meshmode.mesh.generation import generate_icosphere + return generate_icosphere(radius, order=order, + uniform_refinement_rounds=resolution) + else: + raise ValueError("unsupported dimension") + + quad_tag_to_group_factory = {} + if product_tag == "none": + product_tag = None + + if product_tag: + from meshmode.discretization.poly_element import \ + QuadratureSimplexGroupFactory + quad_tag_to_group_factory = { + product_tag: QuadratureSimplexGroupFactory(order=4*order) + } + + from grudge import DGDiscretizationWithBoundaries + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + quad_tag_to_group_factory=quad_tag_to_group_factory) + + # }}} + + # {{{ symbolic operators + + def f_gaussian(x): + return x[0] + + from grudge.models.advection import SurfaceAdvectionOperator + op = SurfaceAdvectionOperator(c, + flux_type=flux_type, + quad_tag=product_tag) + + bound_op = bind(discr, op.sym_operator()) + u = bind(discr, f_gaussian(sym_x))(queue, t=0) + + def rhs(t, u): + return bound_op(queue, t=t, u=u) + + # }}} + + # {{{ time stepping + + # compute time step + h_min = bind(discr, + sym.h_max_from_volume(discr.ambient_dim, dim=discr.dim))(queue) + dt = dt_factor * h_min/order**2 + nsteps = int(final_time // dt) + 1 + dt = final_time/nsteps + 1.0e-15 + + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("u", dt, u, rhs) + plot = Plotter(queue, discr, order, visualize=visualize) + + norm = bind(discr, sym.norm(2, sym.var("u"))) + + step = 0 + norm_u = 0.0 + for event in dt_stepper.run(t_end=final_time): + if not isinstance(event, dt_stepper.StateComputed): + continue + + if step % 1 == 0: + norm_u = norm(queue, u=event.state_component) + plot(event, "fld-surface-%04d" % step) + + step += 1 + logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + + # }}} + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("--dim", default=2, type=int) + parser.add_argument("--qtag", default="none") + args = parser.parse_args() + + logging.basicConfig(level=logging.INFO) + main(cl.create_some_context, + dim=args.dim, + product_tag=args.qtag) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 36ce5c15..ae3d8a33 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -182,4 +182,65 @@ class VariableCoefficientAdvectionOperator(AdvectionOperatorBase): )) # }}} + +# {{{ closed surface advection + +def v_dot_n_tpair(velocity, dd=None): + ambient_dim = len(velocity) + normal = sym.normal(dd, ambient_dim, dim=ambient_dim - 2) + + v_dot_n_i = sym.cse(velocity.dot(normal), "v_dot_normal") + v_dot_n_e = sym.cse(sym.OppositeInteriorFaceSwap()(v_dot_n_i)) + + return sym.TracePair(dd, v_dot_n_i, -v_dot_n_e) + + +def surface_advection_weak_flux(flux_type, u, velocity): + v_dot_n = v_dot_n_tpair(velocity, dd=u.dd) + + flux_type = flux_type.lower() + if flux_type == "central": + return u.avg * v_dot_n.avg + elif flux_type == "lf": + return u.avg * v_dot_n.avg + 0.5 * sym.fabs(v_dot_n) * (u.int - u.ext) + else: + raise ValueError("flux `{}` is not implemented".format(flux_type)) + + +class SurfaceAdvectionOperator(AdvectionOperatorBase): + def __init__(self, v, flux_type="central", quad_tag=None): + super(SurfaceAdvectionOperator, self).__init__( + v, inflow_u=None, flux_type=flux_type) + self.quad_tag = quad_tag + + def flux(self, u): + surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v) + return surface_advection_weak_flux(self.flux_type, u, surf_v) + + def sym_operator(self): + u = sym.var("u") + + def flux(pair): + return sym.interp(pair.dd, face_dd)(self.flux(pair)) + + face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) + + quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag) + to_quad = sym.interp(sym.DD_VOLUME, quad_dd) + stiff_t_op = sym.stiffness_t(self.ambient_dim, + dd_in=quad_dd, dd_out=sym.DD_VOLUME) + + quad_v = to_quad(self.v) + quad_u = to_quad(u) + + return sym.InverseMassOperator()( + sum(stiff_t_op[n](quad_u * quad_v[n]) + for n in range(self.ambient_dim)) + - sym.FaceMassOperator(face_dd, sym.DD_VOLUME)( + flux(sym.int_tpair(u, self.quad_tag)) + ) + ) + +# }}} + # vim: foldmethod=marker -- GitLab From 71a40a838d4e628e0645286ea9548dec367f6e32 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 21 Jun 2020 18:48:05 -0500 Subject: [PATCH 10/29] make surface example actually work properly --- examples/advection/surface.py | 80 ++++++++++++++++++++++++++++++----- grudge/models/advection.py | 27 ++++++++++-- grudge/symbolic/primitives.py | 4 ++ 3 files changed, 96 insertions(+), 15 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 69ab3a90..2ef5d770 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -98,12 +98,12 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): # sphere radius radius = 1.0 # sphere resolution - resolution = 64 if dim == 2 else 1 + resolution = 64 if dim == 2 else 1.0 # cfl dt_factor = 1.0 # final time - final_time = 1.0 + final_time = 2.0 * np.pi # velocity field sym_x = sym.nodes(dim) @@ -112,7 +112,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): ])[:dim] norm_c = sym.sqrt((c**2).sum()) # flux - flux_type = "central" + flux_type = "lf" # }}} @@ -125,9 +125,29 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): np.linspace(0.0, 1.0, resolution + 1), order) elif dim == 3: - from meshmode.mesh.generation import generate_icosphere - return generate_icosphere(radius, order=order, - uniform_refinement_rounds=resolution) + if 0: + from meshmode.mesh.generation import generate_icosphere + mesh = generate_icosphere(radius, order=4 * order) + from meshmode.mesh.refinement import refine_uniformly + mesh = refine_uniformly(mesh, resolution, with_adjacency=False) + else: + from meshmode.mesh.io import ScriptSource + source = ScriptSource(""" + SetFactory("OpenCASCADE"); + Sphere(1) = {0, 0, 0, %g}; + Mesh.Algorithm = 6; + """ % radius, "geo") + + from meshmode.mesh.io import generate_gmsh + mesh = generate_gmsh(source, 2, order=order, + other_options=[ + "-optimize_ho", + "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution + ], + target_unit="MM") + + from meshmode.mesh.processing import perform_flips + mesh = perform_flips(mesh, np.ones(mesh.nelements)) else: raise ValueError("unsupported dimension") @@ -146,10 +166,15 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_tag_to_group_factory=quad_tag_to_group_factory) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + logger.info("nnodes: %d", volume_discr.nnodes) + logger.info("nelements: %d", volume_discr.mesh.nelements) + # }}} # {{{ symbolic operators + def f_gaussian(x): return x[0] @@ -164,6 +189,11 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): def rhs(t, u): return bound_op(queue, t=t, u=u) + # check velocity is tangential + sym_normal = sym.surface_normal(dim, dim=dim - 1, dd=sym.DD_VOLUME).as_vector() + error = la.norm(bind(discr, c.dot(sym_normal))(queue).get(queue)) + logger.info("u_dot_n: %.5e", error) + # }}} # {{{ time stepping @@ -175,23 +205,51 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): nsteps = int(final_time // dt) + 1 dt = final_time/nsteps + 1.0e-15 + logger.info("dt: %.5e", dt) + logger.info("nsteps: %d", nsteps) + from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) plot = Plotter(queue, discr, order, visualize=visualize) norm = bind(discr, sym.norm(2, sym.var("u"))) + norm_u = 0.0 step = 0 - norm_u = 0.0 - for event in dt_stepper.run(t_end=final_time): + + event = dt_stepper.StateComputed(0.0, 0, 0, u) + plot(event, "fld-surface-%04d" % 0) + + if visualize and dim == 3: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + vis.write_vtk_file("fld-surface-velocity.vtu", [ + ("u", bind(discr, c)(queue)), + ("n", bind(discr, sym_normal)(queue)) + ], overwrite=True) + + df = sym.DOFDesc(sym.FACE_RESTR_INTERIOR) + face_discr = discr.connection_from_dds(sym.DD_VOLUME, df).to_discr + + face_normal = bind(discr, sym.normal( + df, face_discr.ambient_dim, dim=face_discr.dim))(queue) + + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(queue, face_discr, vis_order=order, + use_high_order_vtk=False) + vis.write_vtk_file("fld-surface-face-normals.vtu", [ + ("n", face_normal) + ], overwrite=True) + + for event in dt_stepper.run(t_end=final_time, max_steps=nsteps + 1): if not isinstance(event, dt_stepper.StateComputed): continue + step += 1 if step % 1 == 0: norm_u = norm(queue, u=event.state_component) plot(event, "fld-surface-%04d" % step) - step += 1 logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) # }}} @@ -201,8 +259,8 @@ if __name__ == "__main__": import argparse parser = argparse.ArgumentParser() - parser.add_argument("--dim", default=2, type=int) - parser.add_argument("--qtag", default="none") + parser.add_argument("--dim", choices=[2, 3], default=2, type=int) + parser.add_argument("--qtag", choices=["none", "product"], default="none") args = parser.parse_args() logging.basicConfig(level=logging.INFO) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index ae3d8a33..3891c9ca 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -189,24 +189,35 @@ def v_dot_n_tpair(velocity, dd=None): ambient_dim = len(velocity) normal = sym.normal(dd, ambient_dim, dim=ambient_dim - 2) - v_dot_n_i = sym.cse(velocity.dot(normal), "v_dot_normal") + v_dot_n_i = sym.cse(velocity.dot(normal)) v_dot_n_e = sym.cse(sym.OppositeInteriorFaceSwap()(v_dot_n_i)) - return sym.TracePair(dd, v_dot_n_i, -v_dot_n_e) + return sym.TracePair(dd, v_dot_n_i, v_dot_n_e) def surface_advection_weak_flux(flux_type, u, velocity): v_dot_n = v_dot_n_tpair(velocity, dd=u.dd) + v_dot_n = sym.cse(0.5 * v_dot_n.jump, "v_dot_normal") flux_type = flux_type.lower() if flux_type == "central": - return u.avg * v_dot_n.avg + return u.avg * v_dot_n elif flux_type == "lf": - return u.avg * v_dot_n.avg + 0.5 * sym.fabs(v_dot_n) * (u.int - u.ext) + return u.avg * v_dot_n + 0.5 * sym.fabs(v_dot_n) * (u.int - u.ext) else: raise ValueError("flux `{}` is not implemented".format(flux_type)) +def surface_penalty_flux(u, velocity, tau=1.0): + if abs(tau) < 1.0e-14: + return 0.0 + + v_dot_n = v_dot_n_tpair(velocity, dd=u.dd) + return sym.If(sym.Comparison(v_dot_n.avg, ">", 0), + 0.5 * tau * u.avg * v_dot_n.int, + 0.0) + + class SurfaceAdvectionOperator(AdvectionOperatorBase): def __init__(self, v, flux_type="central", quad_tag=None): super(SurfaceAdvectionOperator, self).__init__( @@ -217,12 +228,19 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase): surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v) return surface_advection_weak_flux(self.flux_type, u, surf_v) + def penalty(self, u): + surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v) + return surface_penalty_flux(u, surf_v, tau=0.0) + def sym_operator(self): u = sym.var("u") def flux(pair): return sym.interp(pair.dd, face_dd)(self.flux(pair)) + def penalty(pair): + return sym.interp(pair.dd, face_dd)(self.penalty(pair)) + face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag) @@ -238,6 +256,7 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase): for n in range(self.ambient_dim)) - sym.FaceMassOperator(face_dd, sym.DD_VOLUME)( flux(sym.int_tpair(u, self.quad_tag)) + + penalty(sym.int_tpair(u, self.quad_tag)) ) ) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index ec75087e..ceef7a10 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -858,6 +858,10 @@ class TracePair: def ext(self): return self.exterior + @property + def jump(self): + return self.int - self.ext + @property def avg(self): return 0.5*(self.int + self.ext) -- GitLab From 27d0a8256e5ac865c068bcc997ec2709cf82647b Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 21 Jun 2020 19:29:31 -0500 Subject: [PATCH 11/29] make example work with oversampling --- examples/advection/surface.py | 33 ++++++--------------------------- grudge/models/advection.py | 18 ++++++++++-------- grudge/symbolic/operators.py | 1 + grudge/symbolic/primitives.py | 23 ++++++++++++++--------- 4 files changed, 31 insertions(+), 44 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 2ef5d770..a198dabf 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -50,7 +50,8 @@ class Plotter: volume_discr = discr.discr_from_dd(sym.DD_VOLUME) x = volume_discr.nodes().with_queue(queue) - self.x = (2.0 * np.pi * (x[1] < 0) + cl.clmath.atan2(x[1], x[0])).get(queue) + self.x = (2.0 * np.pi * (x[1] < 0) + + cl.clmath.atan2(x[1], x[0])).get(queue) elif self.ambient_dim == 3: from grudge.shortcuts import make_visualizer self.vis = make_visualizer(discr, vis_order=order) @@ -98,7 +99,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): # sphere radius radius = 1.0 # sphere resolution - resolution = 64 if dim == 2 else 1.0 + resolution = 64 if dim == 2 else 1 # cfl dt_factor = 1.0 @@ -110,7 +111,6 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): c = make_obj_array([ -sym_x[1], sym_x[0], 0.0 ])[:dim] - norm_c = sym.sqrt((c**2).sum()) # flux flux_type = "lf" @@ -125,29 +125,9 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): np.linspace(0.0, 1.0, resolution + 1), order) elif dim == 3: - if 0: - from meshmode.mesh.generation import generate_icosphere - mesh = generate_icosphere(radius, order=4 * order) - from meshmode.mesh.refinement import refine_uniformly - mesh = refine_uniformly(mesh, resolution, with_adjacency=False) - else: - from meshmode.mesh.io import ScriptSource - source = ScriptSource(""" - SetFactory("OpenCASCADE"); - Sphere(1) = {0, 0, 0, %g}; - Mesh.Algorithm = 6; - """ % radius, "geo") - - from meshmode.mesh.io import generate_gmsh - mesh = generate_gmsh(source, 2, order=order, - other_options=[ - "-optimize_ho", - "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution - ], - target_unit="MM") - - from meshmode.mesh.processing import perform_flips - mesh = perform_flips(mesh, np.ones(mesh.nelements)) + from meshmode.mesh.generation import generate_icosphere + mesh = generate_icosphere(radius, order=4 * order, + uniform_refinement_rounds=resolution) else: raise ValueError("unsupported dimension") @@ -174,7 +154,6 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): # {{{ symbolic operators - def f_gaussian(x): return x[0] diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 3891c9ca..88c09668 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -186,13 +186,15 @@ class VariableCoefficientAdvectionOperator(AdvectionOperatorBase): # {{{ closed surface advection def v_dot_n_tpair(velocity, dd=None): - ambient_dim = len(velocity) - normal = sym.normal(dd, ambient_dim, dim=ambient_dim - 2) + if dd is None: + dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR) - v_dot_n_i = sym.cse(velocity.dot(normal)) - v_dot_n_e = sym.cse(sym.OppositeInteriorFaceSwap()(v_dot_n_i)) + ambient_dim = len(velocity) + normal = sym.normal(dd.with_qtag(None), ambient_dim, dim=ambient_dim - 2) - return sym.TracePair(dd, v_dot_n_i, v_dot_n_e) + return sym.int_tpair(velocity.dot(normal), + qtag=dd.quadrature_tag, + from_dd=dd.with_qtag(None)) def surface_advection_weak_flux(flux_type, u, velocity): @@ -225,11 +227,11 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase): self.quad_tag = quad_tag def flux(self, u): - surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v) + surf_v = sym.interp(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v) return surface_advection_weak_flux(self.flux_type, u, surf_v) def penalty(self, u): - surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v) + surf_v = sym.interp(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v) return surface_penalty_flux(u, surf_v, tau=0.0) def sym_operator(self): @@ -242,8 +244,8 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase): return sym.interp(pair.dd, face_dd)(self.penalty(pair)) face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) - quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag) + to_quad = sym.interp(sym.DD_VOLUME, quad_dd) stiff_t_op = sym.stiffness_t(self.ambient_dim, dd_in=quad_dd, dd_out=sym.DD_VOLUME) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 77f735b9..8f7609c0 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -742,6 +742,7 @@ def h_max_from_volume(ambient_dim, dim=None, dd=None): ) )**(1.0/dim) + def h_min_from_volume(ambient_dim, dim=None, dd=None): """Defines a characteristic length based on the volume of the elements. This length may not be representative if the elements have very high diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index ceef7a10..7e26223b 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -867,20 +867,25 @@ class TracePair: return 0.5*(self.int + self.ext) -def int_tpair(expression, qtag=None): +def int_tpair(expression, qtag=None, from_dd=None): from grudge.symbolic.operators import interp, OppositeInteriorFaceSwap - i = interp("vol", "int_faces")(expression) - e = cse(OppositeInteriorFaceSwap()(i)) + if from_dd is None: + from_dd = DD_VOLUME + assert not from_dd.uses_quadrature() - if qtag is not None and qtag != QTAG_NONE: - q_dd = DOFDesc("int_faces", qtag) - i = cse(interp("int_faces", q_dd)(i)) - e = cse(interp("int_faces", q_dd)(e)) + trace_dd = DOFDesc(FACE_RESTR_INTERIOR, qtag) + if from_dd.domain_tag == trace_dd.domain_tag: + i = expression else: - q_dd = "int_faces" + i = interp(from_dd, trace_dd.with_qtag(None))(expression) + e = cse(OppositeInteriorFaceSwap()(i)) + + if trace_dd.uses_quadrature(): + i = cse(interp(trace_dd.with_qtag(None), trace_dd)(i)) + e = cse(interp(trace_dd.with_qtag(None), trace_dd)(e)) - return TracePair(q_dd, i, e) + return TracePair(trace_dd, i, e) def bdry_tpair(dd, interior, exterior): -- GitLab From 2580b49f4415a78da3a4a854a4403b9a15575e58 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 26 Jun 2020 21:21:53 -0500 Subject: [PATCH 12/29] disable visualization in example --- examples/advection/surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index a198dabf..47092477 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -90,7 +90,7 @@ class Plotter: # }}} -def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): +def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) -- GitLab From f250eed315367a65a1777f3e28170bacaaa4d8ed Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 14 Jul 2020 20:45:45 -0500 Subject: [PATCH 13/29] fix visualization --- test/test_grudge.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 7d7ffbcc..fbc1f50f 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -673,16 +673,13 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): eoc_local.add_data_point(h_max, err_local) if visualize: - r = cl.array.to_device(queue, op_local) - r = cl.clmath.log10(cl.clmath.fabs(r) + 1.0e-16) + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=builder.order) - from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, vis_order=builder.order) - - filename = "test_surface_divergence_theorem_error_{:04d}".format(i) + filename = f"surface_divergence_theorem_{mesh_name}_{i:04d}.vtu" vis.write_vtk_file(filename, [ - ("r", r) - ], overwrite=True, legend=False) + ("r", actx.np.log10(op_local)) + ], overwrite=True) # }}} -- GitLab From 8d95e87f1cc930b35a9993e3eba22c130c03f0dc Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 14 Jul 2020 21:24:12 -0500 Subject: [PATCH 14/29] update advection examples to array-context --- examples/advection/surface.py | 51 +++++++++++++++--------------- examples/advection/var-velocity.py | 15 +++++---- examples/advection/weak.py | 15 +++++---- 3 files changed, 42 insertions(+), 39 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 47092477..424e100a 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -22,8 +22,9 @@ import numpy as np import numpy.linalg as la import pyopencl as cl -import pyopencl.array -import pyopencl.clmath + +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flatten from grudge import bind, sym from pytools.obj_array import make_obj_array @@ -35,8 +36,8 @@ logger = logging.getLogger(__name__) # {{{ plotting (keep in sync with `var-velocity.py`) class Plotter: - def __init__(self, queue, discr, order, visualize=True): - self.queue = queue + def __init__(self, actx, discr, order, visualize=True): + self.actx = actx self.ambient_dim = discr.ambient_dim self.dim = discr.dim @@ -48,10 +49,8 @@ class Plotter: import matplotlib.pyplot as pt self.fig = pt.figure(figsize=(8, 8), dpi=300) - volume_discr = discr.discr_from_dd(sym.DD_VOLUME) - x = volume_discr.nodes().with_queue(queue) - self.x = (2.0 * np.pi * (x[1] < 0) - + cl.clmath.atan2(x[1], x[0])).get(queue) + x = thaw(actx, discr.discr_from_dd(sym.DD_VOLUME).nodes()) + self.x = actx.to_numpy(flatten(actx.np.atan2(x[1], x[0]))) elif self.ambient_dim == 3: from grudge.shortcuts import make_visualizer self.vis = make_visualizer(discr, vis_order=order) @@ -63,7 +62,7 @@ class Plotter: return if self.ambient_dim == 2: - u = evt.state_component.get(self.queue) + u = self.actx.to_numpy(flatten(evt.state_component)) filename = "%s.png" % basename if not overwrite and os.path.exists(filename): @@ -90,9 +89,10 @@ class Plotter: # }}} -def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): +def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) # {{{ parameters @@ -143,11 +143,11 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): } from grudge import DGDiscretizationWithBoundaries - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + discr = DGDiscretizationWithBoundaries(actx, mesh, order=order, quad_tag_to_group_factory=quad_tag_to_group_factory) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) - logger.info("nnodes: %d", volume_discr.nnodes) + logger.info("ndofs: %d", volume_discr.ndofs) logger.info("nelements: %d", volume_discr.mesh.nelements) # }}} @@ -163,14 +163,14 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): quad_tag=product_tag) bound_op = bind(discr, op.sym_operator()) - u = bind(discr, f_gaussian(sym_x))(queue, t=0) + u = bind(discr, f_gaussian(sym_x))(actx, t=0) def rhs(t, u): - return bound_op(queue, t=t, u=u) + return bound_op(actx, t=t, u=u) # check velocity is tangential sym_normal = sym.surface_normal(dim, dim=dim - 1, dd=sym.DD_VOLUME).as_vector() - error = la.norm(bind(discr, c.dot(sym_normal))(queue).get(queue)) + error = bind(discr, sym.norm(2, c.dot(sym_normal)))(actx) logger.info("u_dot_n: %.5e", error) # }}} @@ -179,7 +179,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): # compute time step h_min = bind(discr, - sym.h_max_from_volume(discr.ambient_dim, dim=discr.dim))(queue) + sym.h_max_from_volume(discr.ambient_dim, dim=discr.dim))(actx) dt = dt_factor * h_min/order**2 nsteps = int(final_time // dt) + 1 dt = final_time/nsteps + 1.0e-15 @@ -189,10 +189,10 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - plot = Plotter(queue, discr, order, visualize=visualize) + plot = Plotter(actx, discr, order, visualize=visualize) norm = bind(discr, sym.norm(2, sym.var("u"))) - norm_u = 0.0 + norm_u = norm(actx, u=u) step = 0 @@ -203,19 +203,18 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): from grudge.shortcuts import make_visualizer vis = make_visualizer(discr, vis_order=order) vis.write_vtk_file("fld-surface-velocity.vtu", [ - ("u", bind(discr, c)(queue)), - ("n", bind(discr, sym_normal)(queue)) + ("u", bind(discr, c)(actx)), + ("n", bind(discr, sym_normal)(actx)) ], overwrite=True) df = sym.DOFDesc(sym.FACE_RESTR_INTERIOR) face_discr = discr.connection_from_dds(sym.DD_VOLUME, df).to_discr face_normal = bind(discr, sym.normal( - df, face_discr.ambient_dim, dim=face_discr.dim))(queue) + df, face_discr.ambient_dim, dim=face_discr.dim))(actx) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, face_discr, vis_order=order, - use_high_order_vtk=False) + vis = make_visualizer(actx, face_discr, vis_order=order) vis.write_vtk_file("fld-surface-face-normals.vtu", [ ("n", face_normal) ], overwrite=True) @@ -225,12 +224,14 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): continue step += 1 - if step % 1 == 0: - norm_u = norm(queue, u=event.state_component) + if step % 10 == 0: + norm_u = norm(actx, u=event.state_component) plot(event, "fld-surface-%04d" % step) logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + plot(event, "fld-surface-%04d" % step) + # }}} diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 388c4ce5..70235cee 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -28,6 +28,7 @@ import numpy as np import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flatten from grudge import bind, sym from pytools.obj_array import flat_obj_array @@ -39,8 +40,8 @@ logger = logging.getLogger(__name__) # {{{ plotting (keep in sync with `weak.py`) class Plotter: - def __init__(self, queue, discr, order, visualize=True, ylim=None): - self.queue = queue + def __init__(self, actx, discr, order, visualize=True, ylim=None): + self.actx = actx self.dim = discr.ambient_dim self.visualize = visualize @@ -53,7 +54,7 @@ class Plotter: self.ylim = ylim volume_discr = discr.discr_from_dd(sym.DD_VOLUME) - self.x = volume_discr.nodes().get(self.queue) + self.x = actx.to_numpy(flatten(thaw(actx, volume_discr.nodes()[0]))) else: from grudge.shortcuts import make_visualizer self.vis = make_visualizer(discr, vis_order=order) @@ -63,7 +64,7 @@ class Plotter: return if self.dim == 1: - u = evt.state_component.get(self.queue) + u = self.actx.to_numpy(flatten(evt.state_component)) filename = "%s.png" % basename if not overwrite and os.path.exists(filename): @@ -71,8 +72,8 @@ class Plotter: raise FileExistsError("output file '%s' already exists" % filename) ax = self.fig.gca() - ax.plot(self.x[0], u, '-') - ax.plot(self.x[0], u, 'k.') + ax.plot(self.x, u, '-') + ax.plot(self.x, u, 'k.') if self.ylim is not None: ax.set_ylim(self.ylim) @@ -188,7 +189,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - plot = Plotter(queue, discr, order, visualize=visualize, + plot = Plotter(actx, discr, order, visualize=visualize, ylim=[-0.1, 1.1]) step = 0 diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 88703a0e..b0692755 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -24,6 +24,7 @@ import numpy.linalg as la import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flatten from grudge import bind, sym @@ -34,8 +35,8 @@ logger = logging.getLogger(__name__) # {{{ plotting (keep in sync with `var-velocity.py`) class Plotter: - def __init__(self, queue, discr, order, visualize=True, ylim=None): - self.queue = queue + def __init__(self, actx, discr, order, visualize=True, ylim=None): + self.actx = actx self.dim = discr.ambient_dim self.visualize = visualize @@ -48,7 +49,7 @@ class Plotter: self.ylim = ylim volume_discr = discr.discr_from_dd(sym.DD_VOLUME) - self.x = volume_discr.nodes().get(self.queue) + self.x = actx.to_numpy(flatten(thaw(actx, volume_discr.nodes()[0]))) else: from grudge.shortcuts import make_visualizer self.vis = make_visualizer(discr, vis_order=order) @@ -58,7 +59,7 @@ class Plotter: return if self.dim == 1: - u = evt.state_component.get(self.queue) + u = self.actx.to_numpy(flatten(evt.state_component)) filename = "%s.png" % basename if not overwrite and os.path.exists(filename): @@ -66,8 +67,8 @@ class Plotter: raise FileExistsError("output file '%s' already exists" % filename) ax = self.fig.gca() - ax.plot(self.x[0], u, '-') - ax.plot(self.x[0], u, 'k.') + ax.plot(self.x, u, '-') + ax.plot(self.x, u, 'k.') if self.ylim is not None: ax.set_ylim(self.ylim) @@ -154,7 +155,7 @@ def main(ctx_factory, dim=2, order=4, visualize=True): from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - plot = Plotter(queue, discr, order, visualize=visualize, + plot = Plotter(actx, discr, order, visualize=visualize, ylim=[-1.1, 1.1]) norm = bind(discr, sym.norm(2, sym.var("u"))) -- GitLab From 9189cc2358caff196bd02225400f8cda0e5dcbe1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 14 Jul 2020 22:19:19 -0500 Subject: [PATCH 15/29] disable visualization in example --- examples/advection/surface.py | 5 ++--- examples/advection/var-velocity.py | 2 +- examples/advection/weak.py | 2 +- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 424e100a..9fcadb43 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -18,9 +18,8 @@ along with this program. If not, see . """ import os -import numpy as np -import numpy.linalg as la +import numpy as np import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext @@ -89,7 +88,7 @@ class Plotter: # }}} -def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): +def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 70235cee..cc2e62f0 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -90,7 +90,7 @@ class Plotter: # }}} -def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): +def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index b0692755..c2e36fb9 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -86,7 +86,7 @@ class Plotter: # }}} -def main(ctx_factory, dim=2, order=4, visualize=True): +def main(ctx_factory, dim=2, order=4, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) -- GitLab From 76deb09383247cc8404b8263d02d074c670fb031 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 20 Aug 2020 11:14:21 -0500 Subject: [PATCH 16/29] add new primitives to docs --- grudge/symbolic/primitives.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 3e3256c1..bc6c54d3 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -104,10 +104,18 @@ Geometry data .. autofunction:: inverse_metric_derivative .. autofunction:: forward_metric_derivative_mat .. autofunction:: inverse_metric_derivative_mat +.. autofunction:: first_fundamental_form +.. autofunction:: inverse_first_fundamental_form +.. autofunction:: inverse_surface_metric_derivative +.. autofunction:: second_fundamental_form +.. autofunction:: shape_operator .. autofunction:: pseudoscalar .. autofunction:: area_element .. autofunction:: mv_normal .. autofunction:: normal +.. autofunction:: surface_normal +.. autofunction:: summed_curvature +.. autofunction:: mean_curvature """ -- GitLab From cdde5edae65c952b3bc1eeaaeea9f17eb23466fc Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 19:28:27 -0500 Subject: [PATCH 17/29] update assert in test_surface_mass_operator_inverse --- test/test_grudge.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 096882aa..7778105e 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -322,8 +322,9 @@ def test_surface_mass_operator_inverse(ctx_factory, name): logger.info("inverse mass error\n%s", str(eoc)) - assert eoc.max_error() < 5.0e-09 \ - or eoc.order_estimate() > builder.order + # 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 # }}} -- GitLab From e60b3d859e09e51dacaf5ffd2bab74f0132a440c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 19:30:29 -0500 Subject: [PATCH 18/29] remove global basicConfig test_grudge --- test/test_grudge.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 7778105e..d138cc35 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -41,7 +41,6 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) -logging.basicConfig(level=logging.INFO) # {{{ inverse metric -- GitLab From 52160fcc241f9df85d7fabc6775e9c49e0229a26 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 19:35:08 -0500 Subject: [PATCH 19/29] update SurfaceAdvectionOperator --- grudge/models/advection.py | 22 +++------------------- grudge/symbolic/primitives.py | 4 ---- 2 files changed, 3 insertions(+), 23 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 659c10f8..b752a0dc 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -199,7 +199,9 @@ def v_dot_n_tpair(velocity, dd=None): def surface_advection_weak_flux(flux_type, u, velocity): v_dot_n = v_dot_n_tpair(velocity, dd=u.dd) - v_dot_n = sym.cse(0.5 * v_dot_n.jump, "v_dot_normal") + # NOTE: the normals in v_dot_n point to the exterior of their respective + # elements, so this is actually just an average + v_dot_n = sym.cse(0.5 * (v_dot_n.int - v_dot_n.ext), "v_dot_normal") flux_type = flux_type.lower() if flux_type == "central": @@ -210,16 +212,6 @@ def surface_advection_weak_flux(flux_type, u, velocity): raise ValueError("flux `{}` is not implemented".format(flux_type)) -def surface_penalty_flux(u, velocity, tau=1.0): - if abs(tau) < 1.0e-14: - return 0.0 - - v_dot_n = v_dot_n_tpair(velocity, dd=u.dd) - return sym.If(sym.Comparison(v_dot_n.avg, ">", 0), - 0.5 * tau * u.avg * v_dot_n.int, - 0.0) - - class SurfaceAdvectionOperator(AdvectionOperatorBase): def __init__(self, v, flux_type="central", quad_tag=None): super(SurfaceAdvectionOperator, self).__init__( @@ -230,19 +222,12 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase): surf_v = sym.project(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v) return surface_advection_weak_flux(self.flux_type, u, surf_v) - def penalty(self, u): - surf_v = sym.project(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v) - return surface_penalty_flux(u, surf_v, tau=0.0) - def sym_operator(self): u = sym.var("u") def flux(pair): return sym.project(pair.dd, face_dd)(self.flux(pair)) - def penalty(pair): - return sym.project(pair.dd, face_dd)(self.penalty(pair)) - face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag) @@ -258,7 +243,6 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase): for n in range(self.ambient_dim)) - sym.FaceMassOperator(face_dd, sym.DD_VOLUME)( flux(sym.int_tpair(u, self.quad_tag)) - + penalty(sym.int_tpair(u, self.quad_tag)) ) ) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index bc6c54d3..d06ae9f3 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -865,10 +865,6 @@ class TracePair: def ext(self): return self.exterior - @property - def jump(self): - return self.int - self.ext - @property def avg(self): return 0.5*(self.int + self.ext) -- GitLab From f8a69089b1a3f29e25897fc1ba53c11ea9569d0b Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 19:56:56 -0500 Subject: [PATCH 20/29] enable WADG if not all groups are affine --- grudge/execution.py | 4 +--- grudge/symbolic/mappers/__init__.py | 13 +++++-------- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index a651226f..59276ec0 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -755,9 +755,7 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=No )(sym_operator) dumper("before-global-to-reference", sym_operator) - sym_operator = mappers.GlobalToReferenceMapper( - discrwb.ambient_dim, - dim=discrwb.dim)(sym_operator) + sym_operator = mappers.GlobalToReferenceMapper(discrwb)(sym_operator) dumper("before-distributed", sym_operator) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index e8a63c1c..ebc771f9 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -587,18 +587,15 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): reference elements, together with explicit multiplication by geometric factors. """ - def __init__(self, ambient_dim, dim=None): + def __init__(self, discrwb): CSECachingMapperMixin.__init__(self) IdentityMapper.__init__(self) - if dim is None: - dim = ambient_dim + self.ambient_dim = discrwb.ambient_dim + self.dim = discrwb.dim - self.ambient_dim = ambient_dim - self.dim = dim - - # NOTE: only use WADG on surfaces at the moment - self.use_wadg = self.ambient_dim == (self.dim + 1) + volume_discr = discrwb.discr_from_dd(sym.DD_VOLUME) + self.use_wadg = not all(grp.is_affine for grp in volume_discr.groups) map_common_subexpression_uncached = \ IdentityMapper.map_common_subexpression -- GitLab From 225bdcba742a156eb6284b61446fb8962fffe3be Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 20:00:06 -0500 Subject: [PATCH 21/29] remove irrelevant TODO --- grudge/symbolic/mappers/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index ebc771f9..a6914e1f 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -643,7 +643,6 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): elif isinstance(expr.op, op.InverseMassOperator): 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)( -- GitLab From a6cb8d1abb4c72e02313de8022d75c85a8bcf2be Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 20:13:53 -0500 Subject: [PATCH 22/29] add some docs --- grudge/symbolic/primitives.py | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index d06ae9f3..d706314c 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -503,11 +503,24 @@ def mv_nodes(ambient_dim, dd=None): def forward_metric_nth_derivative(xyz_axis, ref_axes, dd=None): r""" - Pointwise metric derivatives representing + Pointwise metric derivatives representing repeated derivatives .. math:: - \frac{\partial x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{rst\_axis}} } + \frac{\partial^n x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{ref\_axes}}} + + where *ref_axes* is a multi-index description. + + :arg ref_axes: a :class:`tuple` of tuples indicating indices of + coordinate axes of the reference element to the number of derivatives + which will be taken. For example, the value ``((0, 2), (1, 1))`` + indicates taking the second derivative with respect to the first + axis and the first derivative with respect to the second + axis. Each axis must occur only once and the tuple must be sorted + by the axis index. + + May also be a singile integer *i*, which is viewed as equivalent + to ``((i, 1),)``. """ if isinstance(ref_axes, int): @@ -546,6 +559,14 @@ def forward_metric_nth_derivative(xyz_axis, ref_axes, dd=None): def forward_metric_derivative(xyz_axis, rst_axis, dd=None): + r""" + Pointwise metric derivatives representing + + .. math:: + + \frac{\partial x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{rst\_axis}}} + """ + return forward_metric_nth_derivative(xyz_axis, rst_axis, dd=dd) @@ -800,6 +821,8 @@ def normal(dd, ambient_dim, dim=None): def summed_curvature(ambient_dim, dim=None, dd=None): + """Sum of the principal curvatures""" + if dim is None: dim = ambient_dim - 1 @@ -814,6 +837,7 @@ def summed_curvature(ambient_dim, dim=None, dd=None): def mean_curvature(ambient_dim, dim=None, dd=None): + """Averaged (by dimension) sum of the principal curvatures.""" return 1.0 / (ambient_dim-1.0) * summed_curvature(ambient_dim, dim=dim, dd=dd) # }}} -- GitLab From 19a2a3e84bd04b75338a2df18a543aebacb22d11 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 25 Aug 2020 21:01:51 -0500 Subject: [PATCH 23/29] add a test for advection on curvilinear mesh --- test/test_grudge.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index d138cc35..25707b21 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -703,6 +703,7 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False): ("disk", [0.1, 0.05]), ("rect2", [4, 8]), ("rect3", [4, 6]), + ("warped2", [4, 8]), ]) @pytest.mark.parametrize("op_type", ["strong", "weak"]) @pytest.mark.parametrize("flux_type", ["central"]) @@ -746,11 +747,22 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type dim = 2 dt_factor = 4 elif mesh_name.startswith("rect"): - dim = int(mesh_name[4:]) + dim = int(mesh_name[-1:]) from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, n=(mesh_par,)*dim, order=4) + if dim == 2: + dt_factor = 4 + elif dim == 3: + dt_factor = 2 + else: + raise ValueError("dt_factor not known for %dd" % dim) + elif mesh_name.startswith("warped"): + dim = int(mesh_name[-1:]) + from meshmode.mesh.generation import generate_warped_rect_mesh + mesh = generate_warped_rect_mesh(dim, order=order, n=mesh_par) + if dim == 2: dt_factor = 4 elif dim == 3: @@ -831,7 +843,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.startswith("warped"): + # NOTE: curvilinear meshes are hard + assert eoc_rec.order_estimate() > order - 0.25 + else: + assert eoc_rec.order_estimate() > order # }}} -- GitLab From bdf0df5a56bf7aff6631dfbfa6a12ea1c2798477 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 26 Aug 2020 20:45:39 -0500 Subject: [PATCH 24/29] remove some useless face normal tests --- test/test_grudge.py | 37 ++----------------------------------- 1 file changed, 2 insertions(+), 35 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index faf4bc05..3ce0f8ae 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -324,14 +324,6 @@ def test_surface_mass_operator_inverse(actx_factory, name): # {{{ surface face normal orthogonality -def _avg_face_normal(x, side=-1): - x_i = x - x_e = sym.OppositeInteriorFaceSwap()(x_i) - x_avg = (x_i + side * x_e) / 2.0 - - return x_avg / sym.sqrt(x_avg.dot(x_avg)) - - @pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) def test_face_normal_surface(actx_factory, mesh_name): """Check that face normals are orthogonal to the surface normal""" @@ -370,14 +362,9 @@ def test_face_normal_surface(actx_factory, mesh_name): ) sym_surf_normal = sym_surf_normal / sym.sqrt(sum(sym_surf_normal**2)) - sym_surf_normal_avg = _avg_face_normal(sym_surf_normal, side=1.0) - sym_face_normal_i = sym.normal(df, ambient_dim, dim=dim - 1) sym_face_normal_e = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_i) - sym_face_normal_avg = _avg_face_normal(sym_face_normal_i) - sym_face_normal_op = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_avg) - if mesh.ambient_dim == 3: # NOTE: there's only one face tangent in 3d sym_face_tangent = ( @@ -394,46 +381,26 @@ def test_face_normal_surface(actx_factory, mesh_name): rtol = 1.0e-14 surf_normal = bind(discr, sym_surf_normal)(actx) - surf_normal_avg = bind(discr, sym_surf_normal_avg)(actx) face_normal_i = bind(discr, sym_face_normal_i)(actx) face_normal_e = bind(discr, sym_face_normal_e)(actx) - face_normal_avg = bind(discr, sym_face_normal_avg)(actx) - face_normal_op = bind(discr, sym_face_normal_op)(actx) - # check interpolated surface normal is orthogonal to face normal error = _eval_error(surf_normal.dot(face_normal_i)) logger.info("error[n_dot_i]: %.5e", error) assert error < rtol - # check averaged ones are also orthogonal - error = _eval_error(surf_normal_avg.dot(face_normal_avg)) - logger.info("error[a_dot_a]: %.5e", error) - assert error < rtol - - # check averaged face normal and interpolated face normal - error = _eval_error(surf_normal.dot(face_normal_avg)) - logger.info("error[n_dot_a]: %.5e", error) - assert error > rtol - # check angle between two neighboring elements error = _eval_error(face_normal_i.dot(face_normal_e) + 1.0) logger.info("error[i_dot_e]: %.5e", error) assert error > rtol - # check uniqueness of normal on the two sides - face_normal_sum = face_normal_avg + face_normal_op - error = _eval_error(face_normal_sum.dot(face_normal_sum)) - logger.info("error[a_plus_o]: %.5e", error) - assert error < rtol - # check orthogonality with face tangent if ambient_dim == 3: face_tangent = bind(discr, sym_face_tangent)(actx) - error = _eval_error(face_tangent.dot(face_normal_avg)) - logger.info("error[t_dot_avg]: %.5e", error) + error = _eval_error(face_tangent.dot(face_normal_i)) + logger.info("error[t_dot_i]: %.5e", error) assert error < 5 * rtol # }}} -- GitLab From e3c6cf5d2537599e0cce67ffc490ba06f18c84e8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 26 Aug 2020 20:48:40 -0500 Subject: [PATCH 25/29] fix typo in comment --- grudge/symbolic/primitives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 8d7d8132..5dab6139 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -799,7 +799,7 @@ def mv_normal(dd, ambient_dim, dim=None): # "normal plane") # # The trick done here is that we take the surface normal, move it to the - # face and then take a cross product with the face normal to get the + # face and then take a cross product with the face tangent to get the # correct exterior face normal vector. assert dim == ambient_dim - 2 -- GitLab From 7642da89e2809da77468418789237bc1e42d7163 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 26 Aug 2020 20:58:20 -0500 Subject: [PATCH 26/29] use only quad_tag_to_group_factory in surface advection example --- examples/advection/surface.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 9fcadb43..00ff58ac 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -134,15 +134,19 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): if product_tag == "none": product_tag = None + from meshmode.discretization.poly_element import \ + PolynomialWarpAndBlendGroupFactory, \ + QuadratureSimplexGroupFactory + + quad_tag_to_group_factory[sym.QTAG_NONE] = \ + PolynomialWarpAndBlendGroupFactory(order) + if product_tag: - from meshmode.discretization.poly_element import \ - QuadratureSimplexGroupFactory - quad_tag_to_group_factory = { - product_tag: QuadratureSimplexGroupFactory(order=4*order) - } + quad_tag_to_group_factory[product_tag] = \ + QuadratureSimplexGroupFactory(order=4*order) from grudge import DGDiscretizationWithBoundaries - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order, + discr = DGDiscretizationWithBoundaries(actx, mesh, quad_tag_to_group_factory=quad_tag_to_group_factory) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) -- GitLab From 8a6497259478c87d48c6dacf2e4001853a9f1e60 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 26 Aug 2020 21:02:09 -0500 Subject: [PATCH 27/29] rename initial condition in surface advection example --- examples/advection/surface.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 00ff58ac..f8580eef 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -157,7 +157,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): # {{{ symbolic operators - def f_gaussian(x): + def f_initial_condition(x): return x[0] from grudge.models.advection import SurfaceAdvectionOperator @@ -166,7 +166,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): quad_tag=product_tag) bound_op = bind(discr, op.sym_operator()) - u = bind(discr, f_gaussian(sym_x))(actx, t=0) + u0 = bind(discr, f_initial_condition(sym_x))(actx, t=0) def rhs(t, u): return bound_op(actx, t=t, u=u) @@ -191,15 +191,15 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): logger.info("nsteps: %d", nsteps) from grudge.shortcuts import set_up_rk4 - dt_stepper = set_up_rk4("u", dt, u, rhs) + dt_stepper = set_up_rk4("u", dt, u0, rhs) plot = Plotter(actx, discr, order, visualize=visualize) norm = bind(discr, sym.norm(2, sym.var("u"))) - norm_u = norm(actx, u=u) + norm_u = norm(actx, u=u0) step = 0 - event = dt_stepper.StateComputed(0.0, 0, 0, u) + event = dt_stepper.StateComputed(0.0, 0, 0, u0) plot(event, "fld-surface-%04d" % 0) if visualize and dim == 3: -- GitLab From 3c1aef428c8a57793f8f2d8c67b3db1b2bf9fec9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 4 Sep 2020 18:55:16 -0500 Subject: [PATCH 28/29] update example license --- examples/advection/surface.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index f8580eef..2d9c8d12 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -1,20 +1,23 @@ -from __future__ import division, absolute_import - __copyright__ = "Copyright (C) 2020 Alexandru Fikl" __license__ = """ -This program is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -This program is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with this program. If not, see . +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 os -- GitLab From 2ab2b1f235646fc01540cea96c55e703b711dfa6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 4 Sep 2020 19:29:17 -0500 Subject: [PATCH 29/29] reduce time to make surface advection example faster --- examples/advection/surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 2d9c8d12..c0bb772c 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -104,9 +104,9 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): resolution = 64 if dim == 2 else 1 # cfl - dt_factor = 1.0 + dt_factor = 2.0 # final time - final_time = 2.0 * np.pi + final_time = np.pi # velocity field sym_x = sym.nodes(dim) -- GitLab