diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index cb51964f6442192ed30ea953fe1dc1be2720bdab..bd399c3a38542280f49cb66b00d6979541085a31 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -523,7 +523,7 @@ def test_stepper_equivalence(ctx_factory, order=4):
fused_steps = fused_stepper.run(ic, t_start, dt, t_end)
- for t_ref, (u_ref, v_ref) in stepper.run(ic, t_start, dt, t_end):
+ for t_ref, (u_ref, _v_ref) in stepper.run(ic, t_start, dt, t_end):
step += 1
logger.debug("step %d/%d", step, nsteps)
t, (u, v) = next(fused_steps)
@@ -620,7 +620,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper):
expr.op, expr.field, profile_data)
else:
- assert False, ("unknown operator: %s" % expr.op)
+ raise AssertionError("unknown operator: %s" % expr.op)
else:
logger.debug("assignment not profiled: %s", expr)
@@ -811,7 +811,7 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
step = 0
nsteps = int(np.ceil((t_end + 1e-9) / dt))
- for (_, _, profile_data) in stepper.run(
+ for (_, _, profile_data) in stepper.run( # noqa: B007
ic, t_start, dt, t_end, return_profile_data=True):
step += 1
logger.info("step %d/%d", step, nsteps)
@@ -984,7 +984,7 @@ def test_stepper_timing(ctx_factory, use_fusion):
import time
t = time.time()
nsteps = int(np.ceil((t_end + 1e-9) / dt))
- for (_, _, profile_data) in stepper.run(
+ for (_, _, profile_data) in stepper.run( # noqa: B007
ic, t_start, dt, t_end, return_profile_data=True):
step += 1
tn = time.time()
@@ -1146,7 +1146,7 @@ def mem_ops_results(actx, dims):
result = {}
- for (_, _, profile_data) in stepper.run(
+ for (_, _, profile_data) in stepper.run( # noqa: B007
ic, t_start, dt, t_end, return_profile_data=True):
pass
@@ -1164,7 +1164,7 @@ def mem_ops_results(actx, dims):
result["nonfused_bytes_read_by_scalar_assignments"] \
+ result["nonfused_bytes_written_by_scalar_assignments"]
- for (_, _, profile_data) in fused_stepper.run(
+ for (_, _, profile_data) in fused_stepper.run( # noqa: B007
ic, t_start, dt, t_end, return_profile_data=True):
pass
diff --git a/grudge/execution.py b/grudge/execution.py
index 44721932737712758606ee455649c42e7c800f35..4b86fcb9e359c828a99f8f6d94a4714a999cd705 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -628,7 +628,7 @@ class BoundOperator:
else:
pass
- for key, val in context.items():
+ for val in context.values():
look_for_array_contexts(val)
if array_contexts:
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 244d324e0749ef48d19a67af9ab0be5c3156aa0c..3f4bc13d42a017fcee36ba7b05237c512ca11451 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -327,7 +327,9 @@ class MaxwellOperator(HyperbolicOperator):
1 / sym.sqrt(self.epsilon * self.mu)
)
- def max_eigenvalue(self, t, fields=None, discr=None, context={}):
+ def max_eigenvalue(self, t, fields=None, discr=None, context=None):
+ if context is None:
+ context = {}
if self.fixed_material:
return self.max_eigenvalue_expr()
else:
diff --git a/grudge/models/poisson.py b/grudge/models/poisson.py
deleted file mode 100644
index 6636b94fe4d7ea90894bf92a86d30189fe93348e..0000000000000000000000000000000000000000
--- a/grudge/models/poisson.py
+++ /dev/null
@@ -1,261 +0,0 @@
-"""Operators for Poisson problems."""
-
-__copyright__ = "Copyright (C) 2007 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.
-"""
-
-import numpy as np
-
-from grudge.models import Operator
-from grudge.second_order import LDGSecondDerivative
-import grudge.data
-import grudge.iterative
-
-
-class LaplacianOperatorBase:
- def sym_operator(self, apply_minv, u=None, dir_bc=None, neu_bc=None):
- """
- :param apply_minv: :class:`bool` specifying whether to compute a complete
- divergence operator. If False, the final application of the inverse
- mass operator is skipped. This is used in :meth:`op` in order to
- reduce the scheme :math:`M^{-1} S u = f` to :math:`S u = M f`, so
- that the mass operator only needs to be applied once, when preparing
- the right hand side in :meth:`prepare_rhs`.
-
- :class:`grudge.models.diffusion.DiffusionOperator` needs this.
- """
-
- from grudge.symbolic import Field, make_sym_vector
- from grudge.second_order import SecondDerivativeTarget
-
- if u is None:
- u = Field("u")
- if dir_bc is None:
- dir_bc = Field("dir_bc")
- if neu_bc is None:
- neu_bc = Field("neu_bc")
-
- # strong_form here allows IPDG to reuse the value of grad u.
- grad_tgt = SecondDerivativeTarget(
- self.dimensions, strong_form=True,
- operand=u)
-
- def grad_bc_getter(tag, expr):
- assert tag == self.dirichlet_tag
- return dir_bc
- self.scheme.grad(grad_tgt,
- bc_getter=grad_bc_getter,
- dirichlet_tags=[self.dirichlet_tag],
- neumann_tags=[self.neumann_tag])
-
- def apply_diff_tensor(v):
- if isinstance(self.diffusion_tensor, np.ndarray):
- sym_diff_tensor = self.diffusion_tensor
- else:
- sym_diff_tensor = (make_sym_vector(
- "diffusion", self.dimensions**2)
- .reshape(self.dimensions, self.dimensions))
-
- return np.dot(sym_diff_tensor, v)
-
- div_tgt = SecondDerivativeTarget(
- self.dimensions, strong_form=False,
- operand=apply_diff_tensor(grad_tgt.minv_all))
-
- def div_bc_getter(tag, expr):
- if tag == self.dirichlet_tag:
- return dir_bc
- elif tag == self.neumann_tag:
- return neu_bc
- else:
- assert False, "divergence bc getter " \
- "asked for '%s' BC for '%s'" % (tag, expr)
-
- self.scheme.div(div_tgt,
- div_bc_getter,
- dirichlet_tags=[self.dirichlet_tag],
- neumann_tags=[self.neumann_tag])
-
- if apply_minv:
- return div_tgt.minv_all
- else:
- return div_tgt.all
-
-
-class PoissonOperator(Operator, LaplacianOperatorBase):
- """Implements the Local Discontinuous Galerkin (LDG) Method for elliptic
- operators.
-
- See P. Castillo et al.,
- Local discontinuous Galerkin methods for elliptic problems",
- Communications in Numerical Methods in Engineering 18, no. 1 (2002): 69-75.
- """
-
- def __init__(self, dimensions, diffusion_tensor=None,
- dirichlet_bc=grudge.data.ConstantGivenFunction(),
- dirichlet_tag="dirichlet",
- neumann_bc=grudge.data.ConstantGivenFunction(),
- neumann_tag="neumann",
- scheme=LDGSecondDerivative()):
- self.dimensions = dimensions
-
- self.scheme = scheme
-
- self.dirichlet_bc = dirichlet_bc
- self.dirichlet_tag = dirichlet_tag
- self.neumann_bc = neumann_bc
- self.neumann_tag = neumann_tag
-
- if diffusion_tensor is None:
- diffusion_tensor = np.eye(dimensions)
- self.diffusion_tensor = diffusion_tensor
-
- # bound operator ----------------------------------------------------------
- def bind(self, discr):
- """Return a :class:`BoundPoissonOperator`."""
-
- assert self.dimensions == discr.dimensions
-
- from grudge.mesh import check_bc_coverage
- check_bc_coverage(discr.mesh, [self.dirichlet_tag, self.neumann_tag])
-
- return BoundPoissonOperator(self, discr)
-
-
-class BoundPoissonOperator(grudge.iterative.OperatorBase):
- """Returned by :meth:`PoissonOperator.bind`."""
-
- def __init__(self, poisson_op, discr):
- grudge.iterative.OperatorBase.__init__(self)
- self.discr = discr
-
- pop = self.poisson_op = poisson_op
-
- op = pop.sym_operator(
- apply_minv=False, dir_bc=0, neu_bc=0)
- bc_op = pop.sym_operator(apply_minv=False)
-
- self.compiled_op = discr.compile(op)
- self.compiled_bc_op = discr.compile(bc_op)
-
- if not isinstance(pop.diffusion_tensor, np.ndarray):
- self.diffusion = pop.diffusion_tensor.volume_interpolant(discr)
-
- # Check whether use of Poincaré mean-value method is required.
- # (for pure Neumann or pure periodic)
-
- from grudge.mesh import BTAG_ALL
- self.poincare_mean_value_hack = (
- len(self.discr.get_boundary(BTAG_ALL).nodes)
- == len(self.discr.get_boundary(poisson_op.neumann_tag).nodes))
-
- @property
- def dtype(self):
- return self.discr.default_scalar_type
-
- @property
- def shape(self):
- nodes = len(self.discr)
- return nodes, nodes
-
- def op(self, u):
- context = {"u": u}
- if not isinstance(self.poisson_op.diffusion_tensor, np.ndarray):
- context["diffusion"] = self.diffusion
-
- result = self.compiled_op(**context)
-
- if self.poincare_mean_value_hack:
- state_int = self.discr.integral(u)
- mean_state = state_int / self.discr.mesh_volume()
- return result - mean_state * self.discr._mass_ones()
- else:
- return result
-
- __call__ = op
-
- def prepare_rhs(self, rhs):
- """Prepare the right-hand side for the linear system op(u)=rhs(f).
-
- In matrix form, LDG looks like this:
-
- .. math::
- Mv = Cu + g
- Mf = Av + Bu + h
-
- where v is the auxiliary vector, u is the argument of the operator, f
- is the result of the grad operator, g and h are inhom boundary data, and
- A,B,C are some operator+lifting matrices.
-
- .. math::
-
- M f = A M^{-1}(Cu + g) + Bu + h
-
- so the linear system looks like
-
- .. math::
-
- M f = A M^{-1} Cu + A M^{-1} g + Bu + h
- M f - A M^{-1} g - h = (A M^{-1} C + B)u (*)
-
- So the right hand side we're putting together here is really
-
- .. math::
-
- M f - A M^{-1} g - h
-
- .. note::
-
- Resist the temptation to left-multiply by the inverse
- mass matrix, as this will result in a non-symmetric
- matrix, which (e.g.) a conjugate gradient Krylov
- solver will not take well.
- """
- pop = self.poisson_op
-
- from grudge.symbolic import MassOperator
- return (MassOperator().apply(self.discr, rhs)
- - self.compiled_bc_op(
- u=self.discr.volume_zeros(),
- dir_bc=pop.dirichlet_bc.boundary_interpolant(
- self.discr, pop.dirichlet_tag),
- neu_bc=pop.neumann_bc.boundary_interpolant(
- self.discr, pop.neumann_tag)))
-
-
-class HelmholtzOperator(PoissonOperator):
- def __init__(self, k, *args, **kwargs):
- PoissonOperator.__init__(self, *args, **kwargs)
- self.k = k
-
- def sym_operator(self, apply_minv, u=None, dir_bc=None, neu_bc=None):
- from grudge.symbolic import Field
- if u is None:
- u = Field("u")
-
- result = PoissonOperator.sym_operator(self,
- apply_minv, u, dir_bc, neu_bc)
-
- if apply_minv:
- return result + self.k**2 * u
- else:
- from grudge.symbolic import MassOperator
- return result + self.k**2 * MassOperator()(u)
diff --git a/grudge/models/second_order.py b/grudge/models/second_order.py
deleted file mode 100644
index 8f7251d329ea1cbbb960e368fe3c272f481ca36a..0000000000000000000000000000000000000000
--- a/grudge/models/second_order.py
+++ /dev/null
@@ -1,505 +0,0 @@
-"""Schemes for second-order derivatives."""
-
-__copyright__ = "Copyright (C) 2009 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.
-"""
-
-import numpy as np
-import grudge.symbolic.mappers as mappers
-from grudge import sym
-
-
-# {{{ stabilization term generator
-
-class StabilizationTermGenerator(mappers.IdentityMapper):
- def __init__(self, flux_args):
- super().__init__()
- self.flux_args = flux_args
- self.flux_arg_lookup = {
- flux_arg: i for i, flux_arg in enumerate(flux_args)}
-
- def get_flux_arg_idx(self, expr, quad_above):
- from grudge.symbolic.mappers import QuadratureDetector
-
- quad_below = QuadratureDetector()(expr)
- if quad_above:
- if quad_below is not None:
- # Both the part of the expression above and below the
- # differentiation operator had quadrature upsamplers in it.
- # Since we're removing the differentiation operator, there are
- # now two layers of quadrature operators. We need to change the
- # inner layer to be the only layer.
-
- from grudge.symbolic.mappers import QuadratureUpsamplerChanger
- expr = QuadratureUpsamplerChanger(quad_above[0])(expr)
- else:
- # Only the part of the expression above the differentiation
- # operator had quadrature. Insert quadrature here, be done.
- expr = quad_above[0](expr)
- else:
- if quad_below is not None:
- # Only the part of the expression below the differentiation
- # operator had quadrature--the stuff above doesn't want it.
- # Get rid of it.
- from grudge.symbolic.mappers import QuadratureUpsamplerRemover
- expr = QuadratureUpsamplerRemover({}, do_warn=False)(expr)
- else:
- # No quadrature, no headaches.
- pass
-
- try:
- return self.flux_arg_lookup[expr]
- except KeyError:
- flux_arg_idx = len(self.flux_args)
- self.flux_arg_lookup[expr] = flux_arg_idx
- self.flux_args.append(expr)
- return flux_arg_idx
-
- def map_operator_binding(self, expr, quad_above=[]):
- from grudge.symbolic.operators import (
- DiffOperatorBase, FluxOperatorBase,
- InverseMassOperator,
- QuadratureInteriorFacesGridUpsampler)
-
- if isinstance(expr.op, DiffOperatorBase):
- flux_arg_idx = self.get_flux_arg_idx(expr.field, quad_above=quad_above)
-
- from grudge.symbolic import \
- WeakFormDiffOperatorBase, \
- StrongFormDiffOperatorBase
- if isinstance(expr.op, WeakFormDiffOperatorBase):
- factor = -1
- elif isinstance(expr.op, StrongFormDiffOperatorBase):
- factor = 1
- else:
- raise RuntimeError("unknown type of differentiation "
- "operator encountered by stab term generator")
-
- from grudge.flux import Normal, FluxScalarPlaceholder
- sph = FluxScalarPlaceholder(flux_arg_idx)
- return (factor
- * Normal(expr.op.xyz_axis)
- * (sph.int - sph.ext))
-
- elif isinstance(expr.op, FluxOperatorBase):
- return 0
- elif isinstance(expr.op, InverseMassOperator):
- return self.rec(expr.field, quad_above)
- elif isinstance(expr.op, QuadratureInteriorFacesGridUpsampler):
- if quad_above:
- raise RuntimeError("double quadrature upsampler found "
- "when generating stabilization term")
- return self.rec(expr.field, [expr.op])
- else:
- from grudge.symbolic.tools import pretty
- raise ValueError("stabilization term generator doesn't know "
- "what to do with '%s'" % pretty(expr))
-
- def map_variable(self, expr, quad_above=[]):
- from grudge.flux import FieldComponent
- return FieldComponent(
- self.get_flux_arg_idx(expr, quad_above),
- is_interior=True)
-
- def map_subscript(self, expr, quad_above=[]):
- from pymbolic.primitives import Variable
- assert isinstance(expr.aggregate, Variable)
- from grudge.flux import FieldComponent
- return FieldComponent(
- self.get_flux_arg_idx(expr, quad_above),
- is_interior=True)
-
-# }}}
-
-
-# {{{ neumann bc generator
-
-class NeumannBCGenerator(mappers.IdentityMapper):
- def __init__(self, tag, bc):
- super().__init__()
- self.tag = tag
- self.bc = bc
-
- def map_operator_binding(self, expr):
- if isinstance(expr.op, sym.DiffOperatorBase):
- from grudge.symbolic import \
- WeakFormDiffOperatorBase, \
- StrongFormDiffOperatorBase
- if isinstance(expr.op, WeakFormDiffOperatorBase):
- factor = -1
- elif isinstance(expr.op, StrongFormDiffOperatorBase):
- factor = 1
- else:
- raise RuntimeError("unknown type of differentiation "
- "operator encountered by stab term generator")
-
- from grudge.symbolic import BoundaryNormalComponent
- return (self.bc * factor
- * BoundaryNormalComponent(self.tag, expr.op.xyz_axis))
-
- elif isinstance(expr.op, sym.FluxOperatorBase):
- return 0
- elif isinstance(expr.op, sym.InverseMassOperator):
- return self.rec(expr.field)
- else:
- raise ValueError("neumann normal direction generator doesn't know "
- "what to do with '%s'" % expr)
-
-# }}}
-
-
-class IPDGDerivativeGenerator(mappers.IdentityMapper):
- def map_operator_binding(self, expr):
- if isinstance(expr.op, sym.DiffOperatorBase):
- from grudge.symbolic import (
- WeakFormDiffOperatorBase,
- StrongFormDiffOperatorBase)
-
- if isinstance(expr.op, WeakFormDiffOperatorBase):
- factor = -1
- elif isinstance(expr.op, StrongFormDiffOperatorBase):
- factor = 1
- else:
- raise RuntimeError("unknown type of differentiation "
- "operator encountered by stab term generator")
-
- from grudge.symbolic import DifferentiationOperator
- return factor*DifferentiationOperator(expr.op.xyz_axis)(expr.field)
-
- elif isinstance(expr.op, sym.FluxOperatorBase):
- return 0
- elif isinstance(expr.op, sym.InverseMassOperator):
- return self.rec(expr.field)
- elif isinstance(expr.op,
- sym.QuadratureInteriorFacesGridUpsampler):
- return super().map_operator_binding(expr)
- else:
- from grudge.symbolic.tools import pretty
- raise ValueError("IPDG derivative generator doesn't know "
- "what to do with '%s'" % pretty(expr))
-
-
-# {{{ second derivative target
-
-class SecondDerivativeTarget:
- def __init__(self, dimensions, strong_form, operand,
- int_flux_operand=None,
- bdry_flux_int_operand=None):
- """
- :param int_flux_operand: if not None, is used as the interior
- argument to the interior fluxes. This is useful e.g. if the boundary
- values are on a quadrature grid--in this case, *bdry_flux_int_operand*
- can be passed to also be on a boundary grid. If it is None, it defaults
- to *operand*.
-
- :param bdry_flux_int_operand: if not None, is used as the interior
- argument to the boundary fluxes. This is useful e.g. if the boundary
- values are on a quadrature grid--in this case, *bdry_flux_int_operand*
- can be passed to also be on a boundary grid. If it is None, it defaults
- to *int_flux_operand*.
- """
- self.dimensions = dimensions
- self.operand = operand
-
- if int_flux_operand is None:
- int_flux_operand = operand
- if bdry_flux_int_operand is None:
- bdry_flux_int_operand = int_flux_operand
-
- self.int_flux_operand = int_flux_operand
- self.bdry_flux_int_operand = bdry_flux_int_operand
-
- self.strong_form = strong_form
- if strong_form:
- self.strong_neg = -1
- else:
- self.strong_neg = 1
-
- self.local_derivatives = 0
- self.inner_fluxes = 0
- self.boundary_fluxes = 0
-
- def add_local_derivatives(self, expr):
- self.local_derivatives = self.local_derivatives \
- + (-self.strong_neg)*expr
-
- def vec_times(self, vec, operand):
- from pytools.obj_array import is_obj_array
-
- if is_obj_array(operand):
- if len(operand) != self.dimensions:
- raise ValueError("operand of vec_times must have %d dimensions"
- % self.dimensions)
-
- return np.dot(vec, operand)
- else:
- return vec*operand
-
- def normal_times_flux(self, flux):
- from grudge.flux import make_normal
- return self.vec_times(make_normal(self.dimensions), flux)
-
- def apply_diff(self, nabla, operand):
- from pytools.obj_array import make_obj_array, is_obj_array
- if is_obj_array(operand):
- if len(operand) != self.dimensions:
- raise ValueError("operand of apply_diff must have %d dimensions"
- % self.dimensions)
-
- return sum(nabla[i](operand[i]) for i in range(self.dimensions))
- else:
- return make_obj_array(
- [nabla[i](operand) for i in range(self.dimensions)])
-
- def _local_nabla(self):
- if self.strong_form:
- from grudge.symbolic import make_stiffness
- return make_stiffness(self.dimensions)
- else:
- from grudge.symbolic import make_stiffness_t
- return make_stiffness_t(self.dimensions)
-
- def add_derivative(self, operand=None):
- if operand is None:
- operand = self.operand
-
- self.add_local_derivatives(
- self.apply_diff(self._local_nabla(), operand))
-
- def add_inner_fluxes(self, flux, expr):
- from grudge.symbolic import get_flux_operator
- self.inner_fluxes = self.inner_fluxes \
- + get_flux_operator(self.strong_neg*flux)(expr)
-
- def add_boundary_flux(self, flux, volume_expr, bdry_expr, tag):
- from grudge.symbolic import BoundaryPair, get_flux_operator
- self.boundary_fluxes = self.boundary_fluxes + \
- get_flux_operator(self.strong_neg*flux)(BoundaryPair(
- volume_expr, bdry_expr, tag))
-
- @property
- def fluxes(self):
- return self.inner_fluxes + self.boundary_fluxes
-
- @property
- def all(self):
- return self.local_derivatives + self.fluxes
-
- @property
- def minv_all(self):
- from grudge.symbolic.primitives import make_common_subexpression as cse
- from grudge.symbolic.operators import InverseMassOperator
- return (cse(InverseMassOperator()(self.local_derivatives), "grad_loc")
- + cse(InverseMassOperator()(self.fluxes), "grad_flux"))
-
-# }}}
-
-
-# {{{ second derivative schemes
-
-class SecondDerivativeBase:
- def grad(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
- """
- :param bc_getter: a function (tag, volume_expr) -> boundary expr.
- *volume_expr* will be None to query the Neumann condition.
- """
- n_times = tgt.normal_times_flux
-
- if tgt.strong_form:
- def adjust_flux(f):
- return n_times(u.int) - f
- else:
- def adjust_flux(f):
- return f
-
- from grudge.flux import FluxScalarPlaceholder
- u = FluxScalarPlaceholder()
-
- tgt.add_derivative()
- tgt.add_inner_fluxes(
- adjust_flux(self.grad_interior_flux(tgt, u)),
- tgt.int_flux_operand)
-
- for tag in dirichlet_tags:
- tgt.add_boundary_flux(
- adjust_flux(n_times(u.ext)),
- tgt.bdry_flux_int_operand, bc_getter(tag, tgt.operand), tag)
-
- for tag in neumann_tags:
- tgt.add_boundary_flux(
- adjust_flux(n_times(u.int)),
- tgt.bdry_flux_int_operand, 0, tag)
-
- def add_div_bcs(self, tgt, bc_getter, dirichlet_tags, neumann_tags,
- stab_term, adjust_flux, flux_v, flux_arg_int,
- grad_flux_arg_count):
- from pytools.obj_array import make_obj_array, join_fields
- n_times = tgt.normal_times_flux
-
- def unwrap_cse(expr):
- from pymbolic.primitives import CommonSubexpression
- if isinstance(expr, CommonSubexpression):
- return expr.child
- else:
- return expr
-
- for tag in dirichlet_tags:
- dir_bc_w = join_fields(
- [0]*grad_flux_arg_count,
- [bc_getter(tag, unwrap_cse(vol_expr)) for vol_expr in
- flux_arg_int[grad_flux_arg_count:]])
- tgt.add_boundary_flux(
- adjust_flux(n_times(flux_v.int-stab_term)),
- flux_arg_int, dir_bc_w, tag)
-
- loc_bc_vec = make_obj_array([0]*len(flux_arg_int))
-
- for tag in neumann_tags:
- neu_bc_w = join_fields(
- NeumannBCGenerator(tag, bc_getter(tag, None))(tgt.operand),
- [0]*len(flux_arg_int))
-
- tgt.add_boundary_flux(
- adjust_flux(n_times(flux_v.ext)),
- loc_bc_vec, neu_bc_w, tag)
-
-
-class LDGSecondDerivative(SecondDerivativeBase):
- def __init__(self, beta_value=0.5, stab_coefficient=1):
- self.beta_value = beta_value
- self.stab_coefficient = stab_coefficient
-
- def beta(self, tgt):
- return np.array([self.beta_value]*tgt.dimensions, dtype=np.float64)
-
- def grad_interior_flux(self, tgt, u):
- from grudge.symbolic.primitives import make_common_subexpression as cse
- n_times = tgt.normal_times_flux
- v_times = tgt.vec_times
-
- return n_times(
- cse(u.avg, "u_avg")
- - v_times(self.beta(tgt), n_times(u.int-u.ext)))
-
- def div(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
- """
- :param bc_getter: a function (tag, volume_expr) -> boundary expr.
- *volume_expr* will be None to query the Neumann condition.
- """
-
- from grudge.symbolic.primitives import make_common_subexpression as cse
- from grudge.flux import FluxVectorPlaceholder, PenaltyTerm
-
- n_times = tgt.normal_times_flux
- v_times = tgt.vec_times
-
- if tgt.strong_form:
- def adjust_flux(f):
- return n_times(flux_v.int) - f
- else:
- def adjust_flux(f):
- return f
-
- flux_v = FluxVectorPlaceholder(tgt.dimensions)
-
- stab_term_generator = StabilizationTermGenerator(
- list(tgt.int_flux_operand))
- stab_term = (self.stab_coefficient * PenaltyTerm()
- * stab_term_generator(tgt.int_flux_operand))
-
- flux = n_times(flux_v.avg
- + v_times(self.beta(tgt),
- cse(n_times(flux_v.int - flux_v.ext), "jump_v"))
- - stab_term)
-
- from pytools.obj_array import make_obj_array
- flux_arg_int = cse(make_obj_array(stab_term_generator.flux_args))
-
- tgt.add_derivative(cse(tgt.operand))
- tgt.add_inner_fluxes(adjust_flux(flux), flux_arg_int)
-
- self.add_div_bcs(tgt, bc_getter, dirichlet_tags, neumann_tags,
- stab_term, adjust_flux, flux_v, flux_arg_int, tgt.dimensions)
-
-
-class StabilizedCentralSecondDerivative(LDGSecondDerivative):
- def __init__(self, stab_coefficient=1):
- LDGSecondDerivative.__init__(self, 0, stab_coefficient=stab_coefficient)
-
-
-class CentralSecondDerivative(LDGSecondDerivative):
- def __init__(self):
- LDGSecondDerivative.__init__(self, 0, 0)
-
-
-class IPDGSecondDerivative(SecondDerivativeBase):
- def __init__(self, stab_coefficient=1):
- self.stab_coefficient = stab_coefficient
-
- def grad_interior_flux(self, tgt, u):
- from grudge.symbolic.primitives import make_common_subexpression as cse
- n_times = tgt.normal_times_flux
- return n_times(cse(u.avg, "u_avg"))
-
- def div(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
- """
- :param bc_getter: a function (tag, volume_expr) -> boundary expr.
- *volume_expr* will be None to query the Neumann condition.
- """
-
- from grudge.symbolic.primitives import make_common_subexpression as cse
- from grudge.flux import FluxVectorPlaceholder, PenaltyTerm
-
- n_times = tgt.normal_times_flux
-
- if tgt.strong_form:
- def adjust_flux(f):
- return n_times(flux_v.int) - f
- else:
- def adjust_flux(f):
- return f
-
- dim = tgt.dimensions
-
- flux_w = FluxVectorPlaceholder(2*tgt.dimensions)
- flux_v = flux_w[:dim]
- pure_diff_v = flux_w[dim:]
- flux_args = (
- list(tgt.int_flux_operand)
- + list(IPDGDerivativeGenerator()(tgt.int_flux_operand)))
-
- stab_term_generator = StabilizationTermGenerator(flux_args)
- stab_term = (self.stab_coefficient * PenaltyTerm()
- * stab_term_generator(tgt.int_flux_operand))
- flux = n_times(pure_diff_v.avg - stab_term)
-
- from pytools.obj_array import make_obj_array
- flux_arg_int = cse(make_obj_array(stab_term_generator.flux_args))
-
- tgt.add_derivative(cse(tgt.operand))
- tgt.add_inner_fluxes(adjust_flux(flux), flux_arg_int)
-
- self.add_div_bcs(tgt, bc_getter, dirichlet_tags, neumann_tags,
- stab_term, adjust_flux, flux_v, flux_arg_int, 2*tgt.dimensions)
-
-# }}}
-
-# vim: fdm=marker
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 742483aed9e02264c5ac65cfb40ee8ef0dbb8ee5..f5f7c66499b2adeadc470b9be416a06d5074048a 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -796,7 +796,7 @@ def aggregate_assignments(inf_mapper, instructions, result,
schedulable.sort()
if schedulable:
- for key, name, expr in schedulable:
+ for _key, name, expr in schedulable:
ordered_names_exprs.append((name, expr))
available_names.add(name)
else:
diff --git a/setup.cfg b/setup.cfg
index e8791cc48f86a67e726046a2f7510f32a39a98a5..38dcd9de2c01e21191ca1612388f36ecbc251b70 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -9,11 +9,12 @@ exclude=
grudge/models/nd_calculus.py,
grudge/dt_finding.py
-
inline-quotes = "
docstring-quotes = """
multiline-quotes = """
+# enable-flake8-bugbear
+
[mypy]
ignore_missing_imports = True