From 625db070a9bbb957bbc4e69965c82ff714a6110c Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Sat, 9 May 2020 20:39:49 -0500
Subject: [PATCH] examples: clean up variable velocity example
---
examples/advection/var-velocity.py | 165 +++++++++++++++++------------
examples/advection/weak.py | 137 ++++++++++++++----------
grudge/models/advection.py | 113 +++++++++-----------
3 files changed, 224 insertions(+), 191 deletions(-)
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index bc720731..c8791efd 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -22,114 +22,141 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-
import numpy as np
-import pyopencl as cl # noqa
-import pyopencl.array # noqa
-import pyopencl.clmath # noqa
-import pytest # noqa
+import pyopencl as cl
+import pyopencl.array
+import pyopencl.clmath
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+from grudge import bind, sym
+from pytools.obj_array import join_fields
import logging
logger = logging.getLogger(__name__)
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-from pytools.obj_array import join_fields
+def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
+ cl_ctx = ctx_factory()
+ queue = cl.CommandQueue(cl_ctx)
-FINAL_TIME = 5
+ # {{{ parameters
+ # domain [0, d]^dim
+ d = 1.0
+ # number of points in each dimension
+ npoints = 25
+ # grid spacing
+ h = d / npoints
-def main(write_output=True, order=4):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
+ # cfl
+ dt_factor = 1.0
+ # finale time
+ final_time = 0.5
+ # time steps
+ dt = dt_factor * h/order**2
+ nsteps = int(final_time // dt) + 1
+ dt = final_time/nsteps + 1.0e-15
- dim = 2
+ # flux
+ flux_type = "upwind"
+ # velocity field
+ sym_x = sym.nodes(dim)
+ if dim == 1:
+ c = sym_x
+ else:
+ # solid body rotation
+ c = join_fields(
+ np.pi * (d/2 - sym_x[1]),
+ np.pi * (sym_x[0] - d/2),
+ 0)[:dim]
- resolution = 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=(resolution, resolution), order=order)
+ # }}}
- dt_factor = 5
- h = 1/resolution
+ # {{{ discretization
- sym_x = sym.nodes(2)
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(0,)*dim, b=(d,)*dim,
+ n=(npoints,)*dim,
+ order=order)
+
+ from meshmode.discretization.poly_element import \
+ QuadratureSimplexGroupFactory
+
+ if product_tag:
+ quad_tag_to_group_factory = {
+ product_tag: QuadratureSimplexGroupFactory(order=4*order)
+ }
+ else:
+ quad_tag_to_group_factory = {}
+
+ from grudge import DGDiscretizationWithBoundaries
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+ quad_tag_to_group_factory=quad_tag_to_group_factory)
- advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
+ # }}}
- flux_type = "upwind"
+ # {{{ symbolic operators
- source_center = np.array([0.1, 0.1])
+ # gaussian parameters
+ source_center = np.array([0.5, 0.75, 0.0])[:dim]
source_width = 0.05
-
- sym_x = sym.nodes(2)
- sym_source_center_dist = sym_x - source_center
+ dist_squared = np.dot(sym_x - source_center, sym_x - source_center)
def f_gaussian(x):
- return sym.exp(
- -np.dot(sym_source_center_dist, sym_source_center_dist)
- / source_width**2)
-
- def f_step(x):
- return sym.If(
- sym.Comparison(
- np.dot(sym_source_center_dist, sym_source_center_dist),
- "<",
- (4*source_width)**2),
- 1, 0)
+ return sym.exp(-dist_squared / source_width**2)
- def u_analytic(x):
- return 0
+ def u_bc(x):
+ return 0.0
from grudge.models.advection import VariableCoefficientAdvectionOperator
- from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # noqa
-
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
- quad_tag_to_group_factory={
- #"product": None,
- "product": QuadratureSimplexGroupFactory(order=4*order)
- })
-
- op = VariableCoefficientAdvectionOperator(2, advec_v,
- u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product",
- flux_type=flux_type)
-
- bound_op = bind(
- discr, op.sym_operator(),
- #debug_flags=["dump_sym_operator_stages"]
- )
+ op = VariableCoefficientAdvectionOperator(
+ c,
+ u_bc(sym.nodes(dim, sym.BTAG_ALL)),
+ quad_tag=product_tag,
+ flux_type=flux_type)
+ bound_op = bind(discr, op.sym_operator())
u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0)
def rhs(t, u):
return bound_op(queue, t=t, u=u)
- dt = dt_factor * h/order**2
- nsteps = (FINAL_TIME // dt) + 1
- dt = FINAL_TIME/nsteps + 1e-15
+ # }}}
+
+ # {{{ time stepping
from grudge.shortcuts import set_up_rk4
dt_stepper = set_up_rk4("u", dt, u, rhs)
- from grudge.shortcuts import make_visualizer
- vis = make_visualizer(discr, vis_order=2*order)
+ from weak import Plotter
+ plot = Plotter(queue, discr, order, visualize=visualize)
step = 0
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+ for event in dt_stepper.run(t_end=final_time):
+ if not isinstance(event, dt_stepper.StateComputed):
+ continue
- for event in dt_stepper.run(t_end=FINAL_TIME):
- if isinstance(event, dt_stepper.StateComputed):
+ if step % 5 == 0:
+ norm_u = norm(queue, u=event.state_component)
- step += 1
- if step % 30 == 0:
- print(step)
+ step += 1
+ logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
+ plot(event, "fld-var-velocity-%04d" % step)
- vis.write_vtk_file("fld-var-velocity-%04d.vtu" % step,
- [("u", event.state_component)])
+ # }}}
if __name__ == "__main__":
- main()
+ import argparse
+
+ parser = argparse.ArgumentParser()
+ parser.add_argument("--dim", default=2, type=int)
+ parser.add_argument("--qtag", default="product")
+ args = parser.parse_args()
+
+ logging.basicConfig(level=logging.INFO)
+ main(cl.create_some_context,
+ dim=args.dim,
+ product_tag=args.qtag)
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 6e30094d..770e9717 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -1,17 +1,21 @@
-# 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/>.
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
+
+__license__ = """
+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 numpy.linalg as la
@@ -23,9 +27,49 @@ import pyopencl.clmath
from grudge import bind, sym
import logging
-
logger = logging.getLogger(__name__)
-logging.basicConfig(level=logging.INFO)
+
+
+class Plotter:
+ def __init__(self, queue, discr, order, visualize=True):
+ volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
+
+ self.queue = queue
+ self.x = volume_discr.nodes().get(self.queue)
+
+ self.visualize = visualize
+ if not self.visualize:
+ return
+
+ if self.dim == 1:
+ import matplotlib.pyplot as pt
+ self.fig = pt.figure(figsize=(8, 8), dpi=300)
+ else:
+ from grudge.shortcuts import make_visualizer
+ self.vis = make_visualizer(discr, vis_order=order)
+
+ @property
+ def dim(self):
+ return len(self.x)
+
+ def __call__(self, evt, basename):
+ if not self.visualize:
+ return
+
+ if self.dim == 1:
+ u = evt.state_component.get(self.queue)
+
+ ax = self.fig.gca()
+ ax.plot(self.x[0], u, 'o')
+ ax.set_xlabel("$x$")
+ ax.set_ylabel("$u$")
+ ax.set_title("t = {:.2f}".format(evt.t))
+ self.fig.savefig("%s.png" % basename)
+ self.fig.clf()
+ else:
+ self.vis.write_vtk_file("%s.vtu" % basename, [
+ ("u", evt.state_component)
+ ], overwrite=True)
def main(ctx_factory, dim=2, order=4, visualize=True):
@@ -34,16 +78,21 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
# {{{ parameters
- # domain side [-d/2, d/2]^dim
+ # domain [-d/2, d/2]^dim
d = 1.0
# number of points in each dimension
npoints = 20
# grid spacing
h = d / npoints
- # cfl?
+
+ # cfl
dt_factor = 2.0
# final time
final_time = 1.0
+ # compute number of steps
+ dt = dt_factor * h/order**2
+ nsteps = int(final_time // dt) + 1
+ dt = final_time/nsteps + 1.0e-15
# velocity field
c = np.array([0.5] * dim)
@@ -51,11 +100,6 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
# flux
flux_type = "central"
- # compute number of steps
- dt = dt_factor * h / order**2
- nsteps = int(final_time // dt) + 1
- dt = final_time/nsteps + 1.0e-15
-
# }}}
# {{{ discretization
@@ -67,18 +111,16 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
from grudge import DGDiscretizationWithBoundaries
discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
- volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
# }}}
- # {{{ solve advection
+ # {{{ symbolic operators
def f(x):
return sym.sin(3 * x)
- def u_analytic(x, t=None):
- if t is None:
- t = sym.var("t", sym.DD_SCALAR)
+ def u_analytic(x):
+ t = sym.var("t", sym.DD_SCALAR)
return f(-np.dot(c, x) / norm_c + t * norm_c)
from grudge.models.advection import WeakAdvectionOperator
@@ -92,37 +134,13 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
def rhs(t, u):
return bound_op(queue, t=t, u=u)
- from grudge.shortcuts import set_up_rk4
- dt_stepper = set_up_rk4("u", dt, u, rhs)
-
- if dim == 1:
- import matplotlib.pyplot as pt
- pt.figure(figsize=(8, 8), dpi=300)
-
- volume_x = volume_discr.nodes().get(queue)
- else:
- from grudge.shortcuts import make_visualizer
- vis = make_visualizer(discr, vis_order=order)
+ # }}}
- def plot_solution(evt):
- if not visualize:
- return
+ # {{{ time stepping
- if dim == 1:
- u = event.state_component.get(queue)
- u_ = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=evt.t).get(queue)
-
- pt.plot(volume_x[0], u, 'o')
- pt.plot(volume_x[0], u_, "k.")
- pt.xlabel("$x$")
- pt.ylabel("$u$")
- pt.title("t = {:.2f}".format(evt.t))
- pt.savefig("fld-weak-%04d.png" % step)
- pt.clf()
- else:
- vis.write_vtk_file("fld-weak-%04d.vtu" % step, [
- ("u", evt.state_component)
- ], overwrite=True)
+ from grudge.shortcuts import set_up_rk4
+ dt_stepper = set_up_rk4("u", dt, u, rhs)
+ plot = Plotter(queue, discr, order, visualize=visualize)
step = 0
norm = bind(discr, sym.norm(2, sym.var("u")))
@@ -133,7 +151,9 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
step += 1
norm_u = norm(queue, u=event.state_component)
logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
- plot_solution(event)
+ plot(event, "fld-weak-%04d" % step)
+
+ # }}}
if __name__ == "__main__":
@@ -143,5 +163,6 @@ if __name__ == "__main__":
parser.add_argument("--dim", default=2, type=int)
args = parser.parse_args()
+ logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
dim=args.dim)
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index c3b142e4..36ce5c15 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -33,6 +33,31 @@ from grudge.models import HyperbolicOperator
from grudge import sym
+# {{{ fluxes
+
+def advection_weak_flux(flux_type, u, velocity):
+ normal = sym.normal(u.dd, len(velocity))
+ v_dot_n = sym.cse(velocity.dot(normal), "v_dot_normal")
+
+ flux_type = flux_type.lower()
+ if flux_type == "central":
+ return u.avg * v_dot_n
+ elif flux_type == "lf":
+ norm_v = sym.sqrt((velocity**2).sum())
+ return u.avg * v_dot_n + 0.5 * norm_v * (u.int - u.ext)
+ elif flux_type == "upwind":
+ u_upwind = sym.If(
+ sym.Comparison(v_dot_n, ">", 0),
+ u.int, # outflow
+ u.ext # inflow
+ )
+ return u_upwind * v_dot_n
+ else:
+ raise ValueError("flux `{}` is not implemented".format(flux_type))
+
+# }}}
+
+
# {{{ constant-coefficient advection
class AdvectionOperatorBase(HyperbolicOperator):
@@ -48,25 +73,11 @@ class AdvectionOperatorBase(HyperbolicOperator):
self.inflow_u = inflow_u
self.flux_type = flux_type
- def weak_flux(self, u):
- normal = sym.normal(u. dd, self.ambient_dim)
+ if flux_type not in self.flux_types:
+ raise ValueError("unknown flux type: {}".format(flux_type))
- v_dot_normal = sym.cse(self.v.dot(normal), "v_dot_normal")
- norm_v = sym.sqrt((self.v**2).sum())
-
- if self.flux_type == "central":
- return u.avg*v_dot_normal
- 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 (
- v_dot_normal * sym.If(
- sym.Comparison(v_dot_normal, ">", 0),
- u.int, # outflow
- u.ext, # inflow
- ))
- else:
- raise ValueError("invalid flux type")
+ def weak_flux(self, u):
+ return advection_weak_flux(self.flux_type, u, self.v)
def max_eigenvalue(self, t=None, fields=None, discr=None):
return la.norm(self.v)
@@ -106,14 +117,13 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
def sym_operator(self):
u = sym.var("u")
- # boundary conditions -------------------------------------------------
- bc_in = self.inflow_u
- # bc_out = sym.interp("vol", self.outflow_tag)(u)
-
def flux(pair):
return sym.interp(pair.dd, "all_faces")(
self.flux(pair))
+ bc_in = self.inflow_u
+ # bc_out = sym.interp("vol", self.outflow_tag)(u)
+
return sym.InverseMassOperator()(
np.dot(
self.v, sym.stiffness_t(self.ambient_dim)*u)
@@ -131,65 +141,40 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
# {{{ variable-coefficient advection
-class VariableCoefficientAdvectionOperator(HyperbolicOperator):
- 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
+class VariableCoefficientAdvectionOperator(AdvectionOperatorBase):
+ def __init__(self, v, inflow_u, flux_type="central", quad_tag="product"):
+ super(VariableCoefficientAdvectionOperator, self).__init__(
+ v, inflow_u, flux_type=flux_type)
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)
-
- 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.avg*v_dot_normal
-
- 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 (
- v_dot_normal * sym.If(
- sym.Comparison(v_dot_normal, ">", 0),
- u.int, # outflow
- u.ext, # inflow
- ))
- else:
- raise ValueError("invalid flux type")
+ surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v)
+ return advection_weak_flux(self.flux_type, u, surf_v)
def sym_operator(self):
u = sym.var("u")
- # boundary conditions -------------------------------------------------
- bc_in = self.inflow_u
-
- all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
- boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
-
def flux(pair):
- return sym.interp(pair.dd, all_faces_dd)(
- self.flux(pair))
+ return sym.interp(pair.dd, face_dd)(self.flux(pair))
- quad_dd = sym.DOFDesc("vol", self.quad_tag)
+ face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
+ boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+ quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag)
- to_quad = sym.interp("vol", quad_dd)
+ 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)
- stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol")
quad_v = to_quad(self.v)
quad_u = to_quad(u)
return sym.InverseMassOperator()(
- (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1]))
- - sym.FaceMassOperator(all_faces_dd, "vol")(
+ sum(stiff_t_op[n](quad_u * quad_v[n])
+ for n in range(self.ambient_dim))
+ - sym.FaceMassOperator(face_dd, sym.DD_VOLUME)(
flux(sym.int_tpair(u, self.quad_tag))
- + flux(sym.bv_tpair(boundary_dd, u, bc_in))
+ + flux(sym.bv_tpair(boundary_dd, u, self.inflow_u))
# FIXME: Add back support for inflow/outflow tags
#+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
--
GitLab