diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 04b4c635ef61a5a1ad96d36b5684812e5b34d4f9..2e0e438352cbdb7930678602f5734d4cc699c9af 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_metric_derivative( + 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)(imd(rst_axis) * rec_field) - 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))) @@ -1197,7 +1201,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 8948c31718126687cf1b6584eba0ab2a8976a01e..3ba44d0803f3b1ed838db4b05916f55dd5598a4f 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -489,7 +489,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 @@ -497,24 +497,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): @@ -605,6 +625,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, 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 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 @@ -677,6 +782,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 d53b104011ec269791d32bddab2492ab3c428a52..5576e99536352db9d8a0c87c1f9aa9bc353143ef 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 c44909803d15f3990993ab1837e6b260a119af1e..d9d5589110f0d0788b60b8a777743b1708ccdc26 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -454,6 +454,88 @@ def test_2d_gauss_theorem(ctx_factory): assert abs(gauss_err) < 1e-13 +@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) +@pytest.mark.parametrize("op_type", ["stiffness", "diff"]) +def test_surface_divergence_theorem(ctx_factory, mesh_name, op_type): + """ + The surface divergence theorem + """ + 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 == "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[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() + + 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 + + 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) + + # }}} + + logger.info("\n%s", str(eoc)) + + 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]),