diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index bc720731a820fd063fda36657218315174f3c5f3..b6beb7a2d493aa77a154b96b8e20f6461fb757cd 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -22,114 +22,199 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-
+import os
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
+# {{{ plotting (keep in sync with `weak.py`)
-FINAL_TIME = 5
+class Plotter:
+ def __init__(self, queue, discr, order, visualize=True, ylim=None):
+ self.queue = queue
+ self.dim = discr.ambient_dim
+ self.visualize = visualize
+ if not self.visualize:
+ return
-def main(write_output=True, order=4):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
+ if self.dim == 1:
+ import matplotlib.pyplot as pt
+ self.fig = pt.figure(figsize=(8, 8), dpi=300)
+ self.ylim = ylim
- dim = 2
+ volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
+ self.x = volume_discr.nodes().get(self.queue)
+ else:
+ from grudge.shortcuts import make_visualizer
+ self.vis = make_visualizer(discr, vis_order=order)
- 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)
+ def __call__(self, evt, basename, overwrite=True):
+ if not self.visualize:
+ return
- dt_factor = 5
- h = 1/resolution
+ if self.dim == 1:
+ u = evt.state_component.get(self.queue)
- sym_x = sym.nodes(2)
+ filename = "%s.png" % basename
+ if not overwrite and os.path.exists(filename):
+ from meshmode import FileExistsError
+ raise FileExistsError("output file '%s' already exists" % filename)
- advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
+ ax = self.fig.gca()
+ ax.plot(self.x[0], u, '-')
+ ax.plot(self.x[0], u, 'k.')
+ if self.ylim is not None:
+ ax.set_ylim(self.ylim)
- flux_type = "upwind"
+ ax.set_xlabel("$x$")
+ ax.set_ylabel("$u$")
+ ax.set_title("t = {:.2f}".format(evt.t))
+ self.fig.savefig(filename)
+ self.fig.clf()
+ else:
+ self.vis.write_vtk_file("%s.vtu" % basename, [
+ ("u", evt.state_component)
+ ], overwrite=overwrite)
- source_center = np.array([0.1, 0.1])
- source_width = 0.05
+# }}}
- sym_x = sym.nodes(2)
- sym_source_center_dist = 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 main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
+ cl_ctx = ctx_factory()
+ queue = cl.CommandQueue(cl_ctx)
- 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)
+ # {{{ parameters
- def u_analytic(x):
- return 0
+ # domain [0, d]^dim
+ d = 1.0
+ # number of points in each dimension
+ npoints = 25
+ # grid spacing
+ h = d / npoints
- from grudge.models.advection import VariableCoefficientAdvectionOperator
- from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # noqa
+ # 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
+
+ # 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]
+
+ # }}}
+
+ # {{{ discretization
+ 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={
- #"product": None,
- "product": QuadratureSimplexGroupFactory(order=4*order)
- })
+ quad_tag_to_group_factory=quad_tag_to_group_factory)
+
+ # }}}
+
+ # {{{ symbolic operators
+
+ # gaussian parameters
+ source_center = np.array([0.5, 0.75, 0.0])[:dim]
+ source_width = 0.05
+ dist_squared = np.dot(sym_x - source_center, sym_x - source_center)
+
+ def f_gaussian(x):
+ return sym.exp(-dist_squared / source_width**2)
- op = VariableCoefficientAdvectionOperator(2, advec_v,
- u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product",
- flux_type=flux_type)
+ def f_step(x):
+ return sym.If(sym.Comparison(
+ dist_squared, "<", (4*source_width) ** 2),
+ 1, 0)
- bound_op = bind(
- discr, op.sym_operator(),
- #debug_flags=["dump_sym_operator_stages"]
- )
+ def u_bc(x):
+ return 0.0
+ from grudge.models.advection import VariableCoefficientAdvectionOperator
+ 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)
+ plot = Plotter(queue, discr, order, visualize=visualize,
+ ylim=[-0.1, 1.1])
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 % 10 == 0:
+ norm_u = norm(queue, u=event.state_component)
+ plot(event, "fld-var-velocity-%04d" % step)
- step += 1
- if step % 30 == 0:
- print(step)
+ step += 1
+ logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
- 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 6e30094d080033f225a9dbc45227e0d0c7665ad0..14ed08e1a757d149c1c96ecb9d84eb60d15400b1 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -1,18 +1,23 @@
-# 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 .
+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 .
+"""
+
+import os
import numpy as np
import numpy.linalg as la
@@ -23,9 +28,61 @@ import pyopencl.clmath
from grudge import bind, sym
import logging
-
logger = logging.getLogger(__name__)
-logging.basicConfig(level=logging.INFO)
+
+
+# {{{ plotting (keep in sync with `var-velocity.py`)
+
+class Plotter:
+ def __init__(self, queue, discr, order, visualize=True, ylim=None):
+ self.queue = queue
+ self.dim = discr.ambient_dim
+
+ 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)
+ self.ylim = ylim
+
+ volume_discr = discr.discr_from_dd(sym.DD_VOLUME)
+ self.x = volume_discr.nodes().get(self.queue)
+ else:
+ from grudge.shortcuts import make_visualizer
+ self.vis = make_visualizer(discr, vis_order=order)
+
+ def __call__(self, evt, basename, overwrite=True):
+ if not self.visualize:
+ return
+
+ if self.dim == 1:
+ u = evt.state_component.get(self.queue)
+
+ filename = "%s.png" % basename
+ if not overwrite and os.path.exists(filename):
+ from meshmode import FileExistsError
+ raise FileExistsError("output file '%s' already exists" % filename)
+
+ ax = self.fig.gca()
+ ax.plot(self.x[0], u, '-')
+ ax.plot(self.x[0], u, 'k.')
+ if self.ylim is not None:
+ ax.set_ylim(self.ylim)
+
+ ax.set_xlabel("$x$")
+ ax.set_ylabel("$u$")
+ ax.set_title("t = {:.2f}".format(evt.t))
+
+ self.fig.savefig(filename)
+ self.fig.clf()
+ else:
+ self.vis.write_vtk_file("%s.vtu" % basename, [
+ ("u", evt.state_component)
+ ], overwrite=overwrite)
+
+# }}}
def main(ctx_factory, dim=2, order=4, visualize=True):
@@ -34,16 +91,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 +113,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 +124,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,48 +147,31 @@ 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)
+ # {{{ time stepping
- def plot_solution(evt):
- if not visualize:
- return
+ from grudge.shortcuts import set_up_rk4
+ dt_stepper = set_up_rk4("u", dt, u, rhs)
+ plot = Plotter(queue, discr, order, visualize=visualize,
+ ylim=[-1.1, 1.1])
- 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)
+ norm = bind(discr, sym.norm(2, sym.var("u")))
step = 0
- norm = bind(discr, sym.norm(2, sym.var("u")))
+ norm_u = 0.0
for event in dt_stepper.run(t_end=final_time):
if not isinstance(event, dt_stepper.StateComputed):
continue
+ if step % 10 == 0:
+ norm_u = norm(queue, u=event.state_component)
+ plot(event, "fld-weak-%04d" % step)
+
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)
+
+ # }}}
if __name__ == "__main__":
@@ -143,5 +181,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 c3b142e411cd5fd7b34b17bce69e529ff565d316..36ce5c155bde62a7314a3ec6f0d828e8968946a6 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))
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 8c2f7fd0fecad5ef4c4dcae4e2d708d317edc3a7..1950d3f7248705af923d1e9148294d7e402a3028 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -307,8 +307,11 @@ class RefStiffnessTOperator(RefDiffOperatorBase):
grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes)
vand_inv_t = np.linalg.inv(vand).T
- weights = in_elem_grp.weights
+ if not isinstance(grad_vand, tuple):
+ # NOTE: special case for 1d
+ grad_vand = (grad_vand,)
+ weights = in_elem_grp.weights
return np.einsum('c,bz,acz->abc', weights, vand_inv_t, grad_vand)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 12b7b116e010006c957d34a210d3455e223c25c4..f068d34d8ea8b3de15053d20187ff13c0ed66223 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -460,7 +460,7 @@ def test_improvement_quadrature(ctx_factory, order):
advec_v = join_fields(-1*sym_nds[1], sym_nds[0])
flux = "upwind"
- op = VariableCoefficientAdvectionOperator(2, advec_v, 0, flux_type=flux)
+ op = VariableCoefficientAdvectionOperator(advec_v, 0, flux_type=flux)
def gaussian_mode():
source_width = 0.1