From 3e06ba6b0bdc0a2f73f6256e4a0e8cc1a7129a1f Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Mon, 20 Apr 2020 19:25:15 -0500
Subject: [PATCH] clean up advection example and doing it in 1d
---
examples/advection/weak.py | 127 +++++++++++++++++++++---------------
grudge/models/advection.py | 4 +-
grudge/symbolic/compiler.py | 5 ++
3 files changed, 81 insertions(+), 55 deletions(-)
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index bd9d61ca..4d3736b1 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -13,106 +13,127 @@
# 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 numpy.linalg as la
-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
import logging
+
logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.INFO)
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-import numpy.linalg as la
+def main(ctx_factory, dim=1, order=4, visualize=True):
+ cl_ctx = ctx_factory()
+ queue = cl.CommandQueue(cl_ctx)
+ # {{{ parameters
+ # domain side [-d/2, d/2]^dim
+ d = 1.0
+ # number of points in each dimension
+ npoints = 20
+ # grid spacing
+ h = d / npoints
+ # cfl?
+ dt_factor = 4
+ # final time
+ final_time = 1.0
-def main(write_output=True, order=4):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
+ # velocity field
+ c = np.array([0.1] * dim)
+ norm_c = la.norm(c)
+ # flux
+ flux_type = "central"
- dim = 2
+ # compute number of steps
+ dt = dt_factor * h / order**2
+ nsteps = int(final_time // dt) + 1
+ dt = final_time/nsteps + 1.0e-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=(20, 20), order=order)
+ # {{{ discretization
- dt_factor = 4
- h = 1/20
+ from meshmode.mesh.generation import generate_box_mesh
+ mesh = generate_box_mesh(
+ [np.linspace(-d/2, d/2, npoints) for _ in range(dim)],
+ order=order)
+ from grudge import DGDiscretizationWithBoundaries
discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
- c = np.array([0.1,0.1])
- norm_c = la.norm(c)
-
+ # }}}
- flux_type = "central"
-
+ # {{{ solve advection
def f(x):
- return sym.sin(3*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)
+ t = sym.var("t", sym.DD_SCALAR)
+ return f(-np.dot(c, x) / norm_c + t * norm_c)
from grudge.models.advection import WeakAdvectionOperator
- from meshmode.mesh import BTAG_ALL, BTAG_NONE
-
- discr = DGDiscretizationWithBoundaries(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)
+ return
+
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)
+ if dim == 1:
+ import matplotlib.pyplot as pt
+ pt.figure(figsize=(8, 8), dpi=300)
+ else:
+ 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-weak-%04d.vtu" % step,
- [ ("u", event.state_component) ])
-
-
+ if not isinstance(event, dt_stepper.StateComputed):
+ continue
+ step += 1
+ norm_u = norm(queue, u=event.state_component)
+ logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
+ if dim == 1:
+ x = discr.discr_from_dd(sym.DTAG_VOLUME_ALL).nodes().get(queue)
+ u = event.state_component.get(queue)
+ pt.plot(x, u)
+ pt.xlabel("$x$")
+ pt.ylabel("$u$")
+ pt.title("t = {:.2f}".format(event.t))
+ pt.savefig("fld-weak-%04d.png" % step)
+ pt.clf()
+ else:
+ vis.write_vtk_file("fld-weak-%04d.vtu" % step, [
+ ("u", event.state_component)
+ ], overwrite=True)
if __name__ == "__main__":
- main()
+ import argparse
+ parser = argparse.ArgumentParser()
+ parser.add_argument("--dim", default=2, type=int)
+ args = parser.parse_args()
+ main(cl.create_some_context,
+ dim=args.dim)
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 091b9870..c3b142e4 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -74,7 +74,7 @@ class AdvectionOperatorBase(HyperbolicOperator):
class StrongAdvectionOperator(AdvectionOperatorBase):
def flux(self, u):
- normal = sym.normal(u. dd, self.ambient_dim)
+ normal = sym.normal(u.dd, self.ambient_dim)
v_dot_normal = sym.cse(self.v.dot(normal), "v_dot_normal")
return u.int * v_dot_normal - self.weak_flux(u)
@@ -142,7 +142,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
def flux(self, u):
- normal = sym.normal(u. dd, self.ambient_dim)
+ normal = sym.normal(u.dd, self.ambient_dim)
surf_v = sym.interp("vol", u.dd)(self.v)
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 51eb5202..4e13caaf 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -1047,6 +1047,11 @@ class ToLoopyInstructionMapper(object):
self.insn_count += 1
+ for expr in insn.exprs:
+ print(expr)
+ print(self.dd_inference_mapper(expr))
+ print()
+
from pytools import single_valued
governing_dd = single_valued(
self.dd_inference_mapper(expr)
--
GitLab