From 7ccdba2729038aeb59739dad1f99df9c56e5645a Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Thu, 23 Apr 2020 21:16:58 -0500
Subject: [PATCH] another attempt at making 1D work nicely
The problem is that in 1D faces are 0D and so carry exactly no
useful information about which side of the element they may be
on.
To hack around that, this commit introduces a hacky operator
that hacks it. Less ambiguously, we just use the connection to
figure out which face we're on on set the normal to `+1` or `-1`
based on that
---
examples/advection/weak.py | 2 --
examples/dagrt-fusion.py | 2 +-
grudge/execution.py | 15 +++++++++
grudge/symbolic/dofdesc_inference.py | 1 +
grudge/symbolic/mappers/__init__.py | 3 ++
grudge/symbolic/operators.py | 7 -----
grudge/symbolic/primitives.py | 46 +++++++++++++++++++++-------
7 files changed, 55 insertions(+), 21 deletions(-)
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 85545303..6e30094d 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -67,9 +67,7 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
from grudge import DGDiscretizationWithBoundaries
discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
-
volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
- faces_discr = discr.discr_from_dd(sym.FACE_RESTR_INTERIOR)
# }}}
diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index 27e35abd..bfb6f761 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -159,7 +159,7 @@ def transcribe_phase(dag, field_var_name, field_components, phase_name,
"<dt>": sym.var("input_dt", sym.DD_SCALAR),
f"<state>{field_var_name}": sym.make_sym_array(
f"input_{field_var_name}", field_components),
- f"<p>residual": sym.make_sym_array(
+ "<p>residual": sym.make_sym_array(
"input_residual", field_components),
}
diff --git a/grudge/execution.py b/grudge/execution.py
index 948be131..fade7fca 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -283,6 +283,21 @@ class ExecutionMapper(mappers.Evaluator,
# }}}
+ def map_signed_face_ones(self, expr):
+ face_discr = self.discrwb.discr_from_dd(expr.dd)
+ assert face_discr.dim == 0
+
+ all_faces_conn = self.discrwb.connection_from_dds("vol", expr.dd)
+ field = face_discr.empty(self.queue, allocator=self.bound_op.allocator)
+ field.fill(1)
+
+ for grp in all_faces_conn.groups:
+ for batch in grp.batches:
+ i = batch.to_element_indices.with_queue(self.queue)
+ field[i] = (2.0 * (batch.to_element_face % 2) - 1.0) * field[i]
+
+ return field
+
# }}}
# {{{ instruction execution functions
diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py
index c44a7940..38c77512 100644
--- a/grudge/symbolic/dofdesc_inference.py
+++ b/grudge/symbolic/dofdesc_inference.py
@@ -183,6 +183,7 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin):
return expr.dd
map_node_coordinate_component = map_ones
+ map_signed_face_ones = map_ones
def map_call(self, expr):
arg_dds = [
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index 0e9b8e02..b277c4c9 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -213,6 +213,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
map_function_symbol = map_grudge_variable
map_ones = map_grudge_variable
+ map_signed_face_ones = map_grudge_variable
map_node_coordinate_component = map_grudge_variable
# }}}
@@ -266,6 +267,7 @@ class DependencyMapper(
return set()
map_ones = _map_leaf
+ map_signed_face_ones = _map_leaf
map_node_coordinate_component = _map_leaf
@@ -1235,6 +1237,7 @@ class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMix
map_function_symbol = map_constant
map_ones = map_grudge_variable
+ map_signed_face_ones = map_grudge_variable
map_node_coordinate_component = map_grudge_variable
map_operator = map_grudge_variable
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 0cbb6eda..f1997b1c 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -536,13 +536,6 @@ class RefFaceMassOperator(ElementwiseLinearOperator):
volgrp.order,
face_vertices)
- if afgrp.dim == 0 and volgrp.dim == 1:
- # NOTE: This complements a choice in `parametrization_derivative`
- # where we don't really get a sign for the 0-dim case, so we add
- # signs here in the hope that it'll only be used for fluxes where
- # this makes sense
- matrix[0, 0, 0] = -matrix[0, 0, 0]
-
# np.set_printoptions(linewidth=200, precision=3)
# matrix[np.abs(matrix) < 1e-13] = 0
# print(matrix)
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 80dd8413..875a7bff 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -426,6 +426,20 @@ class Ones(HasDOFDesc, ExpressionBase):
mapper_method = intern("map_ones")
+class _SignedFaceOnes(HasDOFDesc, ExpressionBase):
+ """Sets the face to :math:`-1` if its odd and :math:`+1` if its even.
+ This is based on the face order of
+ :meth:`meshmode.mesh.MeshElementGroup.face_vertex_indices`.
+ """
+
+ def __init__(self, dd):
+ dd = as_dofdesc(dd)
+ assert dd.is_trace()
+ super(_SignedFaceOnes, self).__init__(dd)
+
+ mapper_method = intern("map_signed_face_ones")
+
+
# {{{ geometry data
class DiscretizationProperty(HasDOFDesc, ExpressionBase):
@@ -507,8 +521,7 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None):
dim = ambient_dim
if dim == 0:
- inner_dd = (DD_VOLUME if dd is None else dd).with_qtag(QTAG_NONE)
- return MultiVector(np.array([Ones(dd=inner_dd)]))
+ return MultiVector(np.array([_SignedFaceOnes(dd)]))
from pytools import product
return product(
@@ -607,18 +620,29 @@ def mv_normal(dd, ambient_dim, dim=None):
if dim is None:
dim = ambient_dim - 1
- # Don't be tempted to add a sign here. As it is, it produces
+ # NOTE: Don't be tempted to add a sign here. As it is, it produces
# exterior normals for positively oriented curves.
- pder = (
- pseudoscalar(ambient_dim, dim, dd=dd)
- / # noqa: W504
- area_element(ambient_dim, dim, dd=dd))
+ pder = pseudoscalar(ambient_dim, dim, dd=dd) \
+ / area_element(ambient_dim, dim, dd=dd)
- return cse(
- # Dorst Section 3.7.2
- pder if dim == 0 else pder << pder.I.inv(),
- "normal", cse_scope.DISCRETIZATION)
+ # Dorst Section 3.7.2
+ mv = pder << pder.I.inv()
+
+ 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())
+
+ interp = _sym().interp(DD_VOLUME, dd)
+ mv = MultiVector(np.array([
+ mv.as_scalar() * interp(t) for t in tangent.as_vector()
+ ]))
+
+ return cse(mv, "normal", cse_scope.DISCRETIZATION)
def normal(dd, ambient_dim, dim=None, quadrature_tag=None):
--
GitLab