From 27d0a8256e5ac865c068bcc997ec2709cf82647b Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Sun, 21 Jun 2020 19:29:31 -0500
Subject: [PATCH] make example work with oversampling
---
examples/advection/surface.py | 33 ++++++---------------------------
grudge/models/advection.py | 18 ++++++++++--------
grudge/symbolic/operators.py | 1 +
grudge/symbolic/primitives.py | 23 ++++++++++++++---------
4 files changed, 31 insertions(+), 44 deletions(-)
diff --git a/examples/advection/surface.py b/examples/advection/surface.py
index 2ef5d770..a198dabf 100644
--- a/examples/advection/surface.py
+++ b/examples/advection/surface.py
@@ -50,7 +50,8 @@ class Plotter:
volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
x = volume_discr.nodes().with_queue(queue)
- self.x = (2.0 * np.pi * (x[1] < 0) + cl.clmath.atan2(x[1], x[0])).get(queue)
+ self.x = (2.0 * np.pi * (x[1] < 0)
+ + cl.clmath.atan2(x[1], x[0])).get(queue)
elif self.ambient_dim == 3:
from grudge.shortcuts import make_visualizer
self.vis = make_visualizer(discr, vis_order=order)
@@ -98,7 +99,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
# sphere radius
radius = 1.0
# sphere resolution
- resolution = 64 if dim == 2 else 1.0
+ resolution = 64 if dim == 2 else 1
# cfl
dt_factor = 1.0
@@ -110,7 +111,6 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
c = make_obj_array([
-sym_x[1], sym_x[0], 0.0
])[:dim]
- norm_c = sym.sqrt((c**2).sum())
# flux
flux_type = "lf"
@@ -125,29 +125,9 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
np.linspace(0.0, 1.0, resolution + 1),
order)
elif dim == 3:
- if 0:
- from meshmode.mesh.generation import generate_icosphere
- mesh = generate_icosphere(radius, order=4 * order)
- from meshmode.mesh.refinement import refine_uniformly
- mesh = refine_uniformly(mesh, resolution, with_adjacency=False)
- else:
- from meshmode.mesh.io import ScriptSource
- source = ScriptSource("""
- SetFactory("OpenCASCADE");
- Sphere(1) = {0, 0, 0, %g};
- Mesh.Algorithm = 6;
- """ % radius, "geo")
-
- from meshmode.mesh.io import generate_gmsh
- mesh = generate_gmsh(source, 2, order=order,
- other_options=[
- "-optimize_ho",
- "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution
- ],
- target_unit="MM")
-
- from meshmode.mesh.processing import perform_flips
- mesh = perform_flips(mesh, np.ones(mesh.nelements))
+ from meshmode.mesh.generation import generate_icosphere
+ mesh = generate_icosphere(radius, order=4 * order,
+ uniform_refinement_rounds=resolution)
else:
raise ValueError("unsupported dimension")
@@ -174,7 +154,6 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
# {{{ symbolic operators
-
def f_gaussian(x):
return x[0]
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 3891c9ca..88c09668 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -186,13 +186,15 @@ class VariableCoefficientAdvectionOperator(AdvectionOperatorBase):
# {{{ closed surface advection
def v_dot_n_tpair(velocity, dd=None):
- ambient_dim = len(velocity)
- normal = sym.normal(dd, ambient_dim, dim=ambient_dim - 2)
+ if dd is None:
+ dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR)
- v_dot_n_i = sym.cse(velocity.dot(normal))
- v_dot_n_e = sym.cse(sym.OppositeInteriorFaceSwap()(v_dot_n_i))
+ ambient_dim = len(velocity)
+ normal = sym.normal(dd.with_qtag(None), ambient_dim, dim=ambient_dim - 2)
- return sym.TracePair(dd, v_dot_n_i, v_dot_n_e)
+ 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):
@@ -225,11 +227,11 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase):
self.quad_tag = quad_tag
def flux(self, u):
- surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v)
+ surf_v = sym.interp(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v)
return surface_advection_weak_flux(self.flux_type, u, surf_v)
def penalty(self, u):
- surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v)
+ surf_v = sym.interp(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v)
return surface_penalty_flux(u, surf_v, tau=0.0)
def sym_operator(self):
@@ -242,8 +244,8 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase):
return sym.interp(pair.dd, face_dd)(self.penalty(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.interp(sym.DD_VOLUME, quad_dd)
stiff_t_op = sym.stiffness_t(self.ambient_dim,
dd_in=quad_dd, dd_out=sym.DD_VOLUME)
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 77f735b9..8f7609c0 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -742,6 +742,7 @@ 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
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index ceef7a10..7e26223b 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -867,20 +867,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 interp, OppositeInteriorFaceSwap
- i = interp("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(interp("int_faces", q_dd)(i))
- e = cse(interp("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 = interp(from_dd, trace_dd.with_qtag(None))(expression)
+ e = cse(OppositeInteriorFaceSwap()(i))
+
+ if trace_dd.uses_quadrature():
+ i = cse(interp(trace_dd.with_qtag(None), trace_dd)(i))
+ e = cse(interp(trace_dd.with_qtag(None), trace_dd)(e))
- return TracePair(q_dd, i, e)
+ return TracePair(trace_dd, i, e)
def bdry_tpair(dd, interior, exterior):
--
GitLab