From 41fa65dfe555581844fb47a5d60395c2c42b4134 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 10:26:20 -0600 Subject: [PATCH 1/7] 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 . +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 . + + +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 From b730db90d10e77aa315409251914d9e4739ceaec Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 10:47:31 -0600 Subject: [PATCH 2/7] Cleaned up models/advection.py --- grudge/models/advection.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index d50a06a3..dda57c6a 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -141,7 +141,6 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): 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) @@ -149,11 +148,11 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): norm_v = sym.sqrt(np.sum(self.v**2)) if self.flux_type == "central": - #return u.avg*v_dot_normal + return u.avg*v_dot_normal # versus?? #return v_dot_normal - return (u.int*v_dot_normal - + u.ext*v_dot_normal) * 0.5 + #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) @@ -172,7 +171,6 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): # 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")( -- GitLab From 79456c15a864935a46f66d1cef8f21807d412403 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 10:54:28 -0600 Subject: [PATCH 3/7] I cleaned the code more, to make the Pipeline happy --- grudge/models/advection.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index dda57c6a..a500918f 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -169,7 +169,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): def sym_operator(self): u = sym.var("u") - # boundary conditions ------------------------------------------------- + # boundary conditions ------------------------------------------------- bc_in = self.inflow_u def flux(pair): @@ -179,20 +179,15 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): 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()( + - sym.FaceMassOperator()( flux(sym.int_tpair(u)) - + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in)) + + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in)) # 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)) )) - - - # }}} - # vim: foldmethod=marker -- GitLab From 5600883d68bc1f3a4caa1545ae2fdc1d4553318b Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 12:58:33 -0600 Subject: [PATCH 4/7] Pleasing the linter --- grudge/models/advection.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index a500918f..707a729d 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -30,7 +30,6 @@ import numpy as np import numpy.linalg as la from grudge.models import HyperbolicOperator -from grudge.models.second_order import CentralSecondDerivative from grudge import sym @@ -134,25 +133,21 @@ class WeakAdvectionOperator(AdvectionOperatorBase): class VariableCoefficientAdvectionOperator(HyperbolicOperator): def __init__(self, dim, v, inflow_u, flux_type="central"): - self.ambient_dim = dim + self.ambient_dim = dim self.v = v self.inflow_u = inflow_u self.flux_type = flux_type - def flux(self, u): + + 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") + + 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 - # 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) @@ -176,7 +171,6 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): return sym.interp(pair.dd, "all_faces")( self.flux(pair)) - return sym.InverseMassOperator()( np.dot( sym.stiffness_t(self.ambient_dim), sym.cse(self.v*u)) @@ -187,7 +181,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): # 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)) - )) + )) # }}} # vim: foldmethod=marker -- GitLab From 7248dc309cafe88efc92d26071caca39233635df Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 13:37:51 -0600 Subject: [PATCH 5/7] Linter --- grudge/models/advection.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 707a729d..68b1a0ba 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -172,16 +172,15 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): self.flux(pair)) return sym.InverseMassOperator()( - np.dot( - 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)) + np.dot(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)) - # 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)) - )) + # 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)) + )) # }}} # vim: foldmethod=marker -- GitLab From ac8c4816f311f059df78ad1eb8b472d252e65bac Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 13:53:40 -0600 Subject: [PATCH 6/7] Linter --- grudge/models/advection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 68b1a0ba..e10fb000 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -174,8 +174,8 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): return sym.InverseMassOperator()( np.dot(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)) + flux(sym.int_tpair(u)) + + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in)) # FIXME: Add back support for inflow/outflow tags #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in)) -- GitLab From 90acbf2ef7cdb13bb740e5f370d2d639a648ef03 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Feb 2017 14:22:15 -0600 Subject: [PATCH 7/7] Linting locally is way easier --- examples/advection/var-velocity.py | 41 +++--------------------------- grudge/models/advection.py | 6 ++--- 2 files changed, 7 insertions(+), 40 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index e022c470..8611b29d 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -16,7 +16,6 @@ import numpy as np -import math import pyopencl as cl # noqa import pyopencl.array # noqa import pyopencl.clmath # noqa @@ -30,22 +29,15 @@ import logging logger = logging.getLogger(__name__) from grudge import sym, bind, Discretization - -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) @@ -55,28 +47,17 @@ def main(write_output=True, order=4): 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 - + advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 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( @@ -86,12 +67,10 @@ def main(write_output=True, order=4): def u_analytic(x): return 0 - from grudge.models.advection import VariableCoefficientAdvectionOperator, WeakAdvectionOperator - from meshmode.mesh import BTAG_ALL, BTAG_NONE + from grudge.models.advection import VariableCoefficientAdvectionOperator - discr = Discretization(cl_ctx, mesh, order=order) - op = VariableCoefficientAdvectionOperator(2,advec_v, + op = VariableCoefficientAdvectionOperator(2, advec_v, inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)), flux_type=flux_type) @@ -108,19 +87,14 @@ def main(write_output=True, order=4): 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): @@ -128,15 +102,8 @@ def main(write_output=True, order=4): #print(step, event.t, norm(queue, u=event.state_component[0])) vis.write_vtk_file("fld-%04d.vtu" % step, - [ ("u", event.state_component) ]) - - - - - + [("u", event.state_component)]) if __name__ == "__main__": main() - - diff --git a/grudge/models/advection.py b/grudge/models/advection.py index e10fb000..1c83ed96 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -177,9 +177,9 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): flux(sym.int_tpair(u)) + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in)) - # 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)) + # 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)) )) # }}} -- GitLab