From 896bd919595c5b6dcde3f609722f67c99db8fa31 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 2 Jun 2019 00:55:56 -0500 Subject: [PATCH 1/2] Towards reviving the Euler operator --- grudge/models/gas_dynamics/__init__.py | 694 +++++++++--------------- grudge/models/gas_dynamics/solutions.py | 88 +++ grudge/symbolic/primitives.py | 2 + grudge/tools.py | 41 ++ setup.cfg | 1 - test/test_gas_dynamics.py | 110 ++++ 6 files changed, 511 insertions(+), 425 deletions(-) create mode 100644 grudge/models/gas_dynamics/solutions.py create mode 100644 test/test_gas_dynamics.py diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py index 4df408d8..84cb76e3 100644 --- a/grudge/models/gas_dynamics/__init__.py +++ b/grudge/models/gas_dynamics/__init__.py @@ -1,11 +1,11 @@ """Operator for compressible Navier-Stokes and Euler equations.""" -from __future__ import division -from __future__ import absolute_import -import six -from six.moves import range +from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2007 Hendrik Riedmann, Andreas Kloeckner" +__copyright__ = """ +Copyright (C) 2007 Hendrik Riedmann, Andreas Kloeckner +Copyright (C) 2019 Andreas Kloeckner +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -29,42 +29,26 @@ THE SOFTWARE. +import six # noqa +from six.moves import range -import numpy -import grudge.tools -import grudge.mesh -import grudge.data -from grudge.models import TimeDependentOperator +import numpy as np from pytools import Record -from grudge.tools import is_zero -from grudge.second_order import ( - StabilizedCentralSecondDerivative, - CentralSecondDerivative, - IPDGSecondDerivative) -from grudge.symbolic.primitives import make_common_subexpression as cse +from grudge import sym + from pytools import memoize_method -from grudge.symbolic.tools import make_sym_vector from pytools.obj_array import make_obj_array, join_fields AXES = ["x", "y", "z", "w"] -from grudge.symbolic.operators import ( - QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler) - -to_vol_quad = QuadratureGridUpsampler("gasdyn_vol") - -# It is recommended (though not required) that these two -# remain the same so that they can be computed together -# by the CUDA backend - -to_int_face_quad = QuadratureInteriorFacesGridUpsampler("gasdyn_face") -to_bdry_quad = QuadratureGridUpsampler("gasdyn_face") - - +# FIXMEs: +# - On which dd should boundary conditions live? +# - Refactor so that flux computation happens exactly in LF flux, +# and not all over the map. # {{{ equations of state + class EquationOfState(object): def q_to_p(self, op, q): raise NotImplementedError @@ -72,6 +56,7 @@ class EquationOfState(object): def p_to_e(self, p, rho, u): raise NotImplementedError + class GammaLawEOS(EquationOfState): # FIXME Shouldn't gamma only occur in the equation of state? # I.e. shouldn't all uses of gamma go through the EOS? @@ -80,25 +65,23 @@ class GammaLawEOS(EquationOfState): self.gamma = gamma def q_to_p(self, op, q): - return (self.gamma-1)*(op.e(q)-0.5*numpy.dot(op.rho_u(q), op.u(q))) + return (self.gamma-1)*(op.e(q)-0.5*np.dot(op.rho_u(q), op.u(q))) def p_to_e(self, p, rho, u): - return p / (self.gamma - 1) + rho / 2 * numpy.dot(u, u) + return p / (self.gamma - 1) + rho / 2 * np.dot(u, u) + class PolytropeEOS(GammaLawEOS): # inverse is same as superclass def q_to_p(self, op, q): - return op.rho(q)**self.gamma + return op.rho(q)**self.gamma # }}} - - - -class GasDynamicsOperator(TimeDependentOperator): - """An nD Navier-Stokes and Euler operator. +class GasDynamicsOperator(object): + r"""An nD Navier-Stokes and Euler operator. see JSH, TW: Nodal Discontinuous Galerkin Methods p.320 and p.206 @@ -125,9 +108,13 @@ class GasDynamicsOperator(TimeDependentOperator): Field order is [rho E rho_u_x rho_u_y ...]. """ - # {{{ initialization ------------------------------------------------------ - def __init__(self, dimensions, - gamma=None, mu=0, + # FIXME: Convert docstring to reST/math + + # {{{ initialization + + def __init__( + self, dimensions, + mu=0, bc_inflow=None, bc_outflow=None, bc_noslip=None, @@ -140,43 +127,10 @@ class GasDynamicsOperator(TimeDependentOperator): wall_tag="wall", supersonic_inflow_tag="supersonic_inflow", supersonic_outflow_tag="supersonic_outflow", - source=None, - second_order_scheme=CentralSecondDerivative(), - artificial_viscosity_mode=None, + #second_order_scheme=CentralSecondDerivative(), + #artificial_viscosity_mode=None, + quad_tag="gas_dyn", ): - """ - :param source: should implement - :class:`grudge.data.IFieldDependentGivenFunction` - or be None. - - :param artificial_viscosity_mode: - """ - from grudge.data import ( - TimeConstantGivenFunction, - ConstantGivenFunction) - - if gamma is not None: - if equation_of_state is not None: - raise ValueError("can only specify one of gamma and equation_of_state") - - from warnings import warn - warn("argument gamma is deprecated in favor of equation_of_state", - DeprecationWarning, stacklevel=2) - - equation_of_state = GammaLawEOS(gamma) - - dull_bc = TimeConstantGivenFunction( - ConstantGivenFunction(make_obj_array( - [1, 1] + [0]*dimensions))) - if bc_inflow is None: - bc_inflow = dull_bc - if bc_outflow is None: - bc_outflow = dull_bc - if bc_noslip is None: - bc_noslip = dull_bc - if bc_supersonic_inflow is None: - bc_supersonic_inflow = dull_bc - self.dimensions = dimensions self.prandtl = prandtl @@ -195,37 +149,35 @@ class GasDynamicsOperator(TimeDependentOperator): self.supersonic_inflow_tag = supersonic_inflow_tag self.supersonic_outflow_tag = supersonic_outflow_tag - self.source = source self.equation_of_state = equation_of_state - self.second_order_scheme = second_order_scheme - - if artificial_viscosity_mode not in [ - "cns", "diffusion", "blended", None]: - raise ValueError("artificial_viscosity_mode has an invalid value") + # self.second_order_scheme = second_order_scheme - self.artificial_viscosity_mode = artificial_viscosity_mode + # if artificial_viscosity_mode not in [ + # "cns", "diffusion", "blended", None]: + # raise ValueError("artificial_viscosity_mode has an invalid value") + # self.artificial_viscosity_mode = artificial_viscosity_mode + self.quad_tag = quad_tag # }}} - # {{{ conversions --------------------------------------------------------- + def to_vol_quad(self, operand): + return sym.interp("vol", sym.DOFDesc("vol", self.quad_tag))(operand) + + # {{{ conversions + def state(self): - return make_sym_vector("q", self.dimensions+2) + return sym.make_sym_array("q", self.dimensions+2) @memoize_method def volq_state(self): - return cse(to_vol_quad(self.state()), "vol_quad_state") - - @memoize_method - def faceq_state(self): - return cse(to_int_face_quad(self.state()), "face_quad_state") + return sym.cse(self.to_vol_quad(self.state()), "vol_quad_state") - @memoize_method - def sensor(self): - from grudge.symbolic.primitives import Field - sensor = Field("sensor") + # @memoize_method + # def faceq_state(self): + # return sym.cse(to_int_face_quad(self.state()), "face_quad_state") def rho(self, q): return q[0] @@ -241,66 +193,68 @@ class GasDynamicsOperator(TimeDependentOperator): rho_u_i/self.rho(q) for rho_u_i in self.rho_u(q)]) - def p(self,q): + def p(self, q): return self.equation_of_state.q_to_p(self, q) def cse_u(self, q): - return cse(self.u(q), "u") + return sym.cse(self.u(q), "u") def cse_rho(self, q): - return cse(self.rho(q), "rho") + return sym.cse(self.rho(q), "rho") def cse_rho_u(self, q): - return cse(self.rho_u(q), "rho_u") + return sym.cse(self.rho_u(q), "rho_u") def cse_p(self, q): - return cse(self.p(q), "p") + return sym.cse(self.p(q), "p") def temperature(self, q): c_v = 1 / (self.equation_of_state.gamma - 1) * self.spec_gas_const - return (self.e(q)/self.rho(q) - 0.5 * numpy.dot(self.u(q), self.u(q))) / c_v + return (self.e(q)/self.rho(q) - 0.5 * np.dot(self.u(q), self.u(q))) / c_v def cse_temperature(self, q): - return cse(self.temperature(q), "temperature") - - def get_mu(self, q, to_quad_op): - """ - :param to_quad_op: If not *None*, represents an operator which transforms - nodal values onto a quadrature grid on which the returned :math:`\mu` - needs to be represented. In that case, *q* is assumed to already be on the - same quadrature grid. - """ - - if to_quad_op is None: - def to_quad_op(x): - return x - - if self.mu == "sutherland": - # Sutherland's law: !!!not tested!!! - t_s = 110.4 - mu_inf = 1.735e-5 - result = cse( - mu_inf * self.cse_temperature(q) ** 1.5 * (1 + t_s) - / (self.cse_temperature(q) + t_s), - "sutherland_mu") - else: - result = self.mu - - if self.artificial_viscosity_mode == "cns": - mapped_sensor = self.sensor() - else: - mapped_sensor = None - - if mapped_sensor is not None: - result = result + cse(to_quad_op(mapped_sensor), "quad_sensor") - - return cse(result, "mu") + return sym.cse(self.temperature(q), "temperature") + + # def get_mu(self, q, to_quad_op): + # r""" + # :param to_quad_op: If not *None*, represents an operator which transforms + # nodal values onto a quadrature grid on which the returned :math:`\mu` + # needs to be represented. In that case, *q* is assumed to already be on + # the same quadrature grid. + # """ + + # if to_quad_op is None: + # def to_quad_op(x): + # return x + + # if self.mu == "sutherland": + # from warnings import warn + # warn("Using Sutherland's law in gas dynamics: untested") + # t_s = 110.4 + # mu_inf = 1.735e-5 + # result = sym.cse( + # mu_inf * self.cse_temperature(q) ** 1.5 * (1 + t_s) + # / (self.cse_temperature(q) + t_s), + # "sutherland_mu") + # else: + # result = self.mu + + # if self.artificial_viscosity_mode == "cns": + # mapped_sensor = self.sensor() + # else: + # mapped_sensor = None + + # if mapped_sensor is not None: + # result = result + sym.cse(to_quad_op(mapped_sensor), "quad_sensor") + + # return sym.cse(result, "mu") def primitive_to_conservative(self, prims, use_cses=True): if not use_cses: from grudge.symbolic.primitives import make_common_subexpression as cse else: - def cse(x, name): return x + def cse(x, name): + return x rho = prims[0] p = prims[1] @@ -312,32 +266,22 @@ class GasDynamicsOperator(TimeDependentOperator): cse(e, "e"), cse(rho * u, "rho_u")) - def conservative_to_primitive(self, q, use_cses=True): - if use_cses: - from grudge.symbolic.primitives import make_common_subexpression as cse - else: - def cse(x, name): return x - + def conservative_to_primitive(self, q): return join_fields( self.rho(q), self.p(q), self.u(q)) - def characteristic_velocity_optemplate(self, state): - from grudge.symbolic.operators import ElementwiseMaxOperator - - from grudge.symbolic.primitives import FunctionSymbol - sqrt = FunctionSymbol("sqrt") - - sound_speed = cse(sqrt( + def sym_characteristic_velocity(self, state): + sound_speed = sym.cse(sym.sqrt( self.equation_of_state.gamma*self.cse_p(state)/self.cse_rho(state)), "sound_speed") u = self.cse_u(state) - speed = cse(sqrt(numpy.dot(u, u)), "norm_u") + sound_speed - return ElementwiseMaxOperator()(speed) + speed = sym.cse(sym.sqrt(np.dot(u, u)), "norm_u") + sound_speed + return sym.ElementwiseMaxOperator(dd)(speed) def bind_characteristic_velocity(self, discr): - state = make_sym_vector("q", self.dimensions+2) + state = sym.make_sym_array("q", self.dimensions+2) compiled = discr.compile( self.characteristic_velocity_optemplate(state)) @@ -352,85 +296,86 @@ class GasDynamicsOperator(TimeDependentOperator): # {{{ helpers for second-order part --------------------------------------- # {{{ compute gradient of state --------------------------------------- - def grad_of(self, var, faceq_var): - from grudge.second_order import SecondDerivativeTarget - grad_tgt = SecondDerivativeTarget( - self.dimensions, strong_form=False, - operand=var, - int_flux_operand=faceq_var, - bdry_flux_int_operand=faceq_var) + # def grad_of(self, var, faceq_var): + # from grudge.second_order import SecondDerivativeTarget + # grad_tgt = SecondDerivativeTarget( + # self.dimensions, strong_form=False, + # operand=var, + # int_flux_operand=faceq_var, + # bdry_flux_int_operand=faceq_var) - self.second_order_scheme.grad(grad_tgt, - bc_getter=self.get_boundary_condition_for, - dirichlet_tags=self.get_boundary_tags(), - neumann_tags=[]) + # self.second_order_scheme.grad(grad_tgt, + # bc_getter=self.get_boundary_condition_for, + # dirichlet_tags=self.get_boundary_tags(), + # neumann_tags=[]) - return grad_tgt.minv_all + # return grad_tgt.minv_all - def grad_of_state(self): - dimensions = self.dimensions + # def grad_of_state(self): + # dimensions = self.dimensions - state = self.state() + # state = self.state() - dq = numpy.zeros((len(state), dimensions), dtype=object) + # dq = np.zeros((len(state), dimensions), dtype=object) - for i in range(len(state)): - dq[i,:] = self.grad_of( - state[i], self.faceq_state()[i]) + # for i in range(len(state)): + # dq[i, :] = self.grad_of( + # state[i], self.faceq_state()[i]) - return dq + # return dq - def grad_of_state_func(self, func, of_what_descr): - return cse(self.grad_of( - func(self.volq_state()), - func(self.faceq_state())), - "grad_"+of_what_descr) + # def grad_of_state_func(self, func, of_what_descr): + # return sym.cse(self.grad_of( + # func(self.volq_state()), + # func(self.faceq_state())), + # "grad_"+of_what_descr) # }}} # {{{ viscous stress tensor - def tau(self, to_quad_op, state, mu=None): - faceq_state = self.faceq_state() + # def tau(self, to_quad_op, state, mu=None): + # dimensions = self.dimensions - dimensions = self.dimensions + # # {{{ compute gradient of u + # # Use the product rule to compute the gradient of + # # u from the gradient of (rho u). This ensures we don't + # # compute the derivatives twice. - # {{{ compute gradient of u --------------------------------------- - # Use the product rule to compute the gradient of - # u from the gradient of (rho u). This ensures we don't - # compute the derivatives twice. + # from pytools.obj_array import with_object_array_or_scalar + # dq = with_object_array_or_scalar( + # to_quad_op, self.grad_of_state()) - from pytools.obj_array import with_object_array_or_scalar - dq = with_object_array_or_scalar( - to_quad_op, self.grad_of_state()) + # q = sym.cse(to_quad_op(state)) - q = cse(to_quad_op(state)) + # du = np.zeros((dimensions, dimensions), dtype=object) + # for i in range(dimensions): + # for j in range(dimensions): + # du[i, j] = sym.cse( + # (dq[i+2, j] - self.cse_u(q)[i] * dq[0, j]) / self.rho(q), + # "du%d_d%s" % (i, AXES[j])) - du = numpy.zeros((dimensions, dimensions), dtype=object) - for i in range(dimensions): - for j in range(dimensions): - du[i,j] = cse( - (dq[i+2,j] - self.cse_u(q)[i] * dq[0,j]) / self.rho(q), - "du%d_d%s" % (i, AXES[j])) + # # }}} - # }}} + # # {{{ put together viscous stress tau - # {{{ put together viscous stress tau ----------------------------- - from pytools import delta + # from pytools import delta - if mu is None: - mu = self.get_mu(q, to_quad_op) + # if mu is None: + # mu = self.get_mu(q, to_quad_op) - tau = numpy.zeros((dimensions, dimensions), dtype=object) - for i in range(dimensions): - for j in range(dimensions): - tau[i,j] = cse(mu * cse(du[i,j] + du[j,i] - - 2/self.dimensions * delta(i,j) * numpy.trace(du)), - "tau_%d%d" % (i, j)) + # tau = np.zeros((dimensions, dimensions), dtype=object) + # for i in range(dimensions): + # for j in range(dimensions): + # tau[i, j] = sym.cse( + # mu * sym.cse( + # du[i, j] + du[j, i] + # - 2/self.dimensions * delta(i, j) * np.trace(du)), + # "tau_%d%d" % (i, j)) - return tau + # return tau - # }}} + # # }}} # }}} @@ -438,21 +383,21 @@ class GasDynamicsOperator(TimeDependentOperator): # {{{ heat conduction - def heat_conduction_coefficient(self, to_quad_op): - mu = self.get_mu(self.state(), to_quad_op) - if self.prandtl is None or numpy.isinf(self.prandtl): - return 0 + # def heat_conduction_coefficient(self, to_quad_op): + # mu = self.get_mu(self.state(), to_quad_op) + # if self.prandtl is None or np.isinf(self.prandtl): + # return 0 - eos = self.equation_of_state - return (mu / self.prandtl) * (eos.gamma / (eos.gamma-1)) + # eos = self.equation_of_state + # return (mu / self.prandtl) * (eos.gamma / (eos.gamma-1)) - def heat_conduction_grad(self, to_quad_op): - grad_p_over_rho = self.grad_of_state_func( - lambda state: self.p(state)/self.rho(state), - "p_over_rho") + # def heat_conduction_grad(self, to_quad_op): + # grad_p_over_rho = self.grad_of_state_func( + # lambda state: self.p(state)/self.rho(state), + # "p_over_rho") - return (self.heat_conduction_coefficient(to_quad_op) - * to_quad_op(grad_p_over_rho)) + # return (self.heat_conduction_coefficient(to_quad_op) + # * to_quad_op(grad_p_over_rho)) # }}} @@ -461,18 +406,19 @@ class GasDynamicsOperator(TimeDependentOperator): def flux(self, q): from pytools import delta - return [ # one entry for each flux direction - cse(join_fields( + return [ + # one entry for each flux direction + sym.cse(join_fields( # flux rho self.rho_u(q)[i], # flux E - cse(self.e(q)+self.cse_p(q))*self.cse_u(q)[i], + sym.cse(self.e(q)+self.cse_p(q))*self.cse_u(q)[i], # flux rho_u make_obj_array([ - self.rho_u(q)[i]*self.cse_u(q)[j] - + delta(i,j) * self.cse_p(q) + self.rho_u(q)[i]*self.cse_u(q)[j] + + delta(i, j) * self.cse_p(q) for j in range(self.dimensions) ]) ), "%s_flux" % AXES[i]) @@ -488,9 +434,9 @@ class GasDynamicsOperator(TimeDependentOperator): BC is linearized. """ if state0 is None: - state0 = make_sym_vector(bc_name, self.dimensions+2) + state0 = sym.make_sym_vector(bc_name, self.dimensions+2) - state0 = cse(to_bdry_quad(state0)) + state0 = sym.cse(to_bdry_quad(state0)) rho0 = self.rho(state0) p0 = self.cse_p(state0) @@ -498,7 +444,6 @@ class GasDynamicsOperator(TimeDependentOperator): c0 = (self.equation_of_state.gamma * p0 / rho0)**0.5 - from grudge.symbolic import RestrictToBoundary bdrize_op = RestrictToBoundary(tag) class SingleBCInfo(Record): @@ -508,11 +453,11 @@ class GasDynamicsOperator(TimeDependentOperator): rho0=rho0, p0=p0, u0=u0, c0=c0, # notation: suffix "m" for "minus", i.e. "interior" - drhom=cse(self.rho(cse(to_bdry_quad(bdrize_op(state)))) + drhom=sym.cse(self.rho(sym.cse(to_bdry_quad(bdrize_op(state)))) - rho0, "drhom"), - dumvec=cse(self.cse_u(cse(to_bdry_quad(bdrize_op(state)))) + dumvec=sym.cse(self.cse_u(sym.cse(to_bdry_quad(bdrize_op(state)))) - u0, "dumvec"), - dpm=cse(self.cse_p(cse(to_bdry_quad(bdrize_op(state)))) + dpm=sym.cse(self.cse_p(sym.cse(to_bdry_quad(bdrize_op(state)))) - p0, "dpm")) def outflow_state(self, state): @@ -523,33 +468,39 @@ class GasDynamicsOperator(TimeDependentOperator): # see grudge/doc/maxima/euler.mac return join_fields( # bc rho - cse(bc.rho0 - + bc.drhom + numpy.dot(normal, bc.dumvec)*bc.rho0/(2*bc.c0) + sym.cse(bc.rho0 + + bc.drhom + np.dot(normal, bc.dumvec)*bc.rho0/(2*bc.c0) - bc.dpm/(2*bc.c0*bc.c0), "bc_rho_outflow"), # bc p - cse(bc.p0 - + bc.c0*bc.rho0*numpy.dot(normal, bc.dumvec)/2 + bc.dpm/2, "bc_p_outflow"), + sym.cse(bc.p0 + + bc.c0*bc.rho0*np.dot(normal, bc.dumvec)/2 + bc.dpm/2, "bc_p_outflow"), # bc u - cse(bc.u0 - + bc.dumvec - normal*numpy.dot(normal, bc.dumvec)/2 + sym.cse(bc.u0 + + bc.dumvec - normal*np.dot(normal, bc.dumvec)/2 + bc.dpm*normal/(2*bc.c0*bc.rho0), "bc_u_outflow")) def inflow_state_inner(self, normal, bc, name): # see grudge/doc/maxima/euler.mac return join_fields( # bc rho - cse(bc.rho0 - + numpy.dot(normal, bc.dumvec)*bc.rho0/(2*bc.c0) + bc.dpm/(2*bc.c0*bc.c0), "bc_rho_"+name), + sym.cse( + bc.rho0 + + np.dot(normal, bc.dumvec)*bc.rho0/(2*bc.c0) + + bc.dpm/(2*bc.c0*bc.c0), + "bc_rho_"+name), # bc p - cse(bc.p0 - + bc.c0*bc.rho0*numpy.dot(normal, bc.dumvec)/2 + bc.dpm/2, "bc_p_"+name), + sym.cse(bc.p0 + + bc.c0*bc.rho0*np.dot(normal, bc.dumvec)/2 + bc.dpm/2, "bc_p_"+name), # bc u - cse(bc.u0 - + normal*numpy.dot(normal, bc.dumvec)/2 + bc.dpm*normal/(2*bc.c0*bc.rho0), "bc_u_"+name)) + sym.cse( + bc.u0 + + normal*np.dot(normal, bc.dumvec)/2 + + bc.dpm*normal/(2*bc.c0*bc.rho0), + "bc_u_"+name)) def inflow_state(self, state): from grudge.symbolic import make_normal @@ -558,19 +509,17 @@ class GasDynamicsOperator(TimeDependentOperator): return self.inflow_state_inner(normal, bc, "inflow") def noslip_state(self, state): - from grudge.symbolic import make_normal state0 = join_fields( - make_sym_vector("bc_q_noslip", 2), + sym.make_sym_vector("bc_q_noslip", 2), [0]*self.dimensions) normal = make_normal(self.noslip_tag, self.dimensions) bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, state0) return self.inflow_state_inner(normal, bc, "noslip") def wall_state(self, state): - from grudge.symbolic import RestrictToBoundary bc = RestrictToBoundary(self.wall_tag)(state) wall_rho = self.rho(bc) - wall_e = self.e(bc) # <3 eve + wall_e = self.e(bc) # <3 eve wall_rho_u = self.rho_u(bc) from grudge.symbolic import make_normal @@ -579,7 +528,7 @@ class GasDynamicsOperator(TimeDependentOperator): return join_fields( wall_rho, wall_e, - wall_rho_u - 2*numpy.dot(wall_rho_u, normal) * normal) + wall_rho_u - 2*np.dot(wall_rho_u, normal) * normal) @memoize_method def get_primitive_boundary_conditions(self): @@ -591,15 +540,13 @@ class GasDynamicsOperator(TimeDependentOperator): self.noslip_tag: self.noslip_state(state) } - @memoize_method def get_conservative_boundary_conditions(self): state = self.state() - from grudge.symbolic import RestrictToBoundary return { self.supersonic_inflow_tag: - make_sym_vector("bc_q_supersonic_in", self.dimensions+2), + sym.make_sym_vector("bc_q_supersonic_in", self.dimensions+2), self.supersonic_outflow_tag: RestrictToBoundary(self.supersonic_outflow_tag)( (state)), @@ -659,14 +606,13 @@ class GasDynamicsOperator(TimeDependentOperator): from pymbolic.primitives import Subscript, Variable if isinstance(expr, Subscript): - assert (isinstance(expr.aggregate, Variable) + assert (isinstance(expr.aggregate, Variable) and expr.aggregate.name == "q") return cbstate[expr.index] - elif isinstance(expr, Variable) and expr.name =="sensor": - from grudge.symbolic import RestrictToBoundary - result = RestrictToBoundary(tag)(self.sensor()) - return cse(to_bdry_quad(result), "bdry_sensor") + # elif isinstance(expr, Variable) and expr.name =="sensor": + # result = RestrictToBoundary(tag)(self.sensor()) + # return sym.cse(to_bdry_quad(result), "bdry_sensor") from grudge.symbolic import SubstitutionMapper return SubstitutionMapper(subst_func)(expr) @@ -674,68 +620,63 @@ class GasDynamicsOperator(TimeDependentOperator): # }}} # {{{ second order part - def div(self, vol_operand, int_face_operand): - from grudge.second_order import SecondDerivativeTarget - div_tgt = SecondDerivativeTarget( - self.dimensions, strong_form=False, - operand=vol_operand, - int_flux_operand=int_face_operand) - - self.second_order_scheme.div(div_tgt, - bc_getter=self.get_boundary_condition_for, - dirichlet_tags=list(self.get_boundary_tags()), - neumann_tags=[]) - return div_tgt.minv_all - - def make_second_order_part(self): - state = self.state() - faceq_state = self.faceq_state() - volq_state = self.volq_state() - - volq_tau_mat = self.tau(to_vol_quad, state) - faceq_tau_mat = self.tau(to_int_face_quad, state) - - return join_fields( - 0, - self.div( - numpy.sum(volq_tau_mat*self.cse_u(volq_state), axis=1) - + self.heat_conduction_grad(to_vol_quad) - , - numpy.sum(faceq_tau_mat*self.cse_u(faceq_state), axis=1) - + self.heat_conduction_grad(to_int_face_quad) - , - ), - [ - self.div(volq_tau_mat[i], faceq_tau_mat[i]) - for i in range(self.dimensions)] - ) + # def div(self, vol_operand, int_face_operand): + # from grudge.second_order import SecondDerivativeTarget + # div_tgt = SecondDerivativeTarget( + # self.dimensions, strong_form=False, + # operand=vol_operand, + # int_flux_operand=int_face_operand) + + # self.second_order_scheme.div(div_tgt, + # bc_getter=self.get_boundary_condition_for, + # dirichlet_tags=list(self.get_boundary_tags()), + # neumann_tags=[]) + + # return div_tgt.minv_all + + # def make_second_order_part(self): + # state = self.state() + # faceq_state = self.faceq_state() + # volq_state = self.volq_state() + + # volq_tau_mat = self.tau(to_vol_quad, state) + # faceq_tau_mat = self.tau(to_int_face_quad, state) + + # return join_fields( + # 0, + # self.div( + # np.sum(volq_tau_mat*self.cse_u(volq_state), axis=1) + # + self.heat_conduction_grad(to_vol_quad) + # , + # np.sum(faceq_tau_mat*self.cse_u(faceq_state), axis=1) + # + self.heat_conduction_grad(to_int_face_quad) + # , + # ), + # [ + # self.div(volq_tau_mat[i], faceq_tau_mat[i]) + # for i in range(self.dimensions)] + # ) # }}} # {{{ operator template --------------------------------------------------- - def make_extra_terms(self): - return 0 def sym_operator(self, sensor_scaling=None, viscosity_only=False): - u = self.cse_u - rho = self.cse_rho - rho_u = self.rho_u - p = self.p - e = self.e - # {{{ artificial diffusion - def make_artificial_diffusion(): - if self.artificial_viscosity_mode not in ["diffusion"]: - return 0 - dq = self.grad_of_state() + # def make_artificial_diffusion(): + # if self.artificial_viscosity_mode not in ["diffusion"]: + # return 0 + + # dq = self.grad_of_state() + + # return make_obj_array([ + # self.div( + # to_vol_quad(self.sensor())*to_vol_quad(dq[i]), + # to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) + # for i in range(dq.shape[0])]) - return make_obj_array([ - self.div( - to_vol_quad(self.sensor())*to_vol_quad(dq[i]), - to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) - for i in range(dq.shape[0])]) # }}} # {{{ state setup @@ -743,25 +684,20 @@ class GasDynamicsOperator(TimeDependentOperator): volq_flux = self.flux(self.volq_state()) faceq_flux = self.flux(self.faceq_state()) - from grudge.symbolic.primitives import FunctionSymbol - sqrt = FunctionSymbol("sqrt") - speed = self.characteristic_velocity_optemplate(self.state()) - has_viscosity = not is_zero(self.get_mu(self.state(), to_quad_op=None)) + # from grudge.tools import is_zero + # has_viscosity = not is_zero(self.get_mu(self.state(), to_quad_op=None)) # }}} - # {{{ operator assembly ----------------------------------------------- + # {{{ operator assembly from grudge.flux.tools import make_lax_friedrichs_flux - from grudge.symbolic.operators import InverseMassOperator - from grudge.symbolic.tools import make_stiffness_t - - primitive_bcs_as_quad_conservative = dict( - (tag, self.primitive_to_conservative(to_bdry_quad(bc))) - for tag, bc in - six.iteritems(self.get_primitive_boundary_conditions())) + # primitive_bcs_as_quad_conservative = dict( + # (tag, self.primitive_to_conservative(to_bdry_quad(bc))) + # for tag, bc in + # six.iteritems(self.get_primitive_boundary_conditions())) def get_bc_tuple(tag): state = self.state() @@ -769,10 +705,10 @@ class GasDynamicsOperator(TimeDependentOperator): self.get_boundary_condition_for(tag, s_i) for s_i in state]) return tag, bc, self.flux(bc) - first_order_part = InverseMassOperator()( - numpy.dot(make_stiffness_t(self.dimensions), volq_flux) + first_order_part = sym.InverseMassOperator()( + np.dot(sym.StiffnessTOperator(self.dimensions), volq_flux) - make_lax_friedrichs_flux( - wave_speed=cse(to_int_face_quad(speed), "emax_c"), + wave_speed=sym.cse(to_int_face_quad(speed), "emax_c"), state=self.faceq_state(), fluxes=faceq_flux, bdry_tags_states_and_fluxes=[ @@ -783,17 +719,10 @@ class GasDynamicsOperator(TimeDependentOperator): first_order_part = 0*first_order_part result = join_fields( - first_order_part - + self.make_second_order_part() - + make_artificial_diffusion() - + self.make_extra_terms(), - speed) - - if self.source is not None: - result = result + join_fields( - make_sym_vector("source_vect", len(self.state())), - # extra field for speed - 0) + first_order_part + # + self.make_second_order_part() + # + make_artificial_diffusion() + ) return result @@ -801,19 +730,9 @@ class GasDynamicsOperator(TimeDependentOperator): # }}} - # {{{ operator binding ---------------------------------------------------- - def bind(self, discr, sensor=None, sensor_scaling=None, viscosity_only=False): - if (sensor is None and - self.artificial_viscosity_mode is not None): - raise ValueError("must specify a sensor if using " - "artificial viscosity") - - bound_op = discr.compile(self.sym_operator( - sensor_scaling=sensor_scaling, - viscosity_only=False)) - - from grudge.mesh import check_bc_coverage - check_bc_coverage(discr.mesh, [ + def check_bc_coverage(self, mesh): + from meshmode.mesh import check_bc_coverage + check_bc_coverage(mesh, [ self.inflow_tag, self.outflow_tag, self.noslip_tag, @@ -822,43 +741,11 @@ class GasDynamicsOperator(TimeDependentOperator): self.supersonic_outflow_tag, ]) - if self.mu == 0 and not discr.get_boundary(self.noslip_tag).is_empty(): - raise RuntimeError("no-slip BCs only make sense for " - "viscous problems") - - def rhs(t, q): - extra_kwargs = {} - if self.source is not None: - extra_kwargs["source_vect"] = self.source.volume_interpolant( - t, q, discr) - - if sensor is not None: - extra_kwargs["sensor"] = sensor(q) - - opt_result = bound_op(q=q, - bc_q_in=self.bc_inflow.boundary_interpolant( - t, discr, self.inflow_tag), - bc_q_out=self.bc_outflow.boundary_interpolant( - t, discr, self.outflow_tag), - bc_q_noslip=self.bc_noslip.boundary_interpolant( - t, discr, self.noslip_tag), - bc_q_supersonic_in=self.bc_supersonic_inflow - .boundary_interpolant(t, discr, - self.supersonic_inflow_tag), - **extra_kwargs - ) - - max_speed = opt_result[-1] - ode_rhs = opt_result[:-1] - return ode_rhs, discr.nodewise_max(max_speed) - - return rhs - - # }}} - - # {{{ timestep estimation ------------------------------------------------- + # if self.mu == 0 and not discr.get_boundary(self.noslip_tag).is_empty(): + # raise RuntimeError("no-slip BCs only make sense for " + # "viscous problems") - def estimate_timestep(self, discr, + def estimate_timestep(self, discr, stepper=None, stepper_class=None, stepper_args=None, t=None, max_eigenvalue=None): u"""Estimate the largest stable timestep, given a time stepper @@ -876,46 +763,5 @@ class GasDynamicsOperator(TimeDependentOperator): return rk4_dt * approximate_rk4_relative_imag_stability_region( stepper, stepper_class, stepper_args) - # }}} - - - - -# {{{ limiter (unfinished, deprecated) -class SlopeLimiter1NEuler: - def __init__(self, discr, gamma, dimensions, op): - """Construct a limiter from Jan's book page 225 - """ - self.discr = discr - self.gamma=gamma - self.dimensions=dimensions - self.op=op - - from grudge.symbolic.operators import AveragingOperator - self.get_average = AveragingOperator().bind(discr) - - def __call__(self, fields): - from grudge.tools import join_fields - - #get conserved fields - rho=self.op.rho(fields) - e=self.op.e(fields) - rho_velocity=self.op.rho_u(fields) - - #get primitive fields - #to do - - #reset field values to cell average - rhoLim=self.get_average(rho) - eLim=self.get_average(e) - temp=join_fields([self.get_average(rho_vel) - for rho_vel in rho_velocity]) - - #should do for primitive fields too - - return join_fields(rhoLim, eLim, temp) - -# }}} - # vim: foldmethod=marker diff --git a/grudge/models/gas_dynamics/solutions.py b/grudge/models/gas_dynamics/solutions.py new file mode 100644 index 00000000..7a26c25d --- /dev/null +++ b/grudge/models/gas_dynamics/solutions.py @@ -0,0 +1,88 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2019 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +import numpy as np # noqa +import numpy.linalg as la # noqa +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from grudge import sym + +__doc__ = """ +.. autoclass:: IsentropicVortex2D +""" + + +# {{{ IsentropicVortex + +class IsentropicVortex2D(object): + """2D Isentropic vortex analytic solution for the Euler equations. Cf. + [zhou03]_, [jsh09]_. + + .. [zhou03] Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159 + # https://doi.org/10.1016/S0021-9991(03)00206-7 + + .. [jsh09] Jan S. Hesthaven, Tim Warburton. Nodal DG Methods, p. 209 + """ + def __init__(self, beta=5, gamma=1.4, center=None, velocity=None): + self.beta = beta + self.gamma = gamma + + self.center = np.array([5, 0]) + self.velocity = np.array([1, 0]) + self.final_time = 0.5 + + self.mu = 0 + self.prandtl = 0.72 + self.spec_gas_const = 287.1 + + def __call__(self, t=sym.var("t")): + nodes = sym.nodes(2) + + vortex_loc = self.center + t*self.velocity + + # coordinates relative to vortex center + x_rel = nodes[0] - vortex_loc[0] + y_rel = nodes[1] - vortex_loc[1] + + r = sym.sqrt(x_rel**2+y_rel**2) + expterm = self.beta*sym.exp(1-r**2) + u = self.velocity[0] - expterm*y_rel/(2*np.pi) + v = self.velocity[1] + expterm*x_rel/(2*np.pi) + rho = ( + 1-(self.gamma-1)/(16*self.gamma*np.pi**2)*expterm**2 + )**(1/(self.gamma-1)) + p = rho**self.gamma + + e = p/(self.gamma-1) + rho/2*(u**2+v**2) + + return sym.join_fields(rho, e, rho*u, rho*v) + +# }}} + +# vim: foldmethod=marker diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index c59f1518..cef35ac8 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -670,6 +670,8 @@ class TracePair: return 0.5*(self.int + self.ext) +# FIXME: qtag should be removed +# only usage site so far is variable-coeff advection def int_tpair(expression, qtag=None): i = _sym().interp("vol", "int_faces")(expression) e = cse(_sym().OppositeInteriorFaceSwap()(i)) diff --git a/grudge/tools.py b/grudge/tools.py index 8a0f5e81..a1410ab9 100644 --- a/grudge/tools.py +++ b/grudge/tools.py @@ -197,3 +197,44 @@ class OrderedSet(MutableSet): return set(self) == set(other) # }}} + + +# {{{ lax_friedrichs_flux + +def lax_friedrichs_flux(dim, dd, strong, wave_speed, state, flux_func): + """ + :arg strong: a flag for whether this flux is for the + "strong-form" version of DG + :arg fluxes: a callable that takes a state vector and returns + a list of fluxes, one per spatial dimension + """ + from grudge import sym + + normal = sym.normal(dd, dim) + + from pymbolic.primitives import Max + penalty = Max( + (wave_speed.int, wave_speed.ext) + )*(state.ext-state.int) + + int_flux = flux_func(state.int) + ext_flux = flux_func(state.ext) + + if not strong: + return 0.5*( + sum( + n_i*(int_flux_i + ext_flux_i) + for n_i, int_flux_i, ext_flux_i + in zip(normal, int_flux, ext_flux)) + - penalty) + else: + return 0.5*( + sum( + n_i*(int_flux_i - ext_flux_i) + for n_i, int_flux_i, ext_flux_i + in zip(normal, int_flux, ext_flux)) + + penalty) + +# }}} + +# vim: foldmethod=marker diff --git a/setup.cfg b/setup.cfg index b785a55c..f2926605 100644 --- a/setup.cfg +++ b/setup.cfg @@ -2,7 +2,6 @@ ignore = E126,E127,E128,E123,E226,E241,E242,E265,W503,E402 max-line-length=85 exclude= - grudge/models/gas_dynamics, grudge/models/em.py, grudge/models/burgers.py, grudge/models/pml.py, diff --git a/test/test_gas_dynamics.py b/test/test_gas_dynamics.py new file mode 100644 index 00000000..e60772f9 --- /dev/null +++ b/test/test_gas_dynamics.py @@ -0,0 +1,110 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2019 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +import numpy as np # noqa +import numpy.linalg as la # noqa +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 bind, DGDiscretizationWithBoundaries + + +@pytest.mark.parametrize("order", [3, 4, 5]) +def test_isentropic_vortex(ctx_factory, order, visualize=False): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + # from pytools.convergence import EOCRecorder + # eoc_rec = EOCRecorder() + + dims = 3 + ns = [9] + for n in ns: + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + (0, -5), + (10, 5), + n=(n,)*dims) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + from grudge.models.gas_dynamics.solutions import IsentropicVortex2D + vortex = IsentropicVortex2D() + + bound_sym_sol = bind(discr, vortex()) + + dt = 1e-2 + + fields = 0 + + def rhs(t, w): + return 0 + + final_t = 1 + + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("w", dt, fields, rhs) + + step = 0 + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + assert event.component_id == "w" + # esc = event.state_component + + step += 1 + + if visualize and step % 10 == 0: + print(step) + sol_data = bound_sym_sol(queue, t=event.t) + + vis.write_vtk_file("vortex-%s-%04d.vtu" % (n, step), [ + ("rho", sol_data[0]), + ("e", sol_data[1]), + ("rho_u", sol_data[2:]), + ]) + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: foldmethod=marker -- GitLab From eb2f8811ab80391aec6d307fcd088691b1a212e5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 2 Jun 2019 01:07:10 -0500 Subject: [PATCH 2/2] Get rid of unported LBM model --- grudge/models/gas_dynamics/lbm.py | 213 ------------------------------ 1 file changed, 213 deletions(-) delete mode 100644 grudge/models/gas_dynamics/lbm.py diff --git a/grudge/models/gas_dynamics/lbm.py b/grudge/models/gas_dynamics/lbm.py deleted file mode 100644 index 4eb65a61..00000000 --- a/grudge/models/gas_dynamics/lbm.py +++ /dev/null @@ -1,213 +0,0 @@ -# -*- coding: utf8 -*- -"""Lattice-Boltzmann operator.""" - -from __future__ import division -from __future__ import absolute_import -from six.moves import range -from six.moves import zip - -__copyright__ = "Copyright (C) 2011 Andreas Kloeckner" - -__license__ = """ -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in -all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -THE SOFTWARE. -""" - -import numpy as np -import numpy.linalg as la -from grudge.models import HyperbolicOperator -from pytools.obj_array import make_obj_array - - - - - -class LBMMethodBase(object): - def __len__(self): - return len(self.direction_vectors) - - def find_opposites(self): - self.opposites = np.zeros(len(self)) - - for alpha in range(len(self)): - if self.opposites[alpha]: - continue - - found = False - for alpha_2 in range(alpha, len(self)): - if la.norm( - self.direction_vectors[alpha] - + self.direction_vectors[alpha_2]) < 1e-12: - self.opposites[alpha] = alpha_2 - self.opposites[alpha_2] = alpha - found = True - - if not found: - raise RuntimeError( - "direction %s had no opposite" - % self.direction_vectors[alpha]) - - - - -class D2Q9LBMMethod(LBMMethodBase): - def __init__(self): - self.dimensions = 2 - - alphas = np.arange(0, 9) - thetas = (alphas-1)*np.pi/2 - thetas[5:9] += np.pi/4 - - direction_vectors = np.vstack([ - np.cos(thetas), np.sin(thetas)]).T - - direction_vectors[0] *= 0 - direction_vectors[5:9] *= np.sqrt(2) - - direction_vectors[np.abs(direction_vectors) < 1e-12] = 0 - - self.direction_vectors = direction_vectors - - self.weights = np.array([4/9] + [1/9]*4 + [1/36]*4) - - self.speed_of_sound = 1/np.sqrt(3) - self.find_opposites() - - def f_equilibrium(self, rho, alpha, u): - e_alpha = self.direction_vectors[alpha] - c_s = self.speed_of_sound - return self.weights[alpha]*rho*( - 1 - + np.dot(e_alpha, u)/c_s**2 - + 1/2*np.dot(e_alpha, u)**2/c_s**4 - - 1/2*np.dot(u, u)/c_s**2) - - - - -class LatticeBoltzmannOperator(HyperbolicOperator): - def __init__(self, method, lbm_delta_t, nu, flux_type="upwind"): - self.method = method - self.lbm_delta_t = lbm_delta_t - self.nu = nu - - self.flux_type = flux_type - - @property - def tau(self): - return (self.nu - / - (self.lbm_delta_t*self.method.speed_of_sound**2)) - - def get_advection_flux(self, velocity): - from grudge.flux import make_normal, FluxScalarPlaceholder - from pymbolic.primitives import IfPositive - - u = FluxScalarPlaceholder(0) - normal = make_normal(self.method.dimensions) - - if self.flux_type == "central": - return u.avg*np.dot(normal, velocity) - elif self.flux_type == "lf": - return u.avg*np.dot(normal, velocity) \ - + 0.5*la.norm(v)*(u.int - u.ext) - elif self.flux_type == "upwind": - return (np.dot(normal, velocity)* - IfPositive(np.dot(normal, velocity), - u.int, # outflow - u.ext, # inflow - )) - else: - raise ValueError("invalid flux type") - - def get_advection_op(self, q, velocity): - from grudge.symbolic import ( - BoundaryPair, - get_flux_operator, - make_stiffness_t, - InverseMassOperator) - - stiff_t = make_stiffness_t(self.method.dimensions) - - flux_op = get_flux_operator(self.get_advection_flux(velocity)) - return InverseMassOperator()( - np.dot(velocity, stiff_t*q) - flux_op(q)) - - def f_bar(self): - from grudge.symbolic import make_sym_vector - return make_sym_vector("f_bar", len(self.method)) - - def rho(self, f_bar): - return sum(f_bar) - - def rho_u(self, f_bar): - return sum( - dv_i * field_i - for dv_i, field_i in - zip(self.method.direction_vectors, f_bar)) - - def stream_rhs(self, f_bar): - return make_obj_array([ - self.get_advection_op(f_bar_alpha, e_alpha) - for e_alpha, f_bar_alpha in - zip(self.method.direction_vectors, f_bar)]) - - def collision_update(self, f_bar): - from grudge.symbolic.primitives import make_common_subexpression as cse - rho = cse(self.rho(f_bar), "rho") - rho_u = self.rho_u(f_bar) - u = cse(rho_u/rho, "u") - - f_eq_func = self.method.f_equilibrium - f_eq = make_obj_array([ - f_eq_func(rho, alpha, u) for alpha in range(len(self.method))]) - - return f_bar - 1/(self.tau+1/2)*(f_bar - f_eq) - - def bind_rhs(self, discr): - compiled_sym_operator = discr.compile( - self.stream_rhs(self.f_bar())) - - #from grudge.mesh import check_bc_coverage, BTAG_ALL - #check_bc_coverage(discr.mesh, [BTAG_ALL]) - - def rhs(t, f_bar): - return compiled_sym_operator(f_bar=f_bar) - - return rhs - - def bind(self, discr, what): - f_bar_sym = self.f_bar() - - from grudge.symbolic.mappers.type_inference import ( - type_info, NodalRepresentation) - - type_hints = dict( - (f_bar_i, type_info.VolumeVector(NodalRepresentation())) - for f_bar_i in f_bar_sym) - - compiled_sym_operator = discr.compile(what(f_bar_sym), type_hints=type_hints) - - def rhs(f_bar): - return compiled_sym_operator(f_bar=f_bar) - - return rhs - - def max_eigenvalue(self, t=None, fields=None, discr=None): - return max( - la.norm(v) for v in self.method.direction_vectors) -- GitLab