From 41fa65dfe555581844fb47a5d60395c2c42b4134 Mon Sep 17 00:00:00 2001
From: Bogdan Enache <enache2@illinois.edu>
Date: Mon, 27 Feb 2017 10:26:20 -0600
Subject: [PATCH] Made variable-velocity work
---
examples/advection/var-velocity.py | 348 ++++++++++-------------------
examples/advection/weak.py | 119 ++++++++++
grudge/models/advection.py | 236 ++++---------------
3 files changed, 278 insertions(+), 425 deletions(-)
create mode 100644 examples/advection/weak.py
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index 83502851..e022c470 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -1,5 +1,5 @@
# grudge - the Hybrid'n'Easy DG Environment
-# Copyright (C) 2009 Andreas Stock
+# 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
@@ -15,246 +15,128 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
+import numpy as np
+import math
+import pyopencl as cl # noqa
+import pyopencl.array # noqa
+import pyopencl.clmath # noqa
+import pytest # noqa
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-
-
-
-def main(write_output=True, flux_type_arg="central", use_quadrature=True,
- final_time=20):
- from math import sin, cos, pi, sqrt
-
- from grudge.backends import guess_run_context
- rcon = guess_run_context()
-
- # mesh setup --------------------------------------------------------------
- if rcon.is_head_rank:
- #from grudge.mesh.generator import make_disk_mesh
- #mesh = make_disk_mesh()
- from grudge.mesh.generator import make_rect_mesh
- mesh = make_rect_mesh(a=(-1,-1),b=(1,1),max_area=0.008)
-
- if rcon.is_head_rank:
- mesh_data = rcon.distribute_mesh(mesh)
- else:
- mesh_data = rcon.receive_mesh()
-
- # space-time-dependent-velocity-field -------------------------------------
- # simple vortex
- class TimeDependentVField:
- """ `TimeDependentVField` is a callable expecting `(x, t)` representing space and time
-
- `x` is of the length of the spatial dimension and `t` is the time."""
- shape = (2,)
-
- def __call__(self, pt, el, t):
- x, y = pt
- # Correction-Factor to make the speed zero on the on the boundary
- #fac = (1-x**2)*(1-y**2)
- fac = 1.
- return numpy.array([-y*fac, x*fac]) * cos(pi*t)
-
- class VField:
- """ `VField` is a callable expecting `(x)` representing space
-
- `x` is of the length of the spatial dimension."""
- shape = (2,)
-
- def __call__(self, pt, el):
- x, y = pt
- # Correction-Factor to make the speed zero on the on the boundary
- #fac = (1-x**2)*(1-y**2)
- fac = 1.
- return numpy.array([-y*fac, x*fac])
-
- # space-time-dependent State BC (optional)-----------------------------------
- class TimeDependentBc_u:
- """ space and time dependent BC for state u"""
- def __call__(self, pt, el, t):
- x, y = pt
- if t <= 0.5:
- if x > 0:
- return 1
- else:
- return 0
- else:
- return 0
-
- class Bc_u:
- """ Only space dependent BC for state u"""
- def __call__(seld, pt, el):
- x, y = pt
- if x > 0:
- return 1
- else:
- return 0
-
-
- # operator setup ----------------------------------------------------------
- # In the operator setup it is possible to switch between a only space
- # dependent velocity field `VField` or a time and space dependent
- # `TimeDependentVField`.
- # For `TimeDependentVField`: advec_v=TimeDependentGivenFunction(VField())
- # For `VField`: advec_v=TimeConstantGivenFunction(GivenFunction(VField()))
- # Same for the Bc_u Function! If you don't define Bc_u then the BC for u = 0.
-
- from grudge.data import \
- ConstantGivenFunction, \
- TimeConstantGivenFunction, \
- TimeDependentGivenFunction, \
- GivenFunction
- from grudge.models.advection import VariableCoefficientAdvectionOperator
- op = VariableCoefficientAdvectionOperator(mesh.dimensions,
- #advec_v=TimeDependentGivenFunction(
- # TimeDependentVField()),
- advec_v=TimeConstantGivenFunction(
- GivenFunction(VField())),
- #bc_u_f=TimeDependentGivenFunction(
- # TimeDependentBc_u()),
- bc_u_f=TimeConstantGivenFunction(
- GivenFunction(Bc_u())),
- flux_type=flux_type_arg)
-
- # discretization setup ----------------------------------------------------
- order = 5
- if use_quadrature:
- quad_min_degrees = {"quad": 3*order}
- else:
- quad_min_degrees = {}
-
- discr = rcon.make_discretization(mesh_data, order=order,
- default_scalar_type=numpy.float64,
- debug=["cuda_no_plan"],
- quad_min_degrees=quad_min_degrees,
- tune_for=op.sym_operator(),
-
- )
- vis_discr = discr
-
- # visualization setup -----------------------------------------------------
- from grudge.visualization import VtkVisualizer
- if write_output:
- vis = VtkVisualizer(vis_discr, rcon, "fld")
-
- # initial condition -------------------------------------------------------
- if True:
- def initial(pt, el):
- # Gauss pulse
- from math import exp
- x = (pt-numpy.array([0.3, 0.5]))*8
- return exp(-numpy.dot(x, x))
- else:
- def initial(pt, el):
- # Rectangle
- x, y = pt
- if abs(x) < 0.5 and abs(y) < 0.2:
- return 2
- else:
- return 1
-
- u = discr.interpolate_volume_function(initial)
-
- # timestep setup ----------------------------------------------------------
- from grudge.timestep.runge_kutta import LSRK4TimeStepper
- stepper = LSRK4TimeStepper(
- vector_primitive_factory=discr.get_vector_primitive_factory())
-
- if rcon.is_head_rank:
- print("%d elements" % len(discr.mesh.elements))
-
- # filter setup-------------------------------------------------------------
- from grudge.discretization import ExponentialFilterResponseFunction
- from grudge.symbolic.operators import FilterOperator
- mode_filter = FilterOperator(
- ExponentialFilterResponseFunction(min_amplification=0.9,order=4))\
- .bind(discr)
-
- # diagnostics setup -------------------------------------------------------
- from pytools.log import LogManager, \
- add_general_quantities, \
- add_simulation_quantities, \
- add_run_info
-
- if write_output:
- log_file_name = "space-dep.dat"
- else:
- log_file_name = None
-
- logmgr = LogManager(log_file_name, "w", rcon.communicator)
- add_run_info(logmgr)
- add_general_quantities(logmgr)
- add_simulation_quantities(logmgr)
- discr.add_instrumentation(logmgr)
-
- stepper.add_instrumentation(logmgr)
-
- from grudge.log import Integral, LpNorm
- u_getter = lambda: u
- logmgr.add_quantity(Integral(u_getter, discr, name="int_u"))
- logmgr.add_quantity(LpNorm(u_getter, discr, p=1, name="l1_u"))
- logmgr.add_quantity(LpNorm(u_getter, discr, name="l2_u"))
-
- logmgr.add_watches(["step.max", "t_sim.max", "l2_u", "t_step.max"])
-
- # Initialize v for data output:
- v = op.advec_v.volume_interpolant(0, discr)
-
- # timestep loop -----------------------------------------------------------
- rhs = op.bind(discr)
- try:
- from grudge.timestep import times_and_steps
- step_it = times_and_steps(
- final_time=final_time, logmgr=logmgr,
- max_dt_getter=lambda t: op.estimate_timestep(discr,
- stepper=stepper, t=t, fields=u))
-
- for step, t, dt in step_it:
- if step % 10 == 0 and write_output:
- visf = vis.make_file("fld-%04d" % step)
- vis.add_data(visf, [
- ("u", discr.convert_volume(u, kind="numpy")),
- ("v", discr.convert_volume(v, kind="numpy"))
- ], time=t, step=step)
- visf.close()
-
- u = stepper(u, t, dt, rhs)
-
- # We're feeding in a discontinuity through the BCs.
- # Quadrature does not help with shock capturing--
- # therefore we do need to filter here, regardless
- # of whether quadrature is enabled.
- u = mode_filter(u)
-
- assert discr.norm(u) < 10
-
- finally:
- if write_output:
- vis.close()
-
- logmgr.close()
- discr.close()
+from pyopencl.tools import ( # noqa
+ pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+import logging
+logger = logging.getLogger(__name__)
+from grudge import sym, bind, Discretization
-if __name__ == "__main__":
- main()
+import numpy.linalg as la
+
+from pytools.obj_array import join_fields
+
+
+
+
+
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dim = 2
+
+
+ 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)
+
+ dt_factor = 4
+ h = 1/20
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+
+ c = np.array([0.1,0.1])
+ norm_c = la.norm(c)
+
+ sym_x = sym.nodes(2)
+
+ #advec_v = sym.If(sym.Comparison(
+ #np.dot(sym_x, sym_x), "<", 0.4**2),
+ #0.4, 0.1)
+
+ advec_v = sym_v = join_fields(-1*sym_x[1],sym_x[0] ) / 2
+ omega = 1/2
+
+
+ flux_type = "central"
+
+ source_center = np.array([0.1, 0.1])
+ source_width = 0.05
+ source_omega = 3
+
+ sym_x = sym.nodes(2)
+ sym_source_center_dist = sym_x - source_center
+
+
+ def f(x):
+ return sym.exp(
+ -np.dot(sym_source_center_dist, sym_source_center_dist)
+ / source_width**2)
+ def u_analytic(x):
+ return 0
+ from grudge.models.advection import VariableCoefficientAdvectionOperator, WeakAdvectionOperator
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+ op = VariableCoefficientAdvectionOperator(2,advec_v,
+ inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
+ flux_type=flux_type)
-# entry points for py.test ----------------------------------------------------
-def test_var_velocity_advection():
- from pytools.test import mark_test
- mark_long = mark_test.long
+ bound_op = bind(discr, op.sym_operator())
+
+ u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
+
+ 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
+ dt = final_time/nsteps + 1e-15
+
+
+ from grudge.shortcuts import set_up_rk4
+ dt_stepper = set_up_rk4("u", dt, u, rhs)
+
+ last_u = None
+
+ from grudge.shortcuts import make_visualizer
+ vis = make_visualizer(discr, vis_order=order)
+
+ step = 0
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ for event in dt_stepper.run(t_end=final_time):
+ if isinstance(event, dt_stepper.StateComputed):
+
+ step += 1
+
+ #print(step, event.t, norm(queue, u=event.state_component[0]))
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [ ("u", event.state_component) ])
+
+
+
+
+
+
+
+if __name__ == "__main__":
+ main()
- for flux_type in ["upwind", "central", "lf"]:
- for use_quadrature in [False, True]:
- descr = "variable-velocity-advection with %s flux" % flux_type
- if use_quadrature:
- descr += " and quadrature"
- yield descr, mark_long(main), False, flux_type, use_quadrature, 1
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
new file mode 100644
index 00000000..de4bc4d3
--- /dev/null
+++ b/examples/advection/weak.py
@@ -0,0 +1,119 @@
+# grudge - the Hybrid'n'Easy DG Environment
+# 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/>.
+
+
+import numpy as np
+import pyopencl as cl # noqa
+import pyopencl.array # noqa
+import pyopencl.clmath # noqa
+
+import pytest # noqa
+
+from pyopencl.tools import ( # noqa
+ pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, Discretization
+
+import numpy.linalg as la
+
+
+
+
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dim = 2
+
+
+ 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=(20, 20), order=order)
+
+ dt_factor = 4
+ h = 1/20
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+
+ c = np.array([0.1,0.1])
+ norm_c = la.norm(c)
+
+
+ flux_type = "central"
+
+
+ def f(x):
+ return sym.sin(3*x)
+
+ def u_analytic(x):
+ return f(-np.dot(c, x)/norm_c+sym.var("t", sym.DD_SCALAR)*norm_c)
+
+ from grudge.models.advection import WeakAdvectionOperator
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+ op = WeakAdvectionOperator(c,
+ inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
+ flux_type=flux_type)
+
+ bound_op = bind(discr, op.sym_operator())
+
+ u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0)
+
+ def rhs(t, u):
+ return bound_op(queue, t=t, u=u)
+
+ final_time = 0.3
+
+ dt = dt_factor * h/order**2
+ nsteps = (final_time // dt) + 1
+ dt = final_time/nsteps + 1e-15
+
+
+ from grudge.shortcuts import set_up_rk4
+ dt_stepper = set_up_rk4("u", dt, u, rhs)
+
+ last_u = None
+
+ from grudge.shortcuts import make_visualizer
+ vis = make_visualizer(discr, vis_order=order)
+
+ step = 0
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ for event in dt_stepper.run(t_end=final_time):
+ if isinstance(event, dt_stepper.StateComputed):
+
+ step += 1
+
+ #print(step, event.t, norm(queue, u=event.state_component[0]))
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [ ("u", event.state_component) ])
+
+
+
+
+
+
+
+if __name__ == "__main__":
+ main()
+
+
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 8fda51a0..d50a06a3 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -133,214 +133,66 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
# {{{ variable-coefficient advection
class VariableCoefficientAdvectionOperator(HyperbolicOperator):
- """A class for space- and time-dependent DG advection operators.
-
- :param advec_v: Adheres to the :class:`grudge.data.ITimeDependentGivenFunction`
- interfacer and is an n-dimensional vector representing the velocity.
- :param bc_u_f: Adheres to the :class:`grudge.data.ITimeDependentGivenFunction`
- interface and is a scalar representing the boundary condition at all
- boundary faces.
-
- Optionally allows diffusion.
- """
-
- flux_types = [
- "central",
- "upwind",
- "lf"
- ]
-
- def __init__(self,
- dimensions,
- advec_v,
- bc_u_f="None",
- flux_type="central",
- diffusion_coeff=None,
- diffusion_scheme=CentralSecondDerivative()):
- self.dimensions = dimensions
- self.advec_v = advec_v
- self.bc_u_f = bc_u_f
+ def __init__(self, dim, v, inflow_u, flux_type="central"):
+ self.ambient_dim = dim
+ self.v = v
+ self.inflow_u = inflow_u
self.flux_type = flux_type
- self.diffusion_coeff = diffusion_coeff
- self.diffusion_scheme = diffusion_scheme
-
- # {{{ flux ----------------------------------------------------------------
- def flux(self):
- from grudge.flux import (
- make_normal,
- FluxVectorPlaceholder,
- flux_max)
- from pymbolic.primitives import IfPositive
-
- d = self.dimensions
-
- w = FluxVectorPlaceholder((1+d)+1)
- u = w[0]
- v = w[1:d+1]
- c = w[1+d]
-
- normal = make_normal(self.dimensions)
-
+ def flux(self, u):
+ normal = sym.normal(u. dd, self.ambient_dim)
+
+ #v = sym.nodes(2,u.dd)#self.v
+ surf_v = sym.interp("vol", u.dd)(self.v)
+
+
+ v_dot_normal = sym.cse(np.dot(surf_v,normal), "v_dot_normal")
+ norm_v = sym.sqrt(np.sum(self.v**2))
+
if self.flux_type == "central":
- return (u.int*np.dot(v.int, normal)
- + u.ext*np.dot(v.ext, normal)) * 0.5
- elif self.flux_type == "lf":
- n_vint = np.dot(normal, v.int)
- n_vext = np.dot(normal, v.ext)
- return 0.5 * (n_vint * u.int + n_vext * u.ext) \
- - 0.5 * (u.ext - u.int) \
- * flux_max(c.int, c.ext)
+ #return u.avg*v_dot_normal
+ # versus??
+ #return v_dot_normal
+ return (u.int*v_dot_normal
+ + u.ext*v_dot_normal) * 0.5
+ elif self.flux_type == "lf":
+ return u.avg*v_dot_normal + 0.5*norm_v*(u.int - u.ext)
elif self.flux_type == "upwind":
return (
- IfPositive(
- np.dot(normal, v.avg),
- np.dot(normal, v.int) * u.int, # outflow
- np.dot(normal, v.ext) * u.ext, # inflow
+ v_dot_normal * sym.If(
+ sym.Comparison(v_dot_normal, ">", 0),
+ u.int, # outflow
+ u.ext, # inflow
))
else:
raise ValueError("invalid flux type")
- # }}}
-
- def bind_characteristic_velocity(self, discr):
- from grudge.symbolic.operators import (
- ElementwiseMaxOperator)
- from grudge.symbolic import make_sym_vector
- velocity_vec = make_sym_vector("v", self.dimensions)
- velocity = ElementwiseMaxOperator()(
- np.dot(velocity_vec, velocity_vec)**0.5)
-
- compiled = discr.compile(velocity)
-
- def do(t, u):
- return compiled(v=self.advec_v.volume_interpolant(t, discr))
-
- return do
-
- def sym_operator(self, with_sensor=False):
- # {{{ operator preliminaries ------------------------------------------
- from grudge.symbolic import (Field, BoundaryPair, get_flux_operator,
- make_stiffness_t, InverseMassOperator, make_sym_vector,
- ElementwiseMaxOperator, RestrictToBoundary)
-
- from grudge.symbolic.primitives import make_common_subexpression as cse
-
- from grudge.symbolic.operators import (
- QuadratureGridUpsampler,
- QuadratureInteriorFacesGridUpsampler)
-
- to_quad = QuadratureGridUpsampler("quad")
- to_if_quad = QuadratureInteriorFacesGridUpsampler("quad")
-
- from grudge.tools import join_fields, \
- ptwise_dot
-
- u = Field("u")
- v = make_sym_vector("v", self.dimensions)
- c = ElementwiseMaxOperator()(ptwise_dot(1, 1, v, v))
- quad_u = cse(to_quad(u))
- quad_v = cse(to_quad(v))
-
- w = join_fields(u, v, c)
- quad_face_w = to_if_quad(w)
- # }}}
-
- # {{{ boundary conditions ---------------------------------------------
-
- from grudge.mesh import BTAG_ALL
- bc_c = to_quad(RestrictToBoundary(BTAG_ALL)(c))
- bc_u = to_quad(Field("bc_u"))
- bc_v = to_quad(RestrictToBoundary(BTAG_ALL)(v))
-
- if self.bc_u_f is "None":
- bc_w = join_fields(0, bc_v, bc_c)
- else:
- bc_w = join_fields(bc_u, bc_v, bc_c)
-
- minv_st = make_stiffness_t(self.dimensions)
- m_inv = InverseMassOperator()
-
- flux_op = get_flux_operator(self.flux())
- # }}}
-
- # {{{ diffusion -------------------------------------------------------
- if with_sensor or (
- self.diffusion_coeff is not None and self.diffusion_coeff != 0):
- if self.diffusion_coeff is None:
- diffusion_coeff = 0
- else:
- diffusion_coeff = self.diffusion_coeff
-
- if with_sensor:
- diffusion_coeff += Field("sensor")
-
- from grudge.second_order import SecondDerivativeTarget
-
- # strong_form here allows IPDG to reuse the value of grad u.
- grad_tgt = SecondDerivativeTarget(
- self.dimensions, strong_form=True,
- operand=u)
-
- self.diffusion_scheme.grad(grad_tgt, bc_getter=None,
- dirichlet_tags=[], neumann_tags=[])
-
- div_tgt = SecondDerivativeTarget(
- self.dimensions, strong_form=False,
- operand=diffusion_coeff*grad_tgt.minv_all)
-
- self.diffusion_scheme.div(div_tgt,
- bc_getter=None,
- dirichlet_tags=[], neumann_tags=[])
-
- diffusion_part = div_tgt.minv_all
- else:
- diffusion_part = 0
-
- # }}}
-
- to_quad = QuadratureGridUpsampler("quad")
- quad_u = cse(to_quad(u))
- quad_v = cse(to_quad(v))
-
- return m_inv(np.dot(minv_st, cse(quad_v*quad_u))
- - (flux_op(quad_face_w)
- + flux_op(BoundaryPair(quad_face_w, bc_w, BTAG_ALL)))) \
- + diffusion_part
-
- def bind(self, discr, sensor=None):
- compiled_sym_operator = discr.compile(
- self.sym_operator(with_sensor=sensor is not None))
+ def sym_operator(self):
+ u = sym.var("u")
- from grudge.mesh import check_bc_coverage, BTAG_ALL
- check_bc_coverage(discr.mesh, [BTAG_ALL])
+ # boundary conditions -------------------------------------------------
+ bc_in = self.inflow_u
+ # bc_out = sym.interp("vol", self.outflow_tag)(u)
- def rhs(t, u):
- kwargs = {}
- if sensor is not None:
- kwargs["sensor"] = sensor(t, u)
+ def flux(pair):
+ return sym.interp(pair.dd, "all_faces")(
+ self.flux(pair))
- if self.bc_u_f is not "None":
- kwargs["bc_u"] = \
- self.bc_u_f.boundary_interpolant(t, discr, tag=BTAG_ALL)
- return compiled_sym_operator(
- u=u,
- v=self.advec_v.volume_interpolant(t, discr),
- **kwargs)
+ return sym.InverseMassOperator()(
+ np.dot(
+ #self.v, sym.stiffness_t(self.ambient_dim)*u)
+ 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))
- return rhs
+ # FIXME: Add back support for inflow/outflow tags
+ #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
+ #+ flux(sym.bv_tpair(self.outflow_tag, u, bc_out))
+ ))
- def max_eigenvalue(self, t, fields=None, discr=None):
- # Gives the max eigenvalue of a vector of eigenvalues.
- # As the velocities of each node is stored in the velocity-vector-field
- # a pointwise dot product of this vector has to be taken to get the
- # magnitude of the velocity at each node. From this vector the maximum
- # values limits the timestep.
- from grudge.tools import ptwise_dot
- v = self.advec_v.volume_interpolant(t, discr)
- return discr.nodewise_max(ptwise_dot(1, 1, v, v)**0.5)
# }}}
--
GitLab