diff --git a/doc/misc.rst b/doc/misc.rst index 3767b3477..3b1d68687 100644 --- a/doc/misc.rst +++ b/doc/misc.rst @@ -71,3 +71,7 @@ References `(DOI) `__ .. [Ihme_2014] Yu Lv and Matthias Ihme (2014) Journal of Computationsl Physics 270 105 \ `(DOI) `__ +.. [Persson_2012] P. Persson and J. Peraire, AIAA 44 \ + `(DOI) `__ +.. [Woodward_1984] Woodward and Colella, Journal of Computational Physics, 54 \ + `(DOI) `__ \ No newline at end of file diff --git a/doc/operators/artificial_viscosity.rst b/doc/operators/artificial_viscosity.rst new file mode 100644 index 000000000..9b583612a --- /dev/null +++ b/doc/operators/artificial_viscosity.rst @@ -0,0 +1,4 @@ +Artificial Viscosity +==================== + +.. automodule:: mirgecom.artificial_viscosity diff --git a/doc/operators/operators.rst b/doc/operators/operators.rst index c0b40c7fe..0907c286d 100644 --- a/doc/operators/operators.rst +++ b/doc/operators/operators.rst @@ -6,4 +6,5 @@ Operators wave-eq diffusion + artificial_viscosity gas-dynamics diff --git a/examples/doublemach-mpi.py b/examples/doublemach-mpi.py new file mode 100644 index 000000000..b144fe3b9 --- /dev/null +++ b/examples/doublemach-mpi.py @@ -0,0 +1,259 @@ +"""Demonstrate double mach reflection.""" + +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + +__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 logging +import pyopencl as cl +import pyopencl.tools as cl_tools +from functools import partial + +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.dof_desc import DTAG_BOUNDARY +from grudge.eager import EagerDGDiscretization +from grudge.shortcuts import make_visualizer + + +from mirgecom.euler import euler_operator +from mirgecom.artificial_viscosity import ( + av_operator, + smoothness_indicator +) +from mirgecom.simutil import ( + inviscid_sim_timestep, + sim_checkpoint, + generate_and_distribute_mesh, + ExactSolutionMismatch, +) +from mirgecom.io import make_init_message +from mirgecom.mpi import mpi_entry_point + +from mirgecom.integrators import rk4_step +from mirgecom.steppers import advance_state +from mirgecom.boundary import AdiabaticSlipBoundary, PrescribedBoundary +from mirgecom.initializers import DoubleMachReflection +from mirgecom.eos import IdealSingleGas + +logger = logging.getLogger(__name__) + + +def get_doublemach_mesh(): + """Generate or import a grid using `gmsh`. + + Input required: + doubleMach.msh (read existing mesh) + + This routine will generate a new grid if it does + not find the grid file (doubleMach.msh). + """ + from meshmode.mesh.io import ( + read_gmsh, + generate_gmsh, + ScriptSource, + ) + import os + meshfile = "doubleMach.msh" + if not os.path.exists(meshfile): + mesh = generate_gmsh( + ScriptSource(""" + x0=1.0/6.0; + setsize=0.025; + Point(1) = {0, 0, 0, setsize}; + Point(2) = {x0,0, 0, setsize}; + Point(3) = {4, 0, 0, setsize}; + Point(4) = {4, 1, 0, setsize}; + Point(5) = {0, 1, 0, setsize}; + Line(1) = {1, 2}; + Line(2) = {2, 3}; + Line(5) = {3, 4}; + Line(6) = {4, 5}; + Line(7) = {5, 1}; + Line Loop(8) = {-5, -6, -7, -1, -2}; + Plane Surface(8) = {8}; + Physical Surface('domain') = {8}; + Physical Curve('ic1') = {6}; + Physical Curve('ic2') = {7}; + Physical Curve('ic3') = {1}; + Physical Curve('wall') = {2}; + Physical Curve('out') = {5}; + """, "geo"), force_ambient_dim=2, dimensions=2, target_unit="M") + else: + mesh = read_gmsh(meshfile, force_ambient_dim=2) + + return mesh + + +@mpi_entry_point +def main(ctx_factory=cl.create_some_context): + """Drive the example.""" + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext( + queue, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)) + ) + + dim = 2 + order = 3 + # Too many steps for CI + # t_final = 1.0e-2 + t_final = 1.0e-3 + current_cfl = 0.1 + current_dt = 1.0e-4 + current_t = 0 + eos = IdealSingleGas() + initializer = DoubleMachReflection() + casename = "doubleMach" + + boundaries = { + DTAG_BOUNDARY("ic1"): PrescribedBoundary(initializer), + DTAG_BOUNDARY("ic2"): PrescribedBoundary(initializer), + DTAG_BOUNDARY("ic3"): PrescribedBoundary(initializer), + DTAG_BOUNDARY("wall"): AdiabaticSlipBoundary(), + DTAG_BOUNDARY("out"): AdiabaticSlipBoundary(), + } + constant_cfl = False + nstatus = 10 + nviz = 100 + rank = 0 + checkpoint_t = current_t + current_step = 0 + timestepper = rk4_step + + s0 = -6.0 + kappa = 1.0 + alpha = 2.0e-2 + from mpi4py import MPI + + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + gen_grid = partial(get_doublemach_mesh) + + local_mesh, global_nelements = generate_and_distribute_mesh(comm, gen_grid) + + local_nelements = local_mesh.nelements + + discr = EagerDGDiscretization(actx, local_mesh, order=order, + mpi_communicator=comm) + + nodes = thaw(actx, discr.nodes()) + current_state = initializer(nodes) + + visualizer = make_visualizer(discr, + discr.order if discr.dim == 2 else discr.order) + initname = initializer.__class__.__name__ + eosname = eos.__class__.__name__ + init_message = make_init_message( + dim=dim, + order=order, + nelements=local_nelements, + global_nelements=global_nelements, + dt=current_dt, + t_final=t_final, + nstatus=nstatus, + nviz=nviz, + cfl=current_cfl, + constant_cfl=constant_cfl, + initname=initname, + eosname=eosname, + casename=casename, + ) + if rank == 0: + logger.info(init_message) + + get_timestep = partial( + inviscid_sim_timestep, + discr=discr, + t=current_t, + dt=current_dt, + cfl=current_cfl, + eos=eos, + t_final=t_final, + constant_cfl=constant_cfl, + ) + + def my_rhs(t, state): + return euler_operator( + discr, cv=state, t=t, boundaries=boundaries, eos=eos + ) + av_operator( + discr, q=state.join(), boundaries=boundaries, + boundary_kwargs={"t": t, "eos": eos}, alpha=alpha, + s0=s0, kappa=kappa + ) + + def my_checkpoint(step, t, dt, state): + tagged_cells = smoothness_indicator(discr, state.mass, s0=s0, kappa=kappa) + viz_fields = [("tagged cells", tagged_cells)] + return sim_checkpoint( + discr, + visualizer, + eos, + cv=state, + vizname=casename, + step=step, + t=t, + dt=dt, + nstatus=nstatus, + nviz=nviz, + constant_cfl=constant_cfl, + comm=comm, + viz_fields=viz_fields, + overwrite=True, + ) + + try: + (current_step, current_t, current_state) = advance_state( + rhs=my_rhs, + timestepper=timestepper, + checkpoint=my_checkpoint, + get_timestep=get_timestep, + state=current_state, + t=current_t, + t_final=t_final, + ) + except ExactSolutionMismatch as ex: + current_step = ex.step + current_t = ex.t + current_state = ex.state + + # if current_t != checkpoint_t: + if rank == 0: + logger.info("Checkpointing final state ...") + my_checkpoint( + current_step, + t=current_t, + dt=(current_t - checkpoint_t), + state=current_state, + ) + + if current_t - t_final < 0: + raise ValueError("Simulation exited abnormally") + + +if __name__ == "__main__": + logging.basicConfig(format="%(message)s", level=logging.INFO) + main() + +# vim: foldmethod=marker diff --git a/examples/heat-source-mpi.py b/examples/heat-source-mpi.py index 196654a61..b683498cd 100644 --- a/examples/heat-source-mpi.py +++ b/examples/heat-source-mpi.py @@ -44,6 +44,7 @@ @mpi_entry_point def main(): + """Drive the example.""" cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue, diff --git a/mirgecom/artificial_viscosity.py b/mirgecom/artificial_viscosity.py new file mode 100644 index 000000000..923f7df03 --- /dev/null +++ b/mirgecom/artificial_viscosity.py @@ -0,0 +1,346 @@ +r""":mod:`mirgecom.artificial_viscosity` Artificial viscosity for hyperbolic systems. + +Consider the following system of conservation laws of the form: + +.. math:: + + \partial_t \mathbf{Q} + \nabla\cdot\mathbf{F} = 0, + +where $\mathbf{Q}$ is the state vector and $\mathbf{F}$ is the vector of +fluxes. This module applies an artificial viscosity term by augmenting +the governing equations in the following way: + +.. math:: + + \partial_t \mathbf{Q} + \nabla\cdot\mathbf{F} = + \nabla\cdot{\varepsilon\nabla\mathbf{Q}}, + +where $\varepsilon$ is the artificial viscosity coefficient. +To evalutate the second order derivative numerically, the problem +is recast as a set of first order problems: + +.. math:: + + \partial_t{\mathbf{Q}} + \nabla\cdot\mathbf{F} &= \nabla\cdot\mathbf{R} \\ + \mathbf{R} &= \varepsilon\nabla\mathbf{Q} + +where $\mathbf{R}$ is an auxiliary variable, and the artificial viscosity +coefficient, $\varepsilon$, is spatially dependent and calculated using the +smoothness indicator of [Persson_2012]_: + +.. math:: + + S_e = \frac{\langle u_{N_p} - u_{N_{p-1}}, u_{N_p} - + u_{N_{p-1}}\rangle_e}{\langle u_{N_p}, u_{N_p} \rangle_e} + +where: + +- $S_e$ is the smoothness indicator +- $u_{N_p}$ is the solution in modal basis at the current polynomial order +- $u_{N_{p-1}}$ is the truncated modal represention to the polynomial order $p-1$ +- The $L_2$ inner product on an element is denoted $\langle \cdot,\cdot \rangle_e$ + +The elementwise viscosity is then calculated: + +.. math:: + + \varepsilon_e = + \begin{cases} + 0, & s_e < s_0 - \kappa \\ + \frac{1}{2}\left( 1 + \sin \frac{\pi(s_e - s_0)}{2 \kappa} \right ), + & s_0-\kappa \le s_e \le s_0+\kappa \\ + 1, & s_e > s_0+\kappa + \end{cases} + +where: + +- $\varepsilon_e$ is the element viscosity +- $s_e = \log_{10}{S_e} \sim 1/p^4$ is the smoothness indicator +- $s_0$ is a reference smoothness value +- $\kappa$ controls the width of the transition between 0 to 1 + +Smoothness Indicator Evaluation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: smoothness_indicator + +AV RHS Evaluation +^^^^^^^^^^^^^^^^^ + +.. autofunction:: av_operator +""" + +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + +__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 + +from pytools import memoize_in, keyed_memoize_in +from pytools.obj_array import obj_array_vectorize +from meshmode.dof_array import thaw, DOFArray +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.eager import interior_trace_pair, cross_rank_trace_pairs +from grudge.symbolic.primitives import TracePair +from grudge.dof_desc import DD_VOLUME_MODAL, DD_VOLUME + + +# FIXME: Remove when get_array_container_context is added to meshmode +def _get_actx(obj): + if isinstance(obj, TracePair): + return _get_actx(obj.int) + if isinstance(obj, np.ndarray): + return _get_actx(obj[0]) + elif isinstance(obj, DOFArray): + return obj.array_context + else: + raise ValueError("Unknown type; can't retrieve array context.") + + +# Tweak the behavior of np.outer to return a lower-dimensional object if either/both +# of the arguments are scalars (np.outer always returns a matrix) +def _outer(a, b): + if isinstance(a, np.ndarray) and isinstance(b, np.ndarray): + return np.outer(a, b) + else: + return a*b + + +def _facial_flux_q(discr, q_tpair): + """Compute facial flux for each scalar component of Q.""" + actx = _get_actx(q_tpair) + + normal = thaw(actx, discr.normal(q_tpair.dd)) + + # This uses a central scalar flux along nhat: + # flux = 1/2 * (Q- + Q+) * nhat + flux_out = _outer(q_tpair.avg, normal) + + return discr.project(q_tpair.dd, "all_faces", flux_out) + + +def _facial_flux_r(discr, r_tpair): + """Compute facial flux for vector component of grad(Q).""" + actx = _get_actx(r_tpair) + + normal = thaw(actx, discr.normal(r_tpair.dd)) + + # This uses a central vector flux along nhat: + # flux = 1/2 * (grad(Q)- + grad(Q)+) .dot. nhat + flux_out = r_tpair.avg @ normal + + return discr.project(r_tpair.dd, "all_faces", flux_out) + + +def av_operator(discr, boundaries, q, alpha, boundary_kwargs=None, **kwargs): + r"""Compute the artificial viscosity right-hand-side. + + Computes the the right-hand-side term for artificial viscosity. + + .. math:: + + \mbox{RHS}_{\mbox{av}} = \nabla\cdot{\varepsilon\nabla\mathbf{Q}} + + Parameters + ---------- + q: :class:`~meshmode.dof_array.DOFArray` or :class:`~numpy.ndarray` + + Single :class:`~meshmode.dof_array.DOFArray` for a scalar or an object array + (:class:`~numpy.ndarray`) for a vector of + :class:`~meshmode.dof_array.DOFArray` on which to operate. + + When used with fluid solvers, *q* is expected to be the fluid state array + of the canonical conserved variables (mass, energy, momentum) + for the fluid along with a vector of species masses for multi-component + fluids. + + boundaries: float + + Dictionary of boundary functions, one for each valid boundary tag + + alpha: float + + The maximum artificial viscosity coefficient to be applied + + boundary_kwargs: :class:`dict` + + dictionary of extra arguments to pass through to the boundary conditions + + Returns + ------- + numpy.ndarray + + The artificial viscosity operator applied to *q*. + """ + if boundary_kwargs is None: + boundary_kwargs = dict() + + # Get smoothness indicator based on first component + indicator_field = q[0] if isinstance(q, np.ndarray) else q + indicator = smoothness_indicator(discr, indicator_field, **kwargs) + + # R=Grad(Q) volume part + if isinstance(q, np.ndarray): + grad_q_vol = np.stack(obj_array_vectorize(discr.weak_grad, q), axis=0) + else: + grad_q_vol = discr.weak_grad(q) + + # R=Grad(Q) Q flux over interior faces + q_flux_int = _facial_flux_q(discr, q_tpair=interior_trace_pair(discr, q)) + # R=Grad(Q) Q flux interior faces on partition boundaries + q_flux_pb = sum(_facial_flux_q(discr, q_tpair=pb_tpair) + for pb_tpair in cross_rank_trace_pairs(discr, q)) + # R=Grad(Q) Q flux domain boundary part (i.e. BCs) + q_flux_db = 0 + for btag in boundaries: + q_tpair = TracePair( + btag, + interior=discr.project("vol", btag, q), + exterior=boundaries[btag].exterior_q(discr, btag=btag, q=q, + **boundary_kwargs)) + q_flux_db = q_flux_db + _facial_flux_q(discr, q_tpair=q_tpair) + # Total Q flux across element boundaries + q_bnd_flux = q_flux_int + q_flux_pb + q_flux_db + + # Compute R + r = discr.inverse_mass( + -alpha * indicator * (grad_q_vol - discr.face_mass(q_bnd_flux)) + ) + + # RHS_av = div(R) volume part + div_r_vol = discr.weak_div(r) + # RHS_av = div(R): grad(Q) flux interior faces part + r_flux_int = _facial_flux_r(discr, r_tpair=interior_trace_pair(discr, r)) + # RHS_av = div(R): grad(Q) flux interior faces on the partition boundaries + r_flux_pb = sum(_facial_flux_r(discr, r_tpair=pb_tpair) + for pb_tpair in cross_rank_trace_pairs(discr, r)) + # RHS_av = div(R): grad(Q) flux domain boundary part (BCs) + r_flux_db = 0 + for btag in boundaries: + r_tpair = TracePair( + btag, + interior=discr.project("vol", btag, r), + exterior=boundaries[btag].exterior_grad_q(discr, btag=btag, grad_q=r, + **boundary_kwargs)) + r_flux_db = r_flux_db + _facial_flux_r(discr, r_tpair=r_tpair) + # Total grad(Q) flux element boundaries + r_flux_bnd = r_flux_int + r_flux_pb + r_flux_db + + # Return the AV RHS term + return discr.inverse_mass(-div_r_vol + discr.face_mass(r_flux_bnd)) + + +def artificial_viscosity(discr, t, eos, boundaries, q, alpha, **kwargs): + """Interface :function:`av_operator` with backwards-compatible API.""" + from warnings import warn + warn("Do not call artificial_viscosity; it is now called av_operator. This" + "function will disappear in 2021", DeprecationWarning, stacklevel=2) + return av_operator(discr=discr, boundaries=boundaries, + boundary_kwargs={"t": t, "eos": eos}, q=q, alpha=alpha, **kwargs) + + +def smoothness_indicator(discr, u, kappa=1.0, s0=-6.0): + r"""Calculate the smoothness indicator. + + Parameters + ---------- + u: meshmode.dof_array.DOFArray + The field that is used to calculate the smoothness indicator. + + kappa + An optional argument that controls the width of the transition from 0 to 1. + s0 + An optional argument that sets the smoothness level to limit + on. Values in the range $(-\infty,0]$ are allowed, where $-\infty$ results in + all cells being tagged and 0 results in none. + + Returns + ------- + meshmode.dof_array.DOFArray + The elementwise constant values between 0 and 1 which indicate the smoothness + of a given element. + """ + if not isinstance(u, DOFArray): + raise ValueError("u argument must be a DOFArray.") + + actx = u.array_context + + @memoize_in(actx, (smoothness_indicator, "smooth_comp_knl")) + def indicator_prg(): + """Compute the smoothness indicator for all elements.""" + from meshmode.array_context import make_loopy_program + return make_loopy_program([ + "{[iel]: 0 <= iel < nelements}", + "{[idof]: 0 <= idof < ndiscr_nodes_in}", + "{[jdof]: 0 <= jdof < ndiscr_nodes_in}", + "{[kdof]: 0 <= kdof < ndiscr_nodes_in}" + ], + """ + result[iel,idof] = sum(kdof, vec[iel, kdof] \ + * vec[iel, kdof] \ + * modes_active_flag[kdof]) / \ + sum(jdof, vec[iel, jdof] \ + * vec[iel, jdof] \ + + 1.0e-12 / ndiscr_nodes_in) + """, + name="smooth_comp", + ) + + @keyed_memoize_in(actx, (smoothness_indicator, + "highest_mode"), + lambda grp: grp.discretization_key()) + def highest_mode(grp): + return actx.from_numpy( + np.asarray([1 if sum(mode_id) == grp.order + else 0 + for mode_id in grp.mode_ids()]) + ) + + # Convert to modal solution representation + modal_map = discr.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) + uhat = modal_map(u) + + # Compute smoothness indicator value + indicator = DOFArray( + actx, + data=tuple( + actx.call_loopy( + indicator_prg(), + vec=uhat[grp.index], + modes_active_flag=highest_mode(grp))["result"] + for grp in discr.discr_from_dd("vol").groups + ) + ) + indicator = actx.np.log10(indicator + 1.0e-12) + + # Compute artificial viscosity percentage based on indicator and set parameters + yesnol = actx.np.greater(indicator, (s0 - kappa)) + yesnou = actx.np.greater(indicator, (s0 + kappa)) + sin_indicator = actx.np.where( + yesnol, + 0.5 * (1.0 + actx.np.sin(np.pi * (indicator - s0) / (2.0 * kappa))), + 0.0 * indicator, + ) + indicator = actx.np.where(yesnou, 1.0 + 0.0 * indicator, sin_indicator) + + return indicator diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 19d37ca87..ac12af480 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -61,13 +61,22 @@ def __init__(self, userfunc): def boundary_pair(self, discr, cv, btag, **kwargs): """Get the interior and exterior solution on the boundary.""" + ext_soln = self.exterior_q(discr, cv, btag, **kwargs) + int_soln = discr.project("vol", btag, cv) + return TracePair(btag, interior=int_soln, exterior=ext_soln) + + def exterior_q(self, discr, cv, btag, **kwargs): + """Get the exterior solution on the boundary.""" actx = cv.array_context boundary_discr = discr.discr_from_dd(btag) nodes = thaw(actx, boundary_discr.nodes()) ext_soln = self._userfunc(nodes, **kwargs) - int_soln = discr.project("vol", btag, cv) - return TracePair(btag, interior=int_soln, exterior=ext_soln) + return ext_soln + + def exterior_grad_q(self, discr, grad_cv, btag, **kwargs): + """Get the exterior solution on the boundary.""" + return discr.project("vol", btag, grad_cv) class DummyBoundary: @@ -78,9 +87,18 @@ class DummyBoundary: def boundary_pair(self, discr, cv, btag, **kwargs): """Get the interior and exterior solution on the boundary.""" - dir_soln = discr.project("vol", btag, cv) + dir_soln = self.exterior_q(discr, cv, btag, **kwargs) return TracePair(btag, interior=dir_soln, exterior=dir_soln) + def exterior_q(self, discr, cv, btag, **kwargs): + """Get the exterior solution on the boundary.""" + dir_soln = discr.project("vol", btag, cv) + return dir_soln + + def exterior_grad_q(self, discr, grad_cv, btag, **kwargs): + """Get the grad_q on the exterior of the boundary.""" + return discr.project("vol", btag, grad_cv) + class AdiabaticSlipBoundary: r"""Boundary condition implementing inviscid slip boundary. @@ -101,7 +119,14 @@ class AdiabaticSlipBoundary: """ def boundary_pair(self, discr, cv, btag, **kwargs): - """Get the interior and exterior solution on the boundary. + """Get the interior and exterior solution on the boundary.""" + bndry_soln = self.exterior_q(discr, cv, btag, **kwargs) + int_soln = discr.project("vol", btag, cv) + + return TracePair(btag, interior=int_soln, exterior=bndry_soln) + + def exterior_q(self, discr, cv, btag, **kwargs): + """Get the exterior solution on the boundary. The exterior solution is set such that there will be vanishing flux through the boundary, preserving mass, momentum (magnitude) and @@ -130,9 +155,31 @@ def boundary_pair(self, discr, cv, btag, **kwargs): ext_mom = int_cv.momentum - 2.0 * wnorm_mom # prescribed ext momentum # Form the external boundary solution with the new momentum - bndry_cv = make_conserved(dim=dim, mass=int_cv.mass, - energy=int_cv.energy, - momentum=ext_mom, - species_mass=int_cv.species_mass) + bndry_soln = make_conserved(dim=dim, mass=int_cv.mass, + energy=int_cv.energy, + momentum=ext_mom, + species_mass=int_cv.species_mass) + + return bndry_soln + + def exterior_grad_q(self, discr, grad_cv, btag, **kwargs): + """Get the exterior grad(Q) on the boundary.""" + # Grab some boundary-relevant data + num_equations, dim = grad_cv.shape + actx = grad_cv.mass[0].array_context + + # Grab a unit normal to the boundary + normal = thaw(actx, discr.normal(btag)) + + # Get the interior soln + gradq_comp = discr.project("vol", btag, grad_cv) + + # Subtract 2*wall-normal component of q + # to enforce q=0 on the wall + s_mom_normcomp = np.outer(normal, np.dot(gradq_comp.momentum, normal)) + s_mom_flux = gradq_comp.momentum - 2*s_mom_normcomp - return TracePair(btag, interior=int_cv, exterior=bndry_cv) + # flip components to set a neumann condition + return make_conserved(dim, mass=-gradq_comp.mass, energy=-gradq_comp.energy, + momentum=-s_mom_flux, + species_mass=-gradq_comp.species_mass) diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index 1174a2de4..a4ae5e65b 100644 --- a/mirgecom/initializers.py +++ b/mirgecom/initializers.py @@ -5,6 +5,7 @@ ^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: Vortex2D .. autoclass:: SodShock1D +.. autoclass:: DoubleMachReflection .. autoclass:: Lump .. autoclass:: MulticomponentLump .. autoclass:: Uniform @@ -266,6 +267,122 @@ def __call__(self, x_vec, *, t=0, eos=None): momentum=mom) +class DoubleMachReflection: + r"""Implement the double shock reflection problem. + + The double shock reflection solution is crafted after [Woodward_1984]_ + and is defined by: + + .. math:: + + {\rho}(x < x_s(y,t)) &= \gamma \rho_j\\ + {\rho}(x > x_s(y,t)) &= \gamma\\ + {\rho}{V_x}(x < x_s(y,t)) &= u_p \cos(\pi/6)\\ + {\rho}{V_x}(x > x_s(y,t)) &= 0\\ + {\rho}{V_y}(x > x_s(y,t)) &= u_p \sin(\pi/6)\\ + {\rho}{V_y}(x > x_s(y,t)) &= 0\\ + {\rho}E(x < x_s(y,t)) &= (\gamma-1)p_j\\ + {\rho}E(x > x_s(y,t)) &= (\gamma-1), + + where the shock position is given, + + .. math:: + + x_s = x_0 + y/\sqrt{3} + 2 u_s t/\sqrt{3} + + and the normal shock jump relations are + + .. math:: + + \rho_j &= \frac{(\gamma + 1) u_s^2}{(\gamma-1) u_s^2 + 2} \\ + p_j &= \frac{2 \gamma u_s^2 - (\gamma - 1)}{\gamma+1} \\ + u_p &= 2 \frac{u_s^2-1}{(\gamma+1) u_s} + + The initial shock location is given by $x_0$ and $u_s$ is the shock speed. + + .. automethod:: __init__ + .. automethod:: __call__ + """ + + def __init__( + self, shock_location=1.0/6.0, shock_speed=4.0 + ): + """Initialize double shock reflection parameters. + + Parameters + ---------- + shock_location: float + initial location of shock + shock_speed: float + shock speed, Mach number + """ + self._shock_location = shock_location + self._shock_speed = shock_speed + + def __call__(self, x_vec, *, t=0, eos=None, **kwargs): + r""" + Create double mach reflection solution at locations *x_vec*. + + At times $t > 0$, calls to this routine create an advanced solution + under the assumption of constant normal shock speed *shock_speed*. + The advanced solution *is not* the exact solution, but is appropriate + for use as an exact boundary solution on the top and upstream (left) + side of the domain. + + Parameters + ---------- + t: float + Time at which to compute the solution + x_vec: numpy.ndarray + Nodal coordinates + eos: :class:`mirgecom.eos.GasEOS` + Equation of state class to be used in construction of soln (if needed) + """ + # Fail if numdim is other than 2 + if(len(x_vec)) != 2: + raise ValueError("Case only defined for 2 dimensions") + if eos is None: + eos = IdealSingleGas() + + gm1 = eos.gamma() - 1.0 + gp1 = eos.gamma() + 1.0 + gmn1 = 1.0 / gm1 + x_rel = x_vec[0] + y_rel = x_vec[1] + actx = x_rel.array_context + + # Normal Shock Relations + shock_speed_2 = self._shock_speed * self._shock_speed + rho_jump = gp1 * shock_speed_2 / (gm1 * shock_speed_2 + 2.) + p_jump = (2. * eos.gamma() * shock_speed_2 - gm1) / gp1 + up = 2. * (shock_speed_2 - 1.) / (gp1 * self._shock_speed) + + rhol = eos.gamma() * rho_jump + rhor = eos.gamma() + ul = up * np.cos(np.pi/6.0) + ur = 0.0 + vl = - up * np.sin(np.pi/6.0) + vr = 0.0 + rhoel = gmn1 * p_jump + rhoer = gmn1 * 1.0 + + xinter = (self._shock_location + y_rel/np.sqrt(3.0) + + 2.0*self._shock_speed*t/np.sqrt(3.0)) + sigma = 0.05 + xtanh = 1.0/sigma*(x_rel-xinter) + mass = rhol/2.0*(actx.np.tanh(-xtanh)+1.0)+rhor/2.0*(actx.np.tanh(xtanh)+1.0) + rhoe = (rhoel/2.0*(actx.np.tanh(-xtanh)+1.0) + + rhoer/2.0*(actx.np.tanh(xtanh)+1.0)) + u = ul/2.0*(actx.np.tanh(-xtanh)+1.0)+ur/2.0*(actx.np.tanh(xtanh)+1.0) + v = vl/2.0*(actx.np.tanh(-xtanh)+1.0)+vr/2.0*(actx.np.tanh(xtanh)+1.0) + + vel = make_obj_array([u, v]) + mom = mass * vel + energy = rhoe + .5*mass*np.dot(vel, vel) + + return make_conserved(dim=2, mass=mass, energy=energy, momentum=mom) + + class Lump: r"""Solution initializer for N-dimensional Gaussian lump of mass. diff --git a/requirements.txt b/requirements.txt index 5a5d4b277..b60e76ad2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,6 +6,7 @@ pyvisfile pymetis importlib-resources psutil +gmsh # The following packages will be git cloned by emirge: --editable git+https://github.com/inducer/pymbolic.git#egg=pymbolic diff --git a/test/test_av.py b/test/test_av.py new file mode 100644 index 000000000..0390bc960 --- /dev/null +++ b/test/test_av.py @@ -0,0 +1,197 @@ +"""Test the artificial viscosity functions.""" + +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + + +__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 logging +import numpy as np +import pyopencl as cl +import pytest +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw +from meshmode.mesh import BTAG_ALL +from mirgecom.artificial_viscosity import ( + av_operator, + smoothness_indicator +) +from mirgecom.boundary import DummyBoundary +from grudge.eager import EagerDGDiscretization +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests, +) + +logger = logging.getLogger(__name__) + + +@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize("order", [1, 5]) +def test_tag_cells(ctx_factory, dim, order): + """Test tag_cells. + + Tests that the cells tagging properly tags cells + given prescribed solutions. + """ + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + nel_1d = 2 + tolerance = 1.e-16 + + def norm_indicator(expected, discr, soln, **kwargs): + return(discr.norm(expected-smoothness_indicator(discr, soln, **kwargs), + np.inf)) + + from meshmode.mesh.generation import generate_regular_rect_mesh + + mesh = generate_regular_rect_mesh( + a=(-1.0, )*dim, b=(1.0, )*dim, n=(nel_1d, ) * dim + ) + + discr = EagerDGDiscretization(actx, mesh, order=order) + nodes = thaw(actx, discr.nodes()) + nele = mesh.nelements + zeros = 0.0*nodes[0] + + # test jump discontinuity + soln = actx.np.where(nodes[0] > 0.0+zeros, 1.0+zeros, zeros) + err = norm_indicator(1.0, discr, soln) + + assert err < tolerance, "Jump discontinuity should trigger indicator (1.0)" + + # get meshmode polynomials + group = discr.discr_from_dd("vol").groups[0] + basis = group.basis() # only one group + unit_nodes = group.unit_nodes + modes = group.mode_ids() + order = group.order + + # loop over modes and check smoothness + for i, mode in enumerate(modes): + ele_soln = basis[i](unit_nodes) + soln[0].set(np.tile(ele_soln, (nele, 1))) + if sum(mode) == order: + expected = 1.0 + else: + expected = 0.0 + err = norm_indicator(expected, discr, soln) + assert err < tolerance, "Only highest modes should trigger indicator (1.0)" + + # Test s0 + s0 = -1. + eps = 1.0e-6 + + phi_n_p = np.sqrt(np.power(10, s0)) + phi_n_pm1 = np.sqrt(1. - np.power(10, s0)) + + # pick a polynomial of order n_p, n_p-1 + n_p = np.array(np.nonzero((np.sum(modes, axis=1) == order))).flat[0] + n_pm1 = np.array(np.nonzero((np.sum(modes, axis=1) == order-1))).flat[0] + + # create test soln perturbed around + # Solution above s0 + ele_soln = ((phi_n_p+eps)*basis[n_p](unit_nodes) + + phi_n_pm1*basis[n_pm1](unit_nodes)) + soln[0].set(np.tile(ele_soln, (nele, 1))) + err = norm_indicator(1.0, discr, soln, s0=s0, kappa=0.0) + assert err < tolerance, ( + "A function with an indicator >s0 should trigger indicator") + + # Solution below s0 + ele_soln = ((phi_n_p-eps)*basis[n_p](unit_nodes) + + phi_n_pm1*basis[n_pm1](unit_nodes)) + soln[0].set(np.tile(ele_soln, (nele, 1))) + err = norm_indicator(0.0, discr, soln, s0=s0, kappa=0.0) + assert err < tolerance, ( + "A function with an indicator tolerance, "s_e>s_0-kappa should trigger indicator" + + # lower bound + err = norm_indicator(1.0, discr, soln, s0=s0-(kappa+shift), kappa=kappa) + assert err < tolerance, "s_e>s_0+kappa should fully trigger indicator (1.0)" + err = norm_indicator(1.0, discr, soln, s0=s0-(kappa-shift), kappa=kappa) + assert err > tolerance, "s_e