diff --git a/examples/advection/surface.py b/examples/advection/surface.py new file mode 100644 index 0000000000000000000000000000000000000000..c0bb772cd4c9b3c811530b2835466c079ef97184 --- /dev/null +++ b/examples/advection/surface.py @@ -0,0 +1,255 @@ +__copyright__ = "Copyright (C) 2020 Alexandru Fikl" + +__license__ = """ +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 + +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 make_obj_array + +import logging +logger = logging.getLogger(__name__) + + +# {{{ plotting (keep in sync with `var-velocity.py`) + +class Plotter: + def __init__(self, actx, discr, order, visualize=True): + self.actx = actx + 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) + + 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) + else: + raise ValueError("unsupported dimension") + + def __call__(self, evt, basename, overwrite=True): + if not self.visualize: + return + + if self.ambient_dim == 2: + u = self.actx.to_numpy(flatten(evt.state_component)) + + 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=False): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + # {{{ parameters + + # sphere radius + radius = 1.0 + # sphere resolution + resolution = 64 if dim == 2 else 1 + + # cfl + dt_factor = 2.0 + # final time + final_time = np.pi + + # velocity field + sym_x = sym.nodes(dim) + c = make_obj_array([ + -sym_x[1], sym_x[0], 0.0 + ])[:dim] + # flux + flux_type = "lf" + + # }}} + + # {{{ 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 + mesh = generate_icosphere(radius, order=4 * order, + uniform_refinement_rounds=resolution) + else: + raise ValueError("unsupported dimension") + + quad_tag_to_group_factory = {} + 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: + quad_tag_to_group_factory[product_tag] = \ + QuadratureSimplexGroupFactory(order=4*order) + + from grudge import DGDiscretizationWithBoundaries + discr = DGDiscretizationWithBoundaries(actx, mesh, + quad_tag_to_group_factory=quad_tag_to_group_factory) + + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + logger.info("ndofs: %d", volume_discr.ndofs) + logger.info("nelements: %d", volume_discr.mesh.nelements) + + # }}} + + # {{{ symbolic operators + + def f_initial_condition(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()) + u0 = bind(discr, f_initial_condition(sym_x))(actx, t=0) + + def rhs(t, 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 = bind(discr, sym.norm(2, c.dot(sym_normal)))(actx) + logger.info("u_dot_n: %.5e", error) + + # }}} + + # {{{ time stepping + + # compute time step + h_min = bind(discr, + 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 + + logger.info("dt: %.5e", dt) + logger.info("nsteps: %d", nsteps) + + from grudge.shortcuts import set_up_rk4 + 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=u0) + + step = 0 + + event = dt_stepper.StateComputed(0.0, 0, 0, u0) + 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)(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))(actx) + + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(actx, face_discr, vis_order=order) + 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 % 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) + + # }}} + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser() + 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) + main(cl.create_some_context, + dim=args.dim, + product_tag=args.qtag) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 388c4ce576fd8e248432831f14efaadb248fb4a1..cc2e62f06125d7b4f16b4c291911d9184f6fafd0 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) @@ -89,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) @@ -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 05983d3cbb1e19f71a5fd7848edc016a4cf187b4..77ce17bc54ed03e84da4f5761d0c251abeb52824 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -29,6 +29,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 @@ -39,8 +40,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 @@ -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) @@ -90,7 +91,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) @@ -159,7 +160,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"))) diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py index fd37f27f867bd9db0bd195371cd768a80a243f25..bab82d545df9a7ba9bf4f6b41937dc10f45191df 100755 --- a/examples/dagrt-fusion.py +++ b/examples/dagrt-fusion.py @@ -611,6 +611,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper): # TODO: Not comprehensive. sym_op.ProjectionOperator, 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 0f6711e8a1467a2a21ac65246d552a03f51afc20..59276ec05aabfcdc0eaee55f93749be07318086f 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -74,7 +74,7 @@ class ExecutionMapper(mappers.Evaluator, result = discr.empty(self.array_context) for grp_ary in result: - grp_ary.fill(1) + grp_ary.fill(1.0) return result def map_node_coordinate_component(self, expr): @@ -755,7 +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)(sym_operator) + sym_operator = mappers.GlobalToReferenceMapper(discrwb)(sym_operator) dumper("before-distributed", sym_operator) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index e160643d5232ca3e6d7c196afe39a6e96b9cf410..b752a0dccea4a9dc0c6036a6432786f16dd32ec1 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -182,4 +182,70 @@ class VariableCoefficientAdvectionOperator(AdvectionOperatorBase): )) # }}} + +# {{{ closed surface advection + +def v_dot_n_tpair(velocity, dd=None): + if dd is None: + dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR) + + ambient_dim = len(velocity) + normal = sym.normal(dd.with_qtag(None), ambient_dim, dim=ambient_dim - 2) + + 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): + v_dot_n = v_dot_n_tpair(velocity, dd=u.dd) + # 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": + return u.avg * v_dot_n + elif flux_type == "lf": + 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)) + + +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.project(sym.DD_VOLUME, u.dd.with_qtag(None))(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.project(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.project(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 diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index b5698c0e750511131180ce60f6177d27a3801998..e5b7f7e872541aacd7e352eb6e247e99915077c4 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -918,7 +918,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 6d060e4691256404e12b600f9dfa6236f1216e7d..a6914e1fd56af9cb43210177a05e2006aed45d34 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -587,15 +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 + 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 @@ -605,44 +605,63 @@ 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) + 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(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: + # 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 +669,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 +824,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 +1196,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/operators.py b/grudge/symbolic/operators.py index df886d4571379570305e46087836d16c88ba28c3..6a5a982ce29cf3a1168b1c3851f1c9a7195f59ec 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -734,7 +734,7 @@ def norm(p, arg, dd=None): return prim.sqrt(norm_squared) - elif p == np.Inf: + elif p == np.inf: result = NodalMax(dd_in=dd)(prim.fabs(arg)) if isinstance(result, np.ndarray): @@ -770,6 +770,27 @@ 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 diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 5376865883c3e397a152f432037d9ac963adbdff..5dab6139b5e360589f10034548b23d5537d8bd0d 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -106,10 +106,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 """ @@ -495,32 +503,73 @@ 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 + 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): + 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 project result = project(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): + 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) def forward_metric_derivative_vector(ambient_dim, rst_axis, dd=None): @@ -539,7 +588,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( @@ -593,6 +644,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 @@ -609,6 +661,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 @@ -627,14 +764,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 @@ -645,28 +776,72 @@ 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, - # project 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 project - project = project(DD_VOLUME, dd) - mv = MultiVector(np.array([ - mv.as_scalar() * project(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 - return cse(mv, "normal", cse_scope.DISCRETIZATION) + if dim == ambient_dim - 1: + return surface_normal(ambient_dim, dim=dim, dd=dd) + + # 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 tangent to get the + # correct exterior face normal vector. + assert dim == ambient_dim - 2 + + from grudge.symbolic.operators import project + volm_normal = ( + surface_normal(ambient_dim, dim=dim + 1, dd=DD_VOLUME) + .map(project(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): return mv_normal(dd, ambient_dim, dim).as_vector() + +def summed_curvature(ambient_dim, dim=None, dd=None): + """Sum of the principal curvatures""" + + 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): + """Averaged (by dimension) sum of the principal curvatures.""" + return 1.0 / (ambient_dim-1.0) * summed_curvature(ambient_dim, dim=dim, dd=dd) + # }}} @@ -721,20 +896,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 project, OppositeInteriorFaceSwap - i = project("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(project("int_faces", q_dd)(i)) - e = cse(project("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 = project(from_dd, trace_dd.with_qtag(None))(expression) + e = cse(OppositeInteriorFaceSwap()(i)) + + if trace_dd.uses_quadrature(): + i = cse(project(trace_dd.with_qtag(None), trace_dd)(i)) + e = cse(project(trace_dd.with_qtag(None), trace_dd)(e)) - return TracePair(q_dd, interior=i, exterior=e) + return TracePair(trace_dd, interior=i, exterior=e) def bdry_tpair(dd, interior, exterior): diff --git a/test/mesh_data.py b/test/mesh_data.py new file mode 100644 index 0000000000000000000000000000000000000000..2ea2cc4782e8157485ed60991e118ab7363d2172 --- /dev/null +++ b/test/mesh_data.py @@ -0,0 +1,132 @@ +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 = 4 + resolutions = [1.0, 0.11, 0.05] + # resolutions = [1.0, 0.11, 0.05, 0.03, 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=[ + "-optimize_ho", + "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution + ], + target_unit="MM") + + 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 25fdbc409227ba01135cb97efb05f4c80945ee42..3ce0f8ae89f00cd6d60fa660f79b76be0d895809 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -22,31 +22,28 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ - import numpy as np import numpy.linalg as la -import pytest # noqa - -from pytools.obj_array import flat_obj_array, make_obj_array import pyopencl as cl - from meshmode.array_context import PyOpenCLArrayContext from meshmode.dof_array import unflatten, flatten, flat_norm +from pytools.obj_array import flat_obj_array, make_obj_array + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import pytest from meshmode.array_context import ( # noqa pytest_generate_tests_for_pyopencl_array_context as pytest_generate_tests) - import logging logger = logging.getLogger(__name__) -logging.basicConfig() -logger.setLevel(logging.INFO) -from grudge import sym, bind, DGDiscretizationWithBoundaries +# {{{ inverse metric @pytest.mark.parametrize("dim", [2, 3]) def test_inverse_metric(actx_factory, dim): @@ -91,6 +88,10 @@ def test_inverse_metric(actx_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"]) @@ -156,6 +157,258 @@ def test_mass_mat_trig(actx_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 + 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(actx_factory, name): + actx = actx_factory() + + # {{{ 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(actx, mesh, order=builder.order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + + logger.info("ndofs: %d", volume_discr.ndofs) + logger.info("nelements: %d", volume_discr.mesh.nelements) + + # {{{ 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)(actx) + + 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, dim=discr.dim, dd=dd))(actx) + 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 + +# }}} + + +# {{{ surface mass inverse + +@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) +def test_surface_mass_operator_inverse(actx_factory, name): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + # {{{ cases + + if name == "2-1-ellipse": + from mesh_data import EllipseMeshBuilder + builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + elif name == "spheroid": + from mesh_data import SpheroidMeshBuilder + builder = SpheroidMeshBuilder() + 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(actx, mesh, order=builder.order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + + logger.info("ndofs: %d", volume_discr.ndofs) + logger.info("nelements: %d", volume_discr.mesh.nelements) + + # {{{ 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)(actx) + f_inv = bind(discr, sym_op)(actx, f=f) + + inv_error = bind(discr, + sym.norm(2, sym.var("x") - sym.var("y")) + / sym.norm(2, sym.var("y")))(actx, x=f_inv, y=f) + + # }}} + + h_max = bind(discr, sym.h_max_from_volume( + discr.ambient_dim, dim=discr.dim, dd=dd))(actx) + eoc.add_data_point(h_max, inv_error) + + # }}} + + logger.info("inverse mass error\n%s", str(eoc)) + + # 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 + +# }}} + + +# {{{ surface face normal orthogonality + +@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""" + actx = actx_factory() + + # {{{ 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[0], builder.mesh_order) + discr = DGDiscretizationWithBoundaries(actx, mesh, order=builder.order) + + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + logger.info("ndofs: %d", volume_discr.ndofs) + logger.info("nelements: %d", volume_discr.mesh.nelements) + + # }}} + + # {{{ symbolic + + dv = sym.DD_VOLUME + df = sym.as_dofdesc(sym.FACE_RESTR_INTERIOR) + + ambient_dim = mesh.ambient_dim + dim = mesh.dim + + sym_surf_normal = sym.project(dv, df)( + 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=dim - 1) + sym_face_normal_e = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_i) + + if mesh.ambient_dim == 3: + # NOTE: there's only one face tangent in 3d + sym_face_tangent = ( + sym.pseudoscalar(ambient_dim, dim - 1, dd=df) + / sym.area_element(ambient_dim, dim - 1, dd=df)).as_vector() + + # }}} + + # {{{ checks + + def _eval_error(x): + return bind(discr, sym.norm(np.inf, sym.var("x", dd=df), dd=df))(actx, x=x) + + rtol = 1.0e-14 + + surf_normal = bind(discr, sym_surf_normal)(actx) + + face_normal_i = bind(discr, sym_face_normal_i)(actx) + face_normal_e = bind(discr, sym_face_normal_e)(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 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 orthogonality with face tangent + if ambient_dim == 3: + face_tangent = bind(discr, sym_face_tangent)(actx) + + error = _eval_error(face_tangent.dot(face_normal_i)) + logger.info("error[t_dot_i]: %.5e", error) + assert error < 5 * rtol + + # }}} + +# }}} + + +# {{{ diff operator @pytest.mark.parametrize("dim", [1, 2, 3]) def test_tri_diff_mat(actx_factory, dim, order=4): @@ -171,7 +424,7 @@ def test_tri_diff_mat(actx_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) @@ -193,8 +446,12 @@ def test_tri_diff_mat(actx_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 + +# }}} + +# {{{ divergence theorem def test_2d_gauss_theorem(actx_factory): """Verify Gauss's theorem explicitly on a mesh""" @@ -237,11 +494,167 @@ def test_2d_gauss_theorem(actx_factory): assert abs(gauss_err) < 1e-13 +@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"]) +def test_surface_divergence_theorem(actx_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). + """ + actx = actx_factory() + + # {{{ 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 flat_obj_array( + 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(actx, mesh, order=builder.order, + quad_tag_to_group_factory={ + "product": QuadratureSimplexGroupFactory(2 * builder.order) + }) + + volume = discr.discr_from_dd(sym.DD_VOLUME) + logger.info("ndofs: %d", volume.ndofs) + 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.project(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 + 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)) + + # evaluate + op_global = bind(discr, sym_op_global)(actx) + op_local = bind(discr, sym_op_local)(actx) + + err_global = abs(op_global) + err_local = bind(discr, sym.norm(np.inf, sym.var("x")))(actx, x=op_local) + 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, dim=discr.dim, dd=dd))(actx) + eoc_global.add_data_point(h_max, err_global) + eoc_local.add_data_point(h_max, err_local) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=builder.order) + + filename = f"surface_divergence_theorem_{mesh_name}_{i:04d}.vtu" + vis.write_vtk_file(filename, [ + ("r", actx.np.log10(op_local)) + ], overwrite=True) + + # }}} + + 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 - 0.5 + + assert eoc_local.max_error() < 1.0e-12 \ + or eoc_local.order_estimate() > order - 0.5 + +# }}} + + +# {{{ models: advection + @pytest.mark.parametrize(("mesh_name", "mesh_pars"), [ ("segment", [8, 16, 32]), ("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"]) @@ -283,7 +696,7 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ 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) @@ -294,7 +707,17 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ 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: + dt_factor = 2 + else: + raise ValueError("dt_factor not known for %dd" % dim) else: raise ValueError("invalid mesh name: " + mesh_name) @@ -369,8 +792,16 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ 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 + +# }}} + +# {{{ models: maxwell @pytest.mark.parametrize("order", [3, 4, 5]) def test_convergence_maxwell(actx_factory, order): @@ -440,6 +871,10 @@ def test_convergence_maxwell(actx_factory, order): assert eoc_rec.order_estimate() > order +# }}} + + +# {{{ models: variable coefficient advection oversampling @pytest.mark.parametrize("order", [2, 3, 4]) def test_improvement_quadrature(actx_factory, order): @@ -508,6 +943,10 @@ def test_improvement_quadrature(actx_factory, order): assert (q_errs < errs).all() assert q_eoc > order +# }}} + + +# {{{ operator collector determinism def test_op_collector_order_determinism(): class TestOperator(sym.Operator): @@ -533,6 +972,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(actx_factory): actx = actx_factory() @@ -561,6 +1004,10 @@ def test_bessel(actx_factory): assert z < 1e-15 +# }}} + + +# {{{ function symbol def test_external_call(actx_factory): actx = actx_factory() @@ -630,6 +1077,8 @@ def test_function_symbol_array(ctx_factory, array_type): norm = bind(discr, sym.norm(2, sym_x))(x=x) assert isinstance(norm, float) +# }}} + @pytest.mark.parametrize("p", [2, np.inf]) def test_norm_obj_array(ctx_factory, p): @@ -699,7 +1148,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