diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index ddc7535288cd82c1bbc56cb8decd9b0644cea73a..00e4f0c8496dced20a21e808e68986ecec1d5be2 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -1,17 +1,26 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# 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 <http://www.gnu.org/licenses/>.
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2017 Bogdan Enache"
+
+__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 numpy as np
@@ -37,20 +46,19 @@ def main(write_output=True, order=4):
dim = 2
+ resolution = 15
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
- n=(15, 15), order=order)
+ n=(resolution, resolution), order=order)
- dt_factor = 4
- h = 1/20
-
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+ dt_factor = 5
+ h = 1/resolution
sym_x = sym.nodes(2)
advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
- flux_type = "central"
+ flux_type = "upwind"
source_center = np.array([0.1, 0.1])
source_width = 0.05
@@ -58,30 +66,46 @@ def main(write_output=True, order=4):
sym_x = sym.nodes(2)
sym_source_center_dist = sym_x - source_center
- def f(x):
+ def f_gaussian(x):
return sym.exp(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2)
+ def f_step(x):
+ return sym.If(
+ sym.Comparison(
+ np.dot(sym_source_center_dist, sym_source_center_dist),
+ "<",
+ (4*source_width)**2),
+ 1, 0)
+
def u_analytic(x):
return 0
from grudge.models.advection import VariableCoefficientAdvectionOperator
+ from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # noqa
+
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+ quad_tag_to_group_factory={
+ #"product": None,
+ "product": QuadratureSimplexGroupFactory(order=4*order)
+ })
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
op = VariableCoefficientAdvectionOperator(2, advec_v,
- inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
+ u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product",
flux_type=flux_type)
- bound_op = bind(discr, op.sym_operator())
+ bound_op = bind(
+ discr, op.sym_operator(),
+ #debug_flags=["dump_sym_operator_stages"]
+ )
- u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
+ u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0)
def rhs(t, u):
return bound_op(queue, t=t, u=u)
- final_time = 2
-
+ final_time = 50
dt = dt_factor * h/order**2
nsteps = (final_time // dt) + 1
dt = final_time/nsteps + 1e-15
@@ -90,7 +114,7 @@ def main(write_output=True, order=4):
dt_stepper = set_up_rk4("u", dt, u, rhs)
from grudge.shortcuts import make_visualizer
- vis = make_visualizer(discr, vis_order=order)
+ vis = make_visualizer(discr, vis_order=2*order)
step = 0
@@ -98,10 +122,11 @@ def main(write_output=True, order=4):
if isinstance(event, dt_stepper.StateComputed):
step += 1
+ if step % 30 == 0:
+ print(step)
- #print(step, event.t, norm(queue, u=event.state_component[0]))
- vis.write_vtk_file("fld-%04d.vtu" % step,
- [("u", event.state_component)])
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [("u", event.state_component)])
if __name__ == "__main__":
diff --git a/grudge/discretization.py b/grudge/discretization.py
index bc59299fefdbc7437a2efd9136c0a83bd4b71c04..74b2060138a888f40f35319df1d983f07e6e019a 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -1,6 +1,6 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -49,18 +49,22 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
.. automethod :: zeros
"""
- def __init__(self, cl_ctx, mesh, order, quad_min_degrees=None,
+ def __init__(self, cl_ctx, mesh, order, quad_tag_to_group_factory=None,
mpi_communicator=None):
"""
- :param quad_min_degrees: A mapping from quadrature tags to the degrees
- to which the desired quadrature is supposed to be exact.
+ :param quad_tag_to_group_factory: A mapping from quadrature tags (typically
+ strings--but may be any hashable/comparable object) to a
+ :class:`meshmode.discretization.ElementGroupFactory` indicating with
+ which quadrature discretization the operations are to be carried out,
+ or *None* to indicate that operations with this quadrature tag should
+ be carried out with the standard volume discretization.
"""
- if quad_min_degrees is None:
- quad_min_degrees = {}
+ if quad_tag_to_group_factory is None:
+ quad_tag_to_group_factory = {}
self.order = order
- self.quad_min_degrees = quad_min_degrees
+ self.quad_tag_to_group_factory = quad_tag_to_group_factory
from meshmode.discretization import Discretization
@@ -129,16 +133,26 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
if dd.is_volume():
if qtag is not sym.QTAG_NONE:
- # FIXME
- raise NotImplementedError("quadrature")
+ return self._quad_volume_discr(qtag)
return self._volume_discr
- elif dd.domain_tag is sym.FACE_RESTR_ALL:
- return self._all_faces_discr(qtag)
+ if qtag is not sym.QTAG_NONE:
+ no_quad_discr = self.discr_from_dd(sym.DOFDesc(dd.domain_tag))
+
+ from meshmode.discretization import Discretization
+ return Discretization(
+ self._volume_discr.cl_context,
+ no_quad_discr.mesh,
+ self.group_factory_for_quadrature_tag(qtag))
+
+ assert qtag is sym.QTAG_NONE
+
+ if dd.domain_tag is sym.FACE_RESTR_ALL:
+ return self._all_faces_volume_connection().to_discr
elif dd.domain_tag is sym.FACE_RESTR_INTERIOR:
- return self._interior_faces_discr(qtag)
+ return self._interior_faces_connection().to_discr
elif dd.is_boundary():
- return self._boundary_discr(dd.domain_tag, qtag)
+ return self._boundary_connection(dd.domain_tag).to_discr
else:
raise ValueError("DOF desc tag not understood: " + str(dd))
@@ -147,18 +161,32 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
from_dd = sym.as_dofdesc(from_dd)
to_dd = sym.as_dofdesc(to_dd)
- if from_dd.quadrature_tag is not sym.QTAG_NONE:
- raise ValueError("cannot interpolate *from* a "
- "(non-interpolatory) quadrature grid")
-
to_qtag = to_dd.quadrature_tag
+ if (
+ not from_dd.is_volume()
+ and from_dd.quadrature_tag == to_dd.quadrature_tag
+ and to_dd.domain_tag is sym.FACE_RESTR_ALL):
+ faces_conn = self.connection_from_dds(
+ sym.DOFDesc("vol"),
+ sym.DOFDesc(from_dd.domain_tag))
+
+ from meshmode.discretization.connection import \
+ make_face_to_all_faces_embedding
+
+ return make_face_to_all_faces_embedding(
+ faces_conn, self.discr_from_dd(to_dd),
+ self.discr_from_dd(from_dd))
+
# {{{ simplify domain + qtag change into chained
- if (from_dd.domain_tag != to_dd.domain_tag
- and from_dd.quadrature_tag != to_dd.quadrature_tag):
+ if (
+ from_dd.domain_tag != to_dd.domain_tag
+ and from_dd.quadrature_tag is sym.QTAG_NONE
+ and to_dd.quadrature_tag is not sym.QTAG_NONE):
- from meshmode.connection import ChainedDiscretizationConnection
+ from meshmode.discretization.connection import \
+ ChainedDiscretizationConnection
intermediate_dd = sym.DOFDesc(to_dd.domain_tag)
return ChainedDiscretizationConnection(
[
@@ -177,8 +205,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
# {{{ generic to-quad
- if (from_dd.domain_tag == to_dd.domain_tag
- and from_dd.quadrature_tag != to_dd.quadrature_tag):
+ if (
+ from_dd.domain_tag == to_dd.domain_tag
+ and from_dd.quadrature_tag is sym.QTAG_NONE
+ and to_dd.quadrature_tag is not sym.QTAG_NONE):
from meshmode.discretization.connection.same_mesh import \
make_same_mesh_connection
@@ -188,38 +218,30 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
# }}}
+ if from_dd.quadrature_tag is not sym.QTAG_NONE:
+ raise ValueError("cannot interpolate *from* a "
+ "(non-interpolatory) quadrature grid")
+
+ assert to_qtag is sym.QTAG_NONE
+
if from_dd.is_volume():
if to_dd.domain_tag is sym.FACE_RESTR_ALL:
- return self._all_faces_volume_connection(to_qtag)
+ return self._all_faces_volume_connection()
if to_dd.domain_tag is sym.FACE_RESTR_INTERIOR:
- return self._interior_faces_connection(to_qtag)
+ return self._interior_faces_connection()
elif to_dd.is_boundary():
assert from_dd.quadrature_tag is sym.QTAG_NONE
return self._boundary_connection(to_dd.domain_tag)
+ elif to_dd.is_volume():
+ from meshmode.discretization.connection import \
+ make_same_mesh_connection
+ to_discr = self._quad_volume_discr(to_dd.quadrature_tag)
+ from_discr = self._volume_discr
+ return make_same_mesh_connection(to_discr, from_discr)
else:
raise ValueError("cannot interpolate from volume to: " + str(to_dd))
- elif from_dd.domain_tag is sym.FACE_RESTR_INTERIOR:
- if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
- return self._all_faces_connection(None)
- else:
- raise ValueError(
- "cannot interpolate from interior faces to: "
- + str(to_dd))
-
- elif from_dd.domain_tag is sym.FACE_RESTR_ALL:
- if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
- return self._all_faces_connection(None)
-
- elif from_dd.is_boundary():
- if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
- return self._all_faces_connection(from_dd.domain_tag)
- else:
- raise ValueError(
- "cannot interpolate from interior faces to: "
- + str(to_dd))
-
else:
raise ValueError("cannot interpolate from: " + str(from_dd))
@@ -232,11 +254,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
quadrature_tag = sym.QTAG_NONE
from meshmode.discretization.poly_element import \
- PolynomialWarpAndBlendGroupFactory, \
- QuadratureSimplexGroupFactory
+ PolynomialWarpAndBlendGroupFactory
if quadrature_tag is not sym.QTAG_NONE:
- return QuadratureSimplexGroupFactory(order=self.order)
+ return self.quad_tag_to_group_factory[quadrature_tag]
else:
return PolynomialWarpAndBlendGroupFactory(order=self.order)
@@ -257,33 +278,17 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
boundary_tag=boundary_tag)
- @memoize_method
- def _boundary_discr(self, boundary_tag, quadrature_tag=None):
- if quadrature_tag is None:
- quadrature_tag = sym.QTAG_NONE
-
- if quadrature_tag is sym.QTAG_NONE:
- return self._boundary_connection(boundary_tag).to_discr
- else:
- no_quad_bdry_discr = self.boundary_discr(boundary_tag, sym.QTAG_NONE)
-
- from meshmode.discretization import Discretization
- return Discretization(
- self._volume_discr.cl_context,
- no_quad_bdry_discr.mesh,
- self.group_factory_for_quadrature_tag(quadrature_tag))
-
# }}}
# {{{ interior faces
@memoize_method
- def _interior_faces_connection(self, quadrature_tag=None):
+ def _interior_faces_connection(self):
from meshmode.discretization.connection import (
make_face_restriction, FACE_RESTR_INTERIOR)
return make_face_restriction(
self._volume_discr,
- self.group_factory_for_quadrature_tag(quadrature_tag),
+ self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
FACE_RESTR_INTERIOR,
# FIXME: This will need to change as soon as we support
@@ -291,32 +296,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
# types.
per_face_groups=False)
- def _interior_faces_discr(self, quadrature_tag=None):
- return self._interior_faces_connection(quadrature_tag).to_discr
-
@memoize_method
- def opposite_face_connection(self, quadrature_tag):
- if quadrature_tag is not sym.QTAG_NONE:
- # FIXME
- raise NotImplementedError("quadrature")
-
+ def opposite_face_connection(self):
from meshmode.discretization.connection import \
make_opposite_face_connection
- return make_opposite_face_connection(
- self._interior_faces_connection(quadrature_tag))
+ return make_opposite_face_connection(self._interior_faces_connection())
# }}}
# {{{ all-faces
@memoize_method
- def _all_faces_volume_connection(self, quadrature_tag=None):
+ def _all_faces_volume_connection(self):
from meshmode.discretization.connection import (
make_face_restriction, FACE_RESTR_ALL)
return make_face_restriction(
self._volume_discr,
- self.group_factory_for_quadrature_tag(quadrature_tag),
+ self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
FACE_RESTR_ALL,
# FIXME: This will need to change as soon as we support
@@ -324,30 +321,6 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
# types.
per_face_groups=False)
- def _all_faces_discr(self, quadrature_tag=None):
- return self._all_faces_volume_connection(quadrature_tag).to_discr
-
- @memoize_method
- def _all_faces_connection(self, boundary_tag):
- """Return a
- :class:`meshmode.discretization.connection.DiscretizationConnection`
- that goes from either
- :meth:`_interior_faces_discr` (if *boundary_tag* is None)
- or
- :meth:`_boundary_discr` (if *boundary_tag* is not None)
- to a discretization containing all the faces of the volume
- discretization.
- """
- from meshmode.discretization.connection import \
- make_face_to_all_faces_embedding
-
- if boundary_tag is None:
- faces_conn = self._interior_faces_connection()
- else:
- faces_conn = self._boundary_connection(boundary_tag)
-
- return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr())
-
# }}}
@property
diff --git a/grudge/execution.py b/grudge/execution.py
index 1fdec1b9aa6ea3942d94e5ce5010996aa6ce2eb6..cc8703dbf75cd84f0908bc2166eca762b1ad5e42 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -1,6 +1,6 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -201,36 +201,37 @@ class ExecutionMapper(mappers.Evaluator,
knl = lp.make_kernel(
"""{[k,i,j]:
0<=k<nelements and
- 0<=i,j<ndiscr_nodes}""",
+ 0<=i<ndiscr_nodes_out and
+ 0<=j<ndiscr_nodes_in}""",
"result[k,i] = sum(j, mat[i, j] * vec[k, j])",
default_offset=lp.auto, name="diff")
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
- discr = self.discrwb.discr_from_dd(op.dd_in)
+ in_discr = self.discrwb.discr_from_dd(op.dd_in)
+ out_discr = self.discrwb.discr_from_dd(op.dd_out)
- # FIXME: This shouldn't really assume that it's dealing with a volume
- # input. What about quadrature? What about boundaries?
- result = discr.empty(
+ result = out_discr.empty(
queue=self.queue,
dtype=field.dtype, allocator=self.bound_op.allocator)
- for grp in discr.groups:
- cache_key = "elwise_linear", grp, op, field.dtype
+ for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
+
+ cache_key = "elwise_linear", in_grp, out_grp, op, field.dtype
try:
matrix = self.bound_op.operator_data_cache[cache_key]
except KeyError:
matrix = (
- cl.array.to_device(
- self.queue,
- np.asarray(op.matrix(grp), dtype=field.dtype))
- .with_queue(None))
+ cl.array.to_device(
+ self.queue,
+ np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype))
+ .with_queue(None))
self.bound_op.operator_data_cache[cache_key] = matrix
- knl()(self.queue, mat=matrix, result=grp.view(result),
- vec=grp.view(field))
+ knl()(self.queue, mat=matrix, result=out_grp.view(result),
+ vec=in_grp.view(field))
return result
@@ -273,14 +274,7 @@ class ExecutionMapper(mappers.Evaluator,
return shuffled_recv_vec
def map_opposite_interior_face_swap(self, op, field_expr):
- dd = op.dd_in
-
- qtag = dd.quadrature_tag
- if qtag is None:
- # FIXME: Remove once proper quadrature support arrives
- qtag = sym.QTAG_NONE
-
- return self.discrwb.opposite_face_connection(qtag)(
+ return self.discrwb.opposite_face_connection()(
self.queue, self.rec(field_expr)).with_queue(self.queue)
# {{{ face mass operator
@@ -316,8 +310,7 @@ class ExecutionMapper(mappers.Evaluator,
assert len(all_faces_discr.groups) == len(vol_discr.groups)
- for cgrp, afgrp, volgrp in zip(all_faces_conn.groups,
- all_faces_discr.groups, vol_discr.groups):
+ for afgrp, volgrp in zip(all_faces_discr.groups, vol_discr.groups):
cache_key = "face_mass", afgrp, op, field.dtype
nfaces = volgrp.mesh_el_group.nfaces
@@ -392,10 +385,6 @@ class ExecutionMapper(mappers.Evaluator,
# This should be unified with map_elementwise_linear, which should
# be extended to support batching.
- discr = self.discrwb.discr_from_dd(repr_op.dd_in)
-
- # FIXME: Enable
- # assert repr_op.dd_in == repr_op.dd_out
assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
@memoize_in(self.discrwb, "reference_derivative_knl")
@@ -404,7 +393,8 @@ class ExecutionMapper(mappers.Evaluator,
"""{[imatrix,k,i,j]:
0<=imatrix<nmatrices and
0<=k<nelements and
- 0<=i,j<nunit_nodes}""",
+ 0<=i<nunit_nodes_out and
+ 0<=j<nunit_nodes_in}""",
"""
result[imatrix, k, i] = sum(
j, diff_mat[imatrix, i, j] * vec[k, j])
@@ -415,26 +405,31 @@ class ExecutionMapper(mappers.Evaluator,
return lp.tag_inames(knl, dict(k="g.0"))
noperators = len(insn.operators)
- result = discr.empty(
+
+ in_discr = self.discrwb.discr_from_dd(repr_op.dd_in)
+ out_discr = self.discrwb.discr_from_dd(repr_op.dd_out)
+
+ result = out_discr.empty(
queue=self.queue,
dtype=field.dtype, extra_dims=(noperators,),
allocator=self.bound_op.allocator)
- for grp in discr.groups:
- if grp.nelements == 0:
+ for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
+
+ if in_grp.nelements == 0:
continue
- matrices = repr_op.matrices(grp)
+ matrices = repr_op.matrices(out_grp, in_grp)
# FIXME: Should transfer matrices to device and cache them
matrices_ary = np.empty((
- noperators, grp.nunit_nodes, grp.nunit_nodes))
+ noperators, out_grp.nunit_nodes, in_grp.nunit_nodes))
for i, op in enumerate(insn.operators):
matrices_ary[i] = matrices[op.rst_axis]
knl()(self.queue,
diff_mat=matrices_ary,
- result=grp.view(result), vec=grp.view(field))
+ result=out_grp.view(result), vec=in_grp.view(field))
return [(name, result[i]) for i, name in enumerate(insn.names)], []
@@ -522,6 +517,10 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
dumper("before-cfold", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
+ dumper("before-qcheck", sym_operator)
+ sym_operator = mappers.QuadratureCheckerAndRemover(
+ discrwb.quad_tag_to_group_factory)(sym_operator)
+
# Work around https://github.com/numpy/numpy/issues/9438
#
# The idea is that we need 1j as an expression to survive
@@ -537,26 +536,6 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
complex_type=discrwb.complex_dtype.type,
)(sym_operator)
- # Ordering restriction:
- #
- # - Must run constant fold before first type inference pass, because zeros,
- # while allowed, violate typing constraints (because they can't be assigned
- # a unique type), and need to be killed before the type inferrer sees them.
-
- # FIXME: Reenable type inference
-
- # from grudge.symbolic.mappers.type_inference import TypeInferrer
- # dumper("before-specializer", sym_operator)
- # sym_operator = mappers.OperatorSpecializer(
- # TypeInferrer()(sym_operator)
- # )(sym_operator)
-
- # Ordering restriction:
- #
- # - Must run OperatorSpecializer before performing the GlobalToReferenceMapper,
- # because otherwise it won't differentiate the type of grids (node or quadrature
- # grids) that the operators will apply on.
-
dumper("before-global-to-reference", sym_operator)
sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator)
@@ -567,12 +546,6 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
connected_parts = get_connected_partitions(volume_mesh)
sym_operator = mappers.DistributedMapper(connected_parts)(sym_operator)
- # Ordering restriction:
- #
- # - Must specialize quadrature operators before performing inverse mass
- # contraction, because there are no inverse-mass-contracted variants of the
- # quadrature operators.
-
dumper("before-imass", sym_operator)
sym_operator = mappers.InverseMassContractor()(sym_operator)
@@ -598,19 +571,20 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
stage = [0]
- def dump_optemplate(name, sym_operator):
- if "dump_optemplate_stages" in debug_flags:
- from grudge.tools import open_unique_debug_file
- from grudge.optemplate import pretty
- open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt").write(
- pretty(sym_operator))
+ def dump_sym_operator(name, sym_operator):
+ if "dump_sym_operator_stages" in debug_flags:
+ from pytools.debug import open_unique_debug_file
+ outf, name = open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt")
+ with outf:
+ outf.write(sym.pretty(sym_operator))
+
stage[0] += 1
sym_operator = process_sym_operator(
discr,
sym_operator,
post_bind_mapper=post_bind_mapper,
- dumper=dump_optemplate)
+ dumper=dump_sym_operator)
from grudge.symbolic.compiler import OperatorCompiler
discr_code, eval_code = OperatorCompiler(discr)(sym_operator)
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 1c83ed9675d5a45c75b22e3b4d5ca7926acf9133..091b9870c4b82a3fde2fc21b84bf92e5ec62ab17 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -3,7 +3,7 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -132,13 +132,16 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
# {{{ variable-coefficient advection
class VariableCoefficientAdvectionOperator(HyperbolicOperator):
- def __init__(self, dim, v, inflow_u, flux_type="central"):
+ def __init__(self, dim, v, inflow_u, flux_type="central", quad_tag="product"):
self.ambient_dim = dim
self.v = v
self.inflow_u = inflow_u
self.flux_type = flux_type
+ self.quad_tag = quad_tag
+
def flux(self, u):
+
normal = sym.normal(u. dd, self.ambient_dim)
surf_v = sym.interp("vol", u.dd)(self.v)
@@ -167,15 +170,26 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
# boundary conditions -------------------------------------------------
bc_in = self.inflow_u
+ all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
+ boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+
def flux(pair):
- return sym.interp(pair.dd, "all_faces")(
+ return sym.interp(pair.dd, all_faces_dd)(
self.flux(pair))
+ quad_dd = sym.DOFDesc("vol", self.quad_tag)
+
+ to_quad = sym.interp("vol", quad_dd)
+
+ stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol")
+ quad_v = to_quad(self.v)
+ quad_u = to_quad(u)
+
return sym.InverseMassOperator()(
- np.dot(sym.stiffness_t(self.ambient_dim), sym.cse(self.v*u))
- - sym.FaceMassOperator()(
- flux(sym.int_tpair(u))
- + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in))
+ (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1]))
+ - sym.FaceMassOperator(all_faces_dd, "vol")(
+ flux(sym.int_tpair(u, self.quad_tag))
+ + flux(sym.bv_tpair(boundary_dd, u, bc_in))
# FIXME: Add back support for inflow/outflow tags
#+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index f6e238083266d4ba0c510243b116b2e49b024631..c555cea0a56e947b4a94e0cd1736e8fa2b97ce63 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -338,7 +338,7 @@ class Code(object):
self.static_schedule_attempts = 5
def dump_dataflow_graph(self):
- from grudge.tools import open_unique_debug_file
+ from pytools.debug import open_unique_debug_file
open_unique_debug_file("dataflow", ".dot")\
.write(dot_dataflow_graph(self, max_node_label_length=None))
@@ -792,7 +792,7 @@ def aggregate_assignments(inf_mapper, instructions, result,
# }}}
-# {{{
+# {{{ to-loopy mapper
def set_once(d, k, v):
try:
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index a810a335b5ee144ea9fe6fe4aca837781fcace5a..2ddd6f5d41f0c4a6be33c4dddb4a744824d47a86 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -798,96 +798,67 @@ class NoCSEStringifyMapper(StringifyMapper):
# {{{ quadrature support
-class QuadratureUpsamplerRemover(CSECachingMapperMixin, IdentityMapper):
- def __init__(self, quad_min_degrees, do_warn=True):
+class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
+ """Checks whether all quadratu
+ """
+
+ def __init__(self, quad_tag_to_group_factory):
IdentityMapper.__init__(self)
CSECachingMapperMixin.__init__(self)
- self.quad_min_degrees = quad_min_degrees
- self.do_warn = do_warn
+ self.quad_tag_to_group_factory = quad_tag_to_group_factory
map_common_subexpression_uncached = \
IdentityMapper.map_common_subexpression
- def map_operator_binding(self, expr):
- if isinstance(expr.op, (op.QuadratureGridUpsampler,
- op.QuadratureInteriorFacesGridUpsampler,
- op.QuadratureBoundaryGridUpsampler)):
- try:
- min_degree = self.quad_min_degrees[expr.op.quadrature_tag]
- except KeyError:
- if self.do_warn:
- from warnings import warn
- warn("No minimum degree for quadrature tag '%s' specified--"
- "falling back to nodal evaluation"
- % expr.op.quadrature_tag)
- return self.rec(expr.field)
- else:
- if min_degree is None:
- return self.rec(expr.field)
- else:
- return expr.op(self.rec(expr.field))
- else:
- return IdentityMapper.map_operator_binding(self, expr)
+ def _process_dd(self, dd, location_descr):
+ from grudge.symbolic.primitives import DOFDesc, QTAG_NONE
+ if dd.quadrature_tag is not QTAG_NONE:
+ if dd.quadrature_tag not in self.quad_tag_to_group_factory:
+ raise ValueError("found unknown quadrature tag '%s' in '%s'"
+ % (dd.quadrature_tag, location_descr))
+ grp_factory = self.quad_tag_to_group_factory[dd.quadrature_tag]
+ if grp_factory is None:
+ dd = DOFDesc(dd.domain_tag, QTAG_NONE)
-class QuadratureDetector(CSECachingMapperMixin, CombineMapper):
- """For a given expression, this mapper returns the upsampling
- operator in effect at the root of the expression, or *None*
- if there isn't one.
- """
- class QuadStatusNotKnown:
- pass
-
- map_common_subexpression_uncached = \
- CombineMapper.map_common_subexpression
+ return dd
- def combine(self, values):
- from pytools import single_valued
- return single_valued([
- v for v in values if v is not self.QuadStatusNotKnown])
+ def map_operator_binding(self, expr):
+ dd_in_orig = dd_in = expr.op.dd_in
+ dd_out_orig = dd_out = expr.op.dd_out
- def map_variable(self, expr):
- return None
+ dd_in = self._process_dd(dd_in, "dd_in of %s" % type(expr.op).__name__)
+ dd_out = self._process_dd(dd_out, "dd_out of %s" % type(expr.op).__name__)
- def map_constant(self, expr):
- return self.QuadStatusNotKnown
+ if dd_in_orig == dd_in and dd_out_orig == dd_out:
+ # unchanged
+ return IdentityMapper.map_operator_binding(self, expr)
- def map_operator_binding(self, expr):
- if isinstance(expr.op, (
- op.DiffOperatorBase, op.FluxOperatorBase,
- op.MassOperatorBase)):
- return None
- elif isinstance(expr.op, (op.QuadratureGridUpsampler,
- op.QuadratureInteriorFacesGridUpsampler)):
- return expr.op
- else:
- return CombineMapper.map_operator_binding(self, expr)
+ import grudge.symbolic.operators as op
+ # changed
+ if dd_in == dd_out and isinstance(expr.op, op.InterpolationOperator):
+ # This used to be to-quad interpolation and has become a no-op.
+ # Remove it.
+ return self.rec(expr.field)
-class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper):
- """This mapper descends the expression tree, down to each
- quadrature-consuming operator (diff, mass) along each branch.
- It then change
- """
- def __init__(self, desired_quad_op):
- IdentityMapper.__init__(self)
- CSECachingMapperMixin.__init__(self)
+ if isinstance(expr.op, op.StiffnessTOperator):
+ new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out)
+ elif isinstance(expr.op, (op.FaceMassOperator, op.InterpolationOperator)):
+ new_op = type(expr.op)(dd_in, dd_out)
+ else:
+ raise NotImplementedError("do not know how to modify dd_in and dd_out "
+ "in %s" % type(expr.op).__name__)
- self.desired_quad_op = desired_quad_op
+ return new_op(self.rec(expr.field))
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
+ def map_ones(self, expr):
+ dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
+ return type(expr)(dd)
- def map_operator_binding(self, expr):
- if isinstance(expr.op, (
- op.DiffOperatorBase, op.FluxOperatorBase,
- op.MassOperatorBase)):
- return expr
- elif isinstance(expr.op, (op.QuadratureGridUpsampler,
- op.QuadratureInteriorFacesGridUpsampler)):
- return self.desired_quad_op(expr.field)
- else:
- return IdentityMapper.map_operator_binding(self, expr)
+ def map_node_coordinate_component(self, expr):
+ dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
+ return type(expr)(expr.axis, dd)
# }}}
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 76bd938c941d5b7e0eea5682f164f188f3e2aff0..294c4374e275ac65666c3ae2d6d1442f81e491a7 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -2,7 +2,7 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -97,6 +97,14 @@ class ElementwiseLinearOperator(Operator):
class InterpolationOperator(Operator):
+ def __init__(self, dd_in, dd_out):
+ official_dd_in = _sym().as_dofdesc(dd_in)
+ official_dd_out = _sym().as_dofdesc(dd_out)
+ if official_dd_in == official_dd_out:
+ raise ValueError("Interpolating from {} to {}"
+ " does not do anything.".format(official_dd_in, official_dd_out))
+
+ super(InterpolationOperator, self).__init__(dd_in, dd_out)
mapper_method = intern("map_interpolation")
@@ -141,8 +149,12 @@ class DiffOperatorBase(Operator):
def __init__(self, xyz_axis, dd_in=None, dd_out=None):
if dd_in is None:
dd_in = _sym().DD_VOLUME
+
if dd_out is None:
dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+ else:
+ dd_out = _sym().as_dofdesc(dd_out)
+
if dd_out.uses_quadrature():
raise ValueError("differentiation outputs are not on "
"quadrature grids")
@@ -215,18 +227,30 @@ class RefDiffOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_diff")
@staticmethod
- def matrices(element_group):
- return element_group.diff_matrices()
+ def matrices(out_element_group, in_element_group):
+ assert in_element_group == out_element_group
+ return in_element_group.diff_matrices()
class RefStiffnessTOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_stiffness_t")
@staticmethod
- def matrices(element_group):
- assert element_group.is_orthogonal_basis()
- mmat = element_group.mass_matrix()
- return [dmat.T.dot(mmat.T) for dmat in element_group.diff_matrices()]
+ def matrices(out_elem_grp, in_elem_grp):
+ if in_elem_grp == out_elem_grp:
+ assert in_elem_grp.is_orthogonal_basis()
+ mmat = in_elem_grp.mass_matrix()
+ return [dmat.T.dot(mmat.T) for dmat in in_elem_grp.diff_matrices()]
+
+ from modepy import vandermonde
+ vand = vandermonde(out_elem_grp.basis(), out_elem_grp.unit_nodes)
+ grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes)
+ vand_inv_t = np.linalg.inv(vand).T
+
+ weights = in_elem_grp.weights
+
+ return np.einsum('c,bz,acz->abc', weights, vand_inv_t, grad_vand)
+
# }}}
@@ -322,8 +346,6 @@ class VandermondeOperator(ElementwiseLinearOperator):
# }}}
-# }}}
-
# {{{ mass operators
@@ -338,9 +360,6 @@ class MassOperatorBase(Operator):
if dd_out is None:
dd_out = dd_in
- if dd_out != dd_in:
- raise ValueError("dd_out and dd_in must be identical")
-
super(MassOperatorBase, self).__init__(dd_in, dd_out)
@@ -358,19 +377,30 @@ class RefMassOperatorBase(ElementwiseLinearOperator):
class RefMassOperator(RefMassOperatorBase):
@staticmethod
- def matrix(element_group):
- return element_group.mass_matrix()
+ def matrix(out_element_group, in_element_group):
+ if out_element_group == in_element_group:
+ return in_element_group.mass_matrix()
+
+ from modepy import vandermonde
+ vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes)
+ o_vand = vandermonde(out_element_group.basis(), in_element_group.unit_nodes)
+ vand_inv_t = np.linalg.inv(vand).T
+
+ weights = in_element_group.weights
+
+ return np.einsum('c,bz,acz->abc', weights, vand_inv_t, o_vand)
mapper_method = intern("map_ref_mass")
class RefInverseMassOperator(RefMassOperatorBase):
@staticmethod
- def matrix(element_group):
+ def matrix(in_element_group, out_element_group):
+ assert in_element_group == out_element_group
import modepy as mp
return mp.inverse_mass_matrix(
- element_group.basis(),
- element_group.unit_nodes)
+ in_element_group.basis(),
+ in_element_group.unit_nodes)
mapper_method = intern("map_ref_inverse_mass")
@@ -426,7 +456,7 @@ class FaceMassOperatorBase(ElementwiseLinearOperator):
if dd_in is None:
dd_in = sym.DOFDesc(sym.FACE_RESTR_ALL, None)
- if dd_out is None:
+ if dd_out is None or dd_out == "vol":
dd_out = sym.DOFDesc("vol", sym.QTAG_NONE)
if dd_out.uses_quadrature():
raise ValueError("face mass operator outputs are not on "
@@ -497,10 +527,10 @@ def stiffness(dim):
[StiffnessOperator(i) for i in range(dim)])
-def stiffness_t(dim):
+def stiffness_t(dim, dd_in=None, dd_out=None):
from pytools.obj_array import make_obj_array
return make_obj_array(
- [StiffnessTOperator(i) for i in range(dim)])
+ [StiffnessTOperator(i, dd_in, dd_out) for i in range(dim)])
def integral(arg, dd=None):
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 3281abef525697ee692aacba76a18ba39a574664..974833c2a1e49587817fdf03c3a129210b9562c7 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -2,7 +2,7 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -657,10 +657,18 @@ class TracePair:
return 0.5*(self.int + self.ext)
-def int_tpair(expression):
- i = cse(_sym().interp("vol", "int_faces")(expression))
+def int_tpair(expression, qtag=None):
+ i = _sym().interp("vol", "int_faces")(expression)
e = cse(_sym().OppositeInteriorFaceSwap()(i))
- return TracePair("int_faces", i, e)
+
+ if qtag is not None and qtag != _sym().QTAG_NONE:
+ q_dd = _sym().DOFDesc("int_faces", qtag)
+ i = cse(_sym().interp("int_faces", q_dd)(i))
+ e = cse(_sym().interp("int_faces", q_dd)(e))
+ else:
+ q_dd = "int_faces"
+
+ return TracePair(q_dd, i, e)
def bdry_tpair(dd, interior, exterior):
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 84c46168ede26efe1a42935501d69fba99b7700a..34532bfc80371e29bde28ea9674c48e36ab889aa 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -334,8 +334,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
@pytest.mark.parametrize("order", [3, 4, 5])
-# test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)'
-def test_convergence_maxwell(ctx_factory, order, visualize=False):
+def test_convergence_maxwell(ctx_factory, order):
"""Test whether 3D maxwells actually converges"""
cl_ctx = cl.create_some_context()
@@ -403,6 +402,73 @@ def test_convergence_maxwell(ctx_factory, order, visualize=False):
assert eoc_rec.order_estimate() > order
+@pytest.mark.parametrize("order", [2, 3, 4])
+def test_improvement_quadrature(ctx_factory, order):
+ """Test whether quadrature improves things and converges"""
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ from grudge.models.advection import VariableCoefficientAdvectionOperator
+ from pytools.convergence import EOCRecorder
+ from pytools.obj_array import join_fields
+ from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 2
+ sym_nds = sym.nodes(dims)
+ advec_v = join_fields(-1*sym_nds[1], sym_nds[0])
+
+ flux = "upwind"
+ op = VariableCoefficientAdvectionOperator(2, advec_v, 0, flux_type=flux)
+
+ def gaussian_mode():
+ source_width = 0.1
+ sym_x = sym.nodes(2)
+ return sym.exp(-np.dot(sym_x, sym_x) / source_width**2)
+
+ def conv_test(descr, use_quad):
+ print("-"*75)
+ print(descr)
+ print("-"*75)
+ eoc_rec = EOCRecorder()
+
+ ns = [20, 25]
+ for n in ns:
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*dims,
+ b=(0.5,)*dims,
+ n=(n,)*dims,
+ order=order)
+
+ if use_quad:
+ quad_tag_to_group_factory = {
+ "product": QuadratureSimplexGroupFactory(order=4*order)
+ }
+ else:
+ quad_tag_to_group_factory = {"product": None}
+
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+ quad_tag_to_group_factory=quad_tag_to_group_factory)
+
+ bound_op = bind(discr, op.sym_operator())
+ fields = bind(discr, gaussian_mode())(queue, t=0)
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ esc = bound_op(queue, u=fields)
+ total_error = norm(queue, u=esc)
+ eoc_rec.add_data_point(1.0/n, total_error)
+
+ print(eoc_rec.pretty_print(abscissa_label="h", error_label="LInf Error"))
+
+ return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history])
+
+ eoc, errs = conv_test("no quadrature", False)
+ q_eoc, q_errs = conv_test("with quadrature", True)
+ assert q_eoc > eoc
+ assert (q_errs < errs).all()
+ assert q_eoc > order
+
+
def test_foreign_points(ctx_factory):
pytest.importorskip("sumpy")
import sumpy.point_calculus as pc