diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index be95db4bb4846b78d02ef7a4a855a364d80bc77c..3b5f5d32d51cf402b15e34b6a086319f2eb083d5 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -26,13 +26,14 @@ THE SOFTWARE.
import numpy as np
+import pyopencl as cl
+from grudge.shortcuts import make_discretization, set_up_rk4
def main(write_output=True):
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5))
- from grudge.shortcuts import make_discretization
discr = make_discretization(mesh, order=4)
#from grudge.visualization import VtkVisualizer
@@ -42,16 +43,19 @@ def main(write_output=True):
source_width = 0.05
source_omega = 3
- import grudge.symbolic as sym
+ from grudge import sym
sym_x = sym.nodes(2)
sym_source_center_dist = sym_x - source_center
+ sym_sin = sym.CFunction("sin")
+ sym_exp = sym.CFunction("sin")
+ sym_t = sym.ScalarParameter("t")
from grudge.models.wave import StrongWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
- op = StrongWaveOperator(-0.1, discr.dimensions,
+ op = StrongWaveOperator(-0.1, discr.dim,
source_f=(
- sym.CFunction("sin")(source_omega*sym.ScalarParameter("t"))
- * sym.CFunction("exp")(
+ sym_sin(source_omega*sym_t)
+ * sym_exp(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2)),
dirichlet_tag=BTAG_NONE,
@@ -59,18 +63,16 @@ def main(write_output=True):
radiation_tag=BTAG_ALL,
flux_type="upwind")
+ queue = cl.CommandQueue(discr.cl_context)
from pytools.obj_array import join_fields
- fields = join_fields(discr.volume_zeros(),
- [discr.volume_zeros() for i in range(discr.dimensions)])
+ fields = join_fields(discr.zeros(queue),
+ [discr.zeros(queue) for i in range(discr.dim)])
- from leap.method.rk import LSRK4TimeStepper
- from leap.vm.codegen import PythonCodeGenerator
+ # FIXME
+ #dt = op.estimate_rk4_timestep(discr, fields=fields)
- dt_method = LSRK4TimeStepper(component_id="y")
- dt_stepper = PythonCodeGenerator.get_class(dt_method.generate())
-
- dt = op.estimate_timestep(discr, fields=fields)
- dt_stepper.set_up(t_start=0, dt_start=dt, context={"y": fields})
+ dt = 0.001
+ dt_stepper = set_up_rk4(dt, fields)
final_t = 10
nsteps = int(final_t/dt)
diff --git a/grudge/dt_finding.py b/grudge/dt_finding.py
new file mode 100644
index 0000000000000000000000000000000000000000..86e8f9a4be9210e825c25388e2c62cb3925de58d
--- /dev/null
+++ b/grudge/dt_finding.py
@@ -0,0 +1,94 @@
+"""Helpers for estimating a stable time step."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__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.
+"""
+
+from pytools import memoize_on_first_arg
+from meshmode.discretization.poly_element import PolynomialWarpAndBlendElementGroup
+import numpy.linalg as la
+
+
+class WarpAndBlendTimestepInfo(object):
+ @staticmethod
+ def dt_non_geometric_factor(discr, grp):
+ if grp.dim == 1:
+ if grp.order == 0:
+ return 1
+ else:
+ unodes = grp.unit_nodes
+ return la.norm(unodes[0] - unodes[1]) * 0.85
+ else:
+ unodes = grp.unit_nodes
+ vertex_indices = grp.vertex_indices()
+ return 2 / 3 * \
+ min(min(min(
+ la.norm(unodes[face_node_index]-unodes[vertex_index])
+ for vertex_index in vertex_indices
+ if vertex_index != face_node_index)
+ for face_node_index in face_indices)
+ for face_indices in self.face_indices())
+
+ @staticmethod
+ def dt_geometric_factor(discr, grp):
+ if grp.dim == 1:
+ return abs(el.map.jacobian())
+
+ elif grp.dim == 2:
+ area = abs(2 * el.map.jacobian())
+ semiperimeter = sum(la.norm(vertices[vi1] - vertices[vi2])
+ for vi1, vi2 in [(0, 1), (1, 2), (2, 0)])/2
+ return area / semiperimeter
+
+ elif grp.dim == 3:
+ result = abs(el.map.jacobian())/max(abs(fj) for fj in el.face_jacobians)
+ if grp.order in [1, 2]:
+ from warnings import warn
+ warn("cowardly halving timestep for order 1 and 2 tets "
+ "to avoid CFL issues")
+ result /= 2
+
+ return result
+
+ else:
+ raise NotImplementedError("cannot estimate timestep for "
+ "%d-dimensional elements" % grp.dim)
+
+
+_GROUP_TYPE_TO_INFO = {
+ PolynomialWarpAndBlendElementGroup: WarpAndBlendTimestepInfo
+ }
+
+
+@memoize_on_first_arg
+def dt_non_geometric_factor(discr):
+ return min(
+ _GROUP_TYPE_TO_INFO[type(grp)].dt_non_geometric_factor(discr, grp)
+ for grp in discr.groups)
+
+
+@memoize_on_first_arg
+def dt_geometric_factor(discr):
+ return min(
+ _GROUP_TYPE_TO_INFO[type(grp)].dt_geometric_factor(discr, grp)
+ for grp in discr.groups)
diff --git a/grudge/models/__init__.py b/grudge/models/__init__.py
index 0b33a6f3d45a6a125311c885fb2f2e2944469eb3..9ccd1ab95a4ace69555c1cb1d030023cdaf85b04 100644
--- a/grudge/models/__init__.py
+++ b/grudge/models/__init__.py
@@ -51,18 +51,13 @@ class TimeDependentOperator(Operator):
class HyperbolicOperator(Operator):
"""A base class for hyperbolic Discontinuous Galerkin operators."""
- def estimate_timestep(self, discr,
- stepper=None, stepper_class=None, stepper_args=None,
- t=None, fields=None):
- u"""Estimate the largest stable timestep, given a time stepper
- `stepper_class`. If none is given, RK4 is assumed.
+ def estimate_rk4_timestep(self, discr, t=None, fields=None):
+ u"""Estimate the largest stable timestep for an RK4 method.
"""
- rk4_dt = 1 / self.max_eigenvalue(t, fields, discr) \
- * (discr.dt_non_geometric_factor()
- * discr.dt_geometric_factor())
-
- from grudge.timestep.stability import \
- approximate_rk4_relative_imag_stability_region
- return rk4_dt * approximate_rk4_relative_imag_stability_region(
- stepper, stepper_class, stepper_args)
+ from grudge.dt_finding import (
+ dt_non_geometric_factor,
+ dt_geometric_factor)
+ return 1 / self.max_eigenvalue(t, fields, discr) \
+ * (dt_non_geometric_factor(discr)
+ * dt_geometric_factor(discr))
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 697a388cb549074cf562943425c293776cdff4ad..f821206d513d7b8f920999259031126e3c23b755 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -29,10 +29,10 @@ THE SOFTWARE.
from pytools import memoize_method
-import grudge.mesh
from grudge.models import HyperbolicOperator
from grudge.symbolic.primitives import make_common_subexpression as cse
-from grudge.tools import make_obj_array
+from meshmode.mesh import BTAG_ALL, BTAG_NONE
+from pytools.obj_array import join_fields, make_obj_array
# TODO: Check PML
@@ -126,7 +126,6 @@ class MaxwellOperator(HyperbolicOperator):
"""
from grudge.flux import (make_normal, FluxVectorPlaceholder,
FluxConstantPlaceholder)
- from grudge.tools import join_fields
normal = make_normal(self.dimensions)
@@ -205,7 +204,6 @@ class MaxwellOperator(HyperbolicOperator):
return self.space_cross_h(nabla, field)
from grudge.symbolic import make_nabla
- from grudge.tools import join_fields
nabla = make_nabla(self.dimensions)
@@ -229,7 +227,6 @@ class MaxwellOperator(HyperbolicOperator):
"Construct part of the flux operator template for PEC boundary conditions"
e, h = self.split_eh(self.field_placeholder(w))
- from grudge.tools import join_fields
from grudge.symbolic import BoundarizeOperator
pec_e = BoundarizeOperator(self.pec_tag)(e)
pec_h = BoundarizeOperator(self.pec_tag)(h)
@@ -240,7 +237,6 @@ class MaxwellOperator(HyperbolicOperator):
"Construct part of the flux operator template for PMC boundary conditions"
e, h = self.split_eh(self.field_placeholder(w))
- from grudge.tools import join_fields
from grudge.symbolic import BoundarizeOperator
pmc_e = BoundarizeOperator(self.pmc_tag)(e)
pmc_h = BoundarizeOperator(self.pmc_tag)(h)
@@ -256,7 +252,6 @@ class MaxwellOperator(HyperbolicOperator):
absorb_normal = normal(self.absorb_tag, self.dimensions)
from grudge.symbolic import BoundarizeOperator, Field
- from grudge.tools import join_fields
e, h = self.split_eh(self.field_placeholder(w))
@@ -315,7 +310,6 @@ class MaxwellOperator(HyperbolicOperator):
Combines the relevant operator templates for spatial
derivatives, flux, boundary conditions etc.
"""
- from grudge.tools import join_fields
w = self.field_placeholder(w)
if self.fixed_material:
@@ -410,7 +404,6 @@ class MaxwellOperator(HyperbolicOperator):
e = default_fld(e, e_components)
h = default_fld(h, h_components)
- from grudge.tools import join_fields
return join_fields(e, h)
assemble_fields = assemble_eh
diff --git a/grudge/models/second_order.py b/grudge/models/second_order.py
index de8404e51654054a53837d9d2336c5e3e18aeaea..24260c0f7012ab7bbd6f7adc0fc23a0248bbee32 100644
--- a/grudge/models/second_order.py
+++ b/grudge/models/second_order.py
@@ -27,14 +27,15 @@ THE SOFTWARE.
import numpy as np
-import grudge.symbolic
+import grudge.symbolic.mappers as mappers
+from grudge import sym
# {{{ stabilization term generator
-class StabilizationTermGenerator(grudge.symbolic.IdentityMapper):
+class StabilizationTermGenerator(mappers.IdentityMapper):
def __init__(self, flux_args):
- grudge.symbolic.IdentityMapper.__init__(self)
+ super(StabilizationTermGenerator, self).__init__()
self.flux_args = flux_args
self.flux_arg_lookup = dict(
(flux_arg, i) for i, flux_arg in enumerate(flux_args))
@@ -135,14 +136,14 @@ class StabilizationTermGenerator(grudge.symbolic.IdentityMapper):
# {{{ neumann bc generator
-class NeumannBCGenerator(grudge.symbolic.IdentityMapper):
+class NeumannBCGenerator(mappers.IdentityMapper):
def __init__(self, tag, bc):
- grudge.symbolic.IdentityMapper.__init__(self)
+ super(NeumannBCGenerator, self).__init__()
self.tag = tag
self.bc = bc
def map_operator_binding(self, expr):
- if isinstance(expr.op, grudge.symbolic.DiffOperatorBase):
+ if isinstance(expr.op, sym.DiffOperatorBase):
from grudge.symbolic import \
WeakFormDiffOperatorBase, \
StrongFormDiffOperatorBase
@@ -158,9 +159,9 @@ class NeumannBCGenerator(grudge.symbolic.IdentityMapper):
return (self.bc * factor *
BoundaryNormalComponent(self.tag, expr.op.xyz_axis))
- elif isinstance(expr.op, grudge.symbolic.FluxOperatorBase):
+ elif isinstance(expr.op, sym.FluxOperatorBase):
return 0
- elif isinstance(expr.op, grudge.symbolic.InverseMassOperator):
+ elif isinstance(expr.op, sym.InverseMassOperator):
return self.rec(expr.field)
else:
raise ValueError("neumann normal direction generator doesn't know "
@@ -169,9 +170,9 @@ class NeumannBCGenerator(grudge.symbolic.IdentityMapper):
# }}}
-class IPDGDerivativeGenerator(grudge.symbolic.IdentityMapper):
+class IPDGDerivativeGenerator(mappers.IdentityMapper):
def map_operator_binding(self, expr):
- if isinstance(expr.op, grudge.symbolic.DiffOperatorBase):
+ if isinstance(expr.op, sym.DiffOperatorBase):
from grudge.symbolic import (
WeakFormDiffOperatorBase,
StrongFormDiffOperatorBase)
@@ -187,14 +188,13 @@ class IPDGDerivativeGenerator(grudge.symbolic.IdentityMapper):
from grudge.symbolic import DifferentiationOperator
return factor*DifferentiationOperator(expr.op.xyz_axis)(expr.field)
- elif isinstance(expr.op, grudge.symbolic.FluxOperatorBase):
+ elif isinstance(expr.op, sym.FluxOperatorBase):
return 0
- elif isinstance(expr.op, grudge.symbolic.InverseMassOperator):
+ elif isinstance(expr.op, sym.InverseMassOperator):
return self.rec(expr.field)
elif isinstance(expr.op,
- grudge.symbolic.QuadratureInteriorFacesGridUpsampler):
- return grudge.symbolic.IdentityMapper.map_operator_binding(
- self, expr)
+ sym.QuadratureInteriorFacesGridUpsampler):
+ return super(IPDGDerivativeGenerator, self).map_operator_binding(expr)
else:
from grudge.symbolic.tools import pretty
raise ValueError("IPDG derivative generator doesn't know "
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index a3a9aaf4448127c0bf4889f92adecef821966d9f..53a0a848c25fb5c3650f9cb5d2f1f6cb525bdca5 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -29,9 +29,10 @@ THE SOFTWARE.
import numpy as np
from grudge.models import HyperbolicOperator
-from grudge.second_order import CentralSecondDerivative
+from grudge.models.second_order import CentralSecondDerivative
from meshmode.mesh import BTAG_ALL, BTAG_NONE
from grudge import sym, sym_flux
+from pytools.obj_array import join_fields
# {{{ constant-velocity
@@ -122,7 +123,6 @@ class StrongWaveOperator(HyperbolicOperator):
v = w[1:]
# boundary conditions -------------------------------------------------
- from grudge.tools import join_fields
# dirichlet BCs -------------------------------------------------------
from grudge.symbolic import normal, Field
@@ -160,7 +160,6 @@ class StrongWaveOperator(HyperbolicOperator):
nabla = make_nabla(d)
flux_op = get_flux_operator(self.flux())
- from grudge.tools import join_fields
result = (
- join_fields(
-self.c*np.dot(nabla, v),
@@ -252,7 +251,6 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
v = w[2:]
normal = make_normal(dim)
- from grudge.tools import join_fields
flux = self.time_sign*1/2*join_fields(
c.ext * np.dot(v.ext, normal)
- c.int * np.dot(v.int, normal),
@@ -299,11 +297,9 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
u = w[0]
v = w[1:]
- from grudge.tools import join_fields
flux_w = join_fields(self.c, w)
# {{{ boundary conditions
- from grudge.tools import join_fields
# Dirichlet
dir_c = BoundarizeOperator(self.dirichlet_tag) * self.c
diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py
index 9b20ee4090f17e6179d0358d91a6ae9b18594bd3..6ed54d515b7bebe3b5d34d52eaa10942e188b0bb 100644
--- a/grudge/shortcuts.py
+++ b/grudge/shortcuts.py
@@ -37,3 +37,15 @@ def make_discretization(mesh, order, **kwargs):
return Discretization(cl_ctx, mesh,
PolynomialWarpAndBlendGroupFactory(order=order))
+
+
+def set_up_rk4(dt, fields, t_start=0):
+ from leap.method.rk import LSRK4Method
+ from leap.vm.codegen import PythonCodeGenerator
+
+ dt_method = LSRK4Method(component_id="y")
+ dt_stepper = PythonCodeGenerator().get_class(dt_method.generate())
+
+ dt_stepper.set_up(t_start=t_start, dt_start=dt, context={"y": fields})
+
+ return dt_stepper
diff --git a/grudge/symbolic/flux/primitives b/grudge/symbolic/flux/primitives.py
similarity index 100%
rename from grudge/symbolic/flux/primitives
rename to grudge/symbolic/flux/primitives.py