From cd1f91fd4e3dc8679484086bb3ceb68f386105fe Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 18 May 2020 13:46:00 -0500 Subject: [PATCH 1/3] add surface tangential derivative --- grudge/symbolic/mappers/__init__.py | 2 +- grudge/symbolic/primitives.py | 51 +++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 1 deletion(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index b6fcc0fe..713ef856 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -673,7 +673,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): rec_field = jac_tag * rec_field def imd(rst): - return sym.inverse_metric_derivative( + return sym.inverse_surface_metric_derivative( rst, expr.op.xyz_axis, ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index d7584539..fd1fdbe3 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -598,6 +598,57 @@ def inverse_metric_derivative_mat(ambient_dim, dim=None, dd=None): return result +def first_fundamental_form(ambient_dim, dim=None, dd=None): + 1/0 + 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, rst] + * forward_metric_derivative(xyz_axis, rst, dd=dd) + for rst in range(dim)) + + return cse(imd, + prefix="ds%d_dx%d" % (rst_axis, xyz_axis), + scope=cse_scope.DISCRETIZATION) + + def pseudoscalar(ambient_dim, dim=None, dd=None): if dim is None: dim = ambient_dim -- GitLab From 208de0bde46ba1736d89d49791177331412d101f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 19 May 2020 10:46:33 -0500 Subject: [PATCH 2/3] tests surface derivatives with surface divergence theorem --- grudge/symbolic/mappers/__init__.py | 4 +- grudge/symbolic/primitives.py | 92 +++++++++++++++++++++++++---- test/mesh_data.py | 8 +-- test/test_grudge.py | 75 +++++++++++++++++++++++ 4 files changed, 163 insertions(+), 16 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 713ef856..d2868f7c 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -679,7 +679,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): dd=dd_in) return sum( - ref_class(rst_axis, dd_in=dd_in)(imd(rst_axis) * rec_field) + 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): @@ -1243,7 +1243,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/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index fd1fdbe3..d1f5d39e 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -482,7 +482,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 @@ -490,24 +490,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): @@ -599,7 +619,6 @@ def inverse_metric_derivative_mat(ambient_dim, dim=None, dd=None): def first_fundamental_form(ambient_dim, dim=None, dd=None): - 1/0 if dim is None: dim = ambient_dim - 1 @@ -649,6 +668,41 @@ def inverse_surface_metric_derivative(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 @@ -721,6 +775,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), "mean_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 d53b1040..5576e995 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -4,14 +4,14 @@ 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) - @property - def mesh_order(self): - return self.order + if self.mesh_order is None: + self.mesh_order = self.order @property def ambient_dim(self): @@ -60,7 +60,7 @@ class SphereMeshBuilder(MeshBuilder): class SpheroidMeshBuilder(MeshBuilder): ambient_dim = 3 - order = 2 + mesh_order = 2 resolutions = [1.0, 0.1, 0.05] @property diff --git a/test/test_grudge.py b/test/test_grudge.py index caf7f3f2..fdd224fc 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -454,6 +454,81 @@ 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): + """ + The surface divergence theorem + """ + cl_ctx = cl.create_some_context() + 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 == "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 div(dim, arg): + return sum(nabla(u) for nabla, u in zip(sym.nabla(dim), arg)) + + def f(x): + return join_fields( + sym.sin(3*x[0]) + sym.cos(3*x[1]), + sym.sin(2*x[0]) + sym.cos(x[1]), + 0.0)[:ambient_dim] + + 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) + logger.info("nnodes: %d", discr._volume_discr.nnodes) + + dd = sym.DD_VOLUME + ambient_dim = discr.ambient_dim + dim = discr.dim + + sym_f = f(sym.nodes(ambient_dim, dd=dd)) + sym_kappa = sym.summed_curvature(ambient_dim, dim=dim, dd=dd) + sym_surf_normal = sym.surface_normal(ambient_dim, dim=dim, dd=dd).as_vector() + + sym_div = sym.integral(div(ambient_dim, sym_f), dd=dd) + sym_k = sym.integral(sym_kappa * sym_f.dot(sym_surf_normal), dd=dd) + sym_op = sym_div - sym_k + + value = la.norm(bind(discr, sym_op)(queue)) + logger.info("integral %.5e", value) + + # compute max element size + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc.add_data_point(h_max, value) + + # }}} + + from pystopt.tools import stringify_eoc + logger.info("\n%s", stringify_eoc(eoc, table_format="latex")) + + order = min(builder.order, builder.mesh_order) - 0.5 + assert eoc.max_error() < 1.0e-12 \ + or eoc.order_estimate() >= order + + @pytest.mark.parametrize(("mesh_name", "mesh_pars"), [ ("segment", [8, 16, 32]), ("disk", [0.1, 0.05]), -- GitLab From b12d72838e5d14ab3b3b61172d822d8d8ab3be02 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 20 May 2020 11:00:56 -0500 Subject: [PATCH 3/3] fix stiffness_t operator --- grudge/symbolic/mappers/__init__.py | 22 +++++++++++++--------- test/test_grudge.py | 23 +++++++++++++++-------- 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 861465e5..2e0e4383 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -622,20 +622,24 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): dd=dd_in.with_qtag(sym.QTAG_NONE)) def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True): - rec_field = self.rec(field) - if with_jacobian: - jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in) - rec_field = jac_tag * rec_field - def imd(rst): return sym.inverse_surface_metric_derivative( rst, expr.op.xyz_axis, ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in) - return sum( - ref_class(rst_axis, dd_in=dd_in)(rec_field) * imd(rst_axis) - for rst_axis in range(self.dim)) + 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( + 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)( @@ -662,7 +666,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): 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))) diff --git a/test/test_grudge.py b/test/test_grudge.py index cbd4d5cb..d9d55891 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -455,11 +455,12 @@ def test_2d_gauss_theorem(ctx_factory): @pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) -def test_surface_divergence_theorem(ctx_factory, mesh_name): +@pytest.mark.parametrize("op_type", ["stiffness", "diff"]) +def test_surface_divergence_theorem(ctx_factory, mesh_name, op_type): """ The surface divergence theorem """ - cl_ctx = cl.create_some_context() + cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) # {{{ cases @@ -483,9 +484,6 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name): # {{{ convergene - def div(dim, arg): - return sum(nabla(u) for nabla, u in zip(sym.nabla(dim), arg)) - def f(x): return join_fields( sym.sin(3*x[0]) + sym.cos(3*x[1]), @@ -508,7 +506,17 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name): sym_kappa = sym.summed_curvature(ambient_dim, dim=dim, dd=dd) sym_surf_normal = sym.surface_normal(ambient_dim, dim=dim, dd=dd).as_vector() - sym_div = sym.integral(div(ambient_dim, sym_f), dd=dd) + if op_type == "diff": + sym_div = sym.integral(sum( + nabla(f) for nabla, f in zip(sym.nabla(ambient_dim), sym_f)), + dd=dd) + elif op_type == "stiffness": + sym_div = sym.NodalSum(dd)(sum( + op(f) for op, f in zip(sym.stiffness(ambient_dim), sym_f) + )) + else: + raise ValueError("unknown op type '{}'".format(op_type)) + sym_k = sym.integral(sym_kappa * sym_f.dot(sym_surf_normal), dd=dd) sym_op = sym_div - sym_k @@ -521,8 +529,7 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name): # }}} - from pystopt.tools import stringify_eoc - logger.info("\n%s", stringify_eoc(eoc, table_format="latex")) + logger.info("\n%s", str(eoc)) order = min(builder.order, builder.mesh_order) - 0.5 assert eoc.max_error() < 1.0e-12 \ -- GitLab