diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index 0a69967e8613adc76a9ff1d93ba02db90fb2dc9a..97e793fd9f2d27acc31e4e4cdb8d04eb57ebce29 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -30,6 +30,8 @@ logger = logging.getLogger(__name__)
from grudge import sym, bind, DGDiscretizationWithBoundaries
from pytools.obj_array import join_fields
+import meshmode.discretization.poly_element as poly_element
+
def main(write_output=True, order=4):
cl_ctx = cl.create_some_context()
@@ -37,15 +39,15 @@ def main(write_output=True, order=4):
dim = 2
- n_divs = 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=(n_divs, n_divs), order=order)
+ n=(10, 10), order=order)
dt_factor = 5
- h = 1/n_divs
+ h = 1/10
+ final_time = 5
+
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
sym_x = sym.nodes(2)
@@ -69,10 +71,11 @@ def main(write_output=True, order=4):
from grudge.models.advection import VariableCoefficientAdvectionOperator
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2*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)
+ #op = sym.interp("vol", sym.DOFDesc("vol", "product"))(sym.var("u"))
bound_op = bind(discr, op.sym_operator())
@@ -81,7 +84,6 @@ def main(write_output=True, order=4):
def rhs(t, u):
return bound_op(queue, t=t, u=u)
- final_time = 2
dt = dt_factor * h/order**2
nsteps = (final_time // dt) + 1
@@ -99,6 +101,7 @@ def main(write_output=True, order=4):
if isinstance(event, dt_stepper.StateComputed):
step += 1
+ print(step)
#print(step, event.t, norm(queue, u=event.state_component[0]))
vis.write_vtk_file("fld-%04d.vtu" % step,
diff --git a/grudge/execution.py b/grudge/execution.py
index 4a51477c5af3d05f278f4e41f3a39dd242e078e2..b950a961d3605070c83c554f78efa4976a433000 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -252,8 +252,8 @@ class ExecutionMapper(mappers.Evaluator,
return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
def map_opposite_interior_face_swap(self, op, field_expr):
- if op.dd_in.quadrature_tag is not sym.QTAG_NONE:
- raise ValueError("face swap unavailable on quadrature meshes")
+ #if op.dd_in.quadrature_tag is not sym.QTAG_NONE:
+ #raise ValueError("face swap unavailable on quadrature meshes")
return self.discrwb.opposite_face_connection()(
self.queue, self.rec(field_expr)).with_queue(self.queue)
@@ -291,8 +291,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
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index e0203378214983cceb11ecb46d6230a2dfc7574e..4d67883203653c12664a27caa3d98ec1e4155793 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -132,13 +132,18 @@ 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
+ if quad_tag is None:
+ self.quad_tag = sym.QTAG_NONE
+ else:
+ 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 +172,25 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
# boundary conditions -------------------------------------------------
bc_in = self.inflow_u
+ all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_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)
+
+ if self.quad_tag == sym.QTAG_NONE:
+ to_quad = lambda x : x
+ else:
+ to_quad = sym.interp("vol", quad_dd)
+
+ from_quad_stiffness_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol")
+
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))
+ np.dot(from_quad_stiffness_t, to_quad(self.v)*to_quad(u))
+ - sym.FaceMassOperator(all_faces_dd, "vol")(
+ flux(sym.int_tpair(u, self.quad_tag))
+ + flux(sym.bv_tpair(sym.DOFDesc(sym.BTAG_ALL, self.quad_tag), 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/operators.py b/grudge/symbolic/operators.py
index 3426b8e3c9ef0cdeda1bdbe3196553b62e91ff55..ff909f7f15f6f790ad93100c7276351f55c049d3 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -97,12 +97,15 @@ class ElementwiseLinearOperator(Operator):
class InterpolationOperator(Operator):
- #def __init__(self, dd_in, dd_out):
- #print("ayy")
- #print(dd_in)
- #print(dd_out)
-
- #super().__init__(dd_in, dd_out)
+ 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:
+ from warnings import warn
+ raise ValueError("Interpolating from {} to {}"
+ " does not do anything.".format(official_dd_in, official_dd_out))
+
+ super().__init__(dd_in, dd_out)
mapper_method = intern("map_interpolation")
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 91122a89e215ba26c6ea22b39f7efd21ea49f47e..9234c225ad953f7ad52aeee6252c29d6de827d0f 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -654,10 +654,17 @@ class TracePair:
def int_tpair(expression, qtag=None):
- dd = _sym().DOFDesc("int_faces", qtag)
- i = cse(_sym().interp("vol", dd)(expression))
+ i = _sym().interp("vol", "int_faces")(expression)
e = cse(_sym().OppositeInteriorFaceSwap()(i))
- return TracePair(dd, 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)
#i = cse(_sym().interp("vol", "int_faces")(expression))
#e = cse(_sym().OppositeInteriorFaceSwap()(i))