Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • inducer/grudge
  • jdsteve2/grudge
  • eshoag2/grudge
  • mattwala/grudge
  • kaushikcfd/grudge
  • fikl2/grudge
6 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1417 additions and 1288 deletions
doc/dof_desc.rst 0 → 100644
View file @ 49b9b7c7
Degree of freedom (DOF) descriptions
====================================
.. automodule:: grudge.dof_desc
doc/geometry.rst 0 → 100644
View file @ 49b9b7c7
Metric terms and transformations
================================
.. automodule:: grudge.geometry.metrics
doc/index.rst
View file @ 49b9b7c7
.. grudge documentation master file, created by
sphinx-quickstart on Sun Sep 27 13:08:30 2015.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
Welcome to grudge's documentation!
Welcome to grudge's Documentation!
==================================
Contents:
Here's an example to solve the PDE
.. math::
\begin{cases}
u_t + 2\pi u_x = 0, \\
u(0, t) = -\sin(2\pi t), \\
u(x, 0) = \sin(x),
\end{cases}
on the domain :math:`x \in [0, 2\pi]`. We closely follow Chapter 3 of
[Hesthaven_2008]_.
.. toctree::
:maxdepth: 2
.. literalinclude:: ../examples/hello-grudge.py
:start-after: BEGINEXAMPLE
:end-before: ENDEXAMPLE
Plotting numerical solution ``uh`` in results in
Indices and tables
.. plot:: ../examples/hello-grudge.py
Contents:
.. toctree::
:maxdepth: 2
discretization
dof_desc
geometry
operators
utils
models
references
misc
🚀 Github <https://github.com/inducer/grudge>
💾 Download Releases <https://pypi.org/project/grudge>
Indices and Tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`
doc/misc.rst 0 → 100644
View file @ 49b9b7c7
.. _installation:
Installation
============
Installing :mod:`grudge`
------------------------
This set of instructions is intended for 64-bit Linux computers.
MacOS support is in the works.
#. Make sure your system has the basics to build software.
On Debian derivatives (Ubuntu and many more),
installing ``build-essential`` should do the trick.
Everywhere else, just making sure you have the ``g++`` package should be
enough.
#. Install `miniforge <https://github.com/conda-forge/miniforge>`_::
curl -L -O https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh
# then run
bash ./Miniforge3-*.sh
#. ``export CONDA=/WHERE/YOU/INSTALLED/miniforge3``
If you accepted the default location, this should work:
``export CONDA=$HOME/miniforge3``
#. ``$CONDA/bin/conda create -n dgfem``
#. ``source $CONDA/bin/activate dgfem``
#. ``conda install git pip pocl islpy pyopencl``
#. Type the following command::
hash -r; for i in pymbolic cgen genpy modepy pyvisfile loopy arraycontext meshmode dagrt leap grudge; do python -m pip install --editable "git+https://github.com/inducer/$i.git#egg=$i"; done
.. note::
In each case, you may leave out the ``--editable`` flag if you would not like
an editable copy of the source code checked out in a subfolder.
Next time you want to use `grudge`, just run the following command::
source /WHERE/YOU/INSTALLED/miniforge3/bin/activate dgfem
You may also like to add this to a startup file (like :file:`$HOME/.bashrc`) or create an alias for it.
After this, you should be able to run the `tests <https://gitlab.tiker.net/inducer/grudge/tree/main/test>`_
or `examples <https://gitlab.tiker.net/inducer/grudge/tree/main/examples>`_.
Troubleshooting the Installation
--------------------------------
/usr/bin/ld: cannot find -lstdc++
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Try::
sudo apt-get install libstdc++-6-dev
to install the missing C++ development package.
No CL platforms found/unknown error -1001
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
If you get::
pyopencl.cffi_cl.LogicError: clGetPlatformIDs failed: <unknown error -1001>
try::
conda update ocl-icd pocl
(This indicates that the OpenCL driver loader didn't find any drivers, or the
drivers were themselves missing dependencies.)
Assertion 'error == 0'
~~~~~~~~~~~~~~~~~~~~~~~
If you get::
/opt/conda/conda-bld/home_1484016200338/work/pocl-0.13/lib/CL/devices/common.c:108:
llvm_codegen: Assertion 'error == 0 ' failed. Aborted (core dumped)
then you're likely out of memory.
User-visible Changes
====================
Version 2016.1
--------------
.. note::
This version is currently under development. You can get snapshots from
grudge's `git repository <https://github.com/inducer/grudge>`_
Licensing
=========
:mod:`grudge` is licensed to you under the MIT/X Consortium license:
Copyright (c) 2014-21 Andreas Klöckner and Contributors.
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.
Acknowledgments
===============
Work on grudge was supported in part by
* the Department of Energy, National Nuclear Security Administration,
under Award Number DE-NA0003963,
* the US Navy ONR, under grant number N00014-14-1-0117, and
* the US National Science Foundation under grant numbers CCF-1524433,
and OAC-1931577.
AK also gratefully acknowledges a hardware gift from Nvidia Corporation.
The views and opinions expressed herein do not necessarily reflect those of the
funding agencies.
doc/models.rst 0 → 100644
View file @ 49b9b7c7
Discontinuous Galerkin Models
=============================
Compressible Euler Equations
----------------------------
.. automodule:: grudge.models.euler
doc/operators.rst 0 → 100644
View file @ 49b9b7c7
Discontinuous Galerkin operators
================================
.. automodule:: grudge.op
.. automodule:: grudge.trace_pair
Transferring data between discretizations
=========================================
.. automodule:: grudge.projection
Reductions
==========
.. automodule:: grudge.reductions
doc/references.rst 0 → 100644
View file @ 49b9b7c7
References
==========
..
When adding references here, please use the demonstrated format:
[FirstAuthor_pubyear]
.. [Shewchuk_2002] J. R. Shewchuk (2002),
*What Is a Good Linear Element? - Interpolation, Conditioning, Quality Measures*,
In Proceedings of the 11th International Meshing Roundtable.
`(link) <https://hal.science/hal-04614934/>`__
`(slides) <https://people.eecs.berkeley.edu/~jrs/papers/elemtalk.pdf>`__
.. [Hesthaven_2008] J. S. Hesthaven and T. Warburton (2008),
*Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications*,
Springer.
`(doi) <https://doi.org/10.1007/978-0-387-72067-8>`__
.. [Chan_2016] J. Chan, R. J. Hewett, and T. Warburton (2016),
*Weight-Adjusted Discontinuous Galerkin Methods: Curvilinear Meshes*,
SIAM Journal on Scientific Computing.
`(doi) <https://doi.org/10.1137/16M1089198>`__
doc/upload-docs.sh
View file @ 49b9b7c7
#! /bin/sh
rsync --progress --verbose --archive --delete _build/html/* doc-upload:doc/grudge
rsync --verbose --archive --delete _build/html/ doc-upload:doc/grudge
doc/utils.rst 0 → 100644
View file @ 49b9b7c7
Helper functions
================
Estimating stable time-steps
----------------------------
.. automodule:: grudge.dt_utils
Array contexts
--------------
.. automodule:: grudge.array_context
Miscellaneous tools
-------------------
.. automodule:: grudge.tools
examples/advection/advection.py deleted 100644 → 0
View file @ 519e3e71
# 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 <http://www.gnu.org/licenses/>.
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
import numpy.linalg as la
def main(write_output=True, flux_type_arg="upwind"):
from grudge.tools import mem_checkpoint
from math import sin, cos, pi, sqrt
from math import floor
from grudge.backends import guess_run_context
rcon = guess_run_context()
def f(x):
return sin(pi*x)
def u_analytic(x, el, t):
return f((-numpy.dot(v, x)/norm_v+t*norm_v))
def boundary_tagger(vertices, el, face_nr, all_v):
if numpy.dot(el.face_normals[face_nr], v) < 0:
return ["inflow"]
else:
return ["outflow"]
dim = 2
if dim == 1:
v = numpy.array([1])
if rcon.is_head_rank:
from grudge.mesh.generator import make_uniform_1d_mesh
mesh = make_uniform_1d_mesh(0, 2, 10, periodic=True)
elif dim == 2:
v = numpy.array([2,0])
if rcon.is_head_rank:
from grudge.mesh.generator import make_disk_mesh
mesh = make_disk_mesh(boundary_tagger=boundary_tagger)
elif dim == 3:
v = numpy.array([0,0,1])
if rcon.is_head_rank:
from grudge.mesh.generator import make_cylinder_mesh, make_ball_mesh, make_box_mesh
mesh = make_cylinder_mesh(max_volume=0.04, height=2, boundary_tagger=boundary_tagger,
periodic=False, radial_subdivisions=32)
else:
raise RuntimeError("bad number of dimensions")
norm_v = la.norm(v)
if rcon.is_head_rank:
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
if dim != 1:
mesh_data = mesh_data.reordered_by("cuthill")
discr = rcon.make_discretization(mesh_data, order=4)
vis_discr = discr
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(vis_discr, rcon, "fld")
# operator setup ----------------------------------------------------------
from grudge.data import \
ConstantGivenFunction, \
TimeConstantGivenFunction, \
TimeDependentGivenFunction
from grudge.models.advection import StrongAdvectionOperator, WeakAdvectionOperator
op = WeakAdvectionOperator(v,
inflow_u=TimeDependentGivenFunction(u_analytic),
flux_type=flux_type_arg)
u = discr.interpolate_volume_function(lambda x, el: u_analytic(x, el, 0))
# timestep setup ----------------------------------------------------------
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper()
if rcon.is_head_rank:
print("%d elements" % len(discr.mesh.elements))
# diagnostics setup -------------------------------------------------------
from pytools.log import LogManager, \
add_general_quantities, \
add_simulation_quantities, \
add_run_info
if write_output:
log_file_name = "advection.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"])
# timestep loop -----------------------------------------------------------
rhs = op.bind(discr)
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=3, 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 % 5 == 0 and write_output:
visf = vis.make_file("fld-%04d" % step)
vis.add_data(visf, [
("u", discr.convert_volume(u, kind="numpy")),
], time=t, step=step)
visf.close()
u = stepper(u, t, dt, rhs)
true_u = discr.interpolate_volume_function(lambda x, el: u_analytic(x, el, t))
print(discr.norm(u-true_u))
assert discr.norm(u-true_u) < 1e-2
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
def test_advection():
from pytools.test import mark_test
mark_long = mark_test.long
for flux_type in ["upwind", "central", "lf"]:
yield "advection with %s flux" % flux_type, \
mark_long(main), False, flux_type
examples/advection/surface.py 0 → 100644
View file @ 49b9b7c7
"""Minimal example of a grudge driver for DG on surfaces."""
__copyright__ = """
Copyright (C) 2020 Alexandru Fikl
Copyright (C) 2021 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 os
import numpy as np
import pyopencl as cl
import pyopencl.tools as cl_tools
from arraycontext import flatten
from meshmode.discretization.connection import FACE_RESTR_INTERIOR
from pytools.obj_array import make_obj_array
import grudge.dof_desc as dof_desc
import grudge.geometry as geo
import grudge.op as op
from grudge.array_context import PyOpenCLArrayContext
logger = logging.getLogger(__name__)
# {{{ plotting (keep in sync with `var-velocity.py`)
class Plotter:
def __init__(self, actx, dcoll, order, visualize=True):
self.actx = actx
self.ambient_dim = dcoll.ambient_dim
self.dim = dcoll.dim
self.visualize = visualize
if not self.visualize:
return
if self.ambient_dim == 2:
import matplotlib.pyplot as pt
self.fig = pt.figure(figsize=(8, 8), dpi=300)
x = actx.thaw(dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL).nodes())
self.x = actx.to_numpy(flatten(actx.np.arctan2(x[1], x[0]), self.actx))
elif self.ambient_dim == 3:
from grudge.shortcuts import make_visualizer
self.vis = make_visualizer(dcoll)
else:
raise ValueError("unsupported dimension")
def __call__(self, evt, basename, overwrite=True):
if not self.visualize:
return
if self.ambient_dim == 2:
u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
filename = f"{basename}.png"
if not overwrite and os.path.exists(filename):
from meshmode import FileExistsError
raise FileExistsError(f"output file '{filename}' already exists")
ax = self.fig.gca()
ax.grid()
ax.plot(self.x, u, "-")
ax.plot(self.x, u, "k.")
ax.set_xlabel(r"$\theta$")
ax.set_ylabel("$u$")
ax.set_title(f"t = {evt.t:.2f}")
self.fig.savefig(filename)
self.fig.clf()
elif self.ambient_dim == 3:
self.vis.write_vtk_file(f"{basename}.vtu", [
("u", evt.state_component)
], overwrite=overwrite)
else:
raise ValueError("unsupported dimension")
# }}}
def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
actx = PyOpenCLArrayContext(queue, allocator=allocator)
# {{{ parameters
# sphere radius
radius = 1.0
# sphere resolution
resolution = 64 if dim == 2 else 1
# final time
final_time = np.pi
# flux
flux_type = "lf"
# }}}
# {{{ discretization
if dim == 2:
from meshmode.mesh.generation import ellipse, make_curve_mesh
mesh = make_curve_mesh(
lambda t: radius * ellipse(1.0, t),
np.linspace(0.0, 1.0, resolution + 1),
order)
elif dim == 3:
from meshmode.mesh.generation import generate_icosphere
mesh = generate_icosphere(radius, order=4 * order,
uniform_refinement_rounds=resolution)
else:
raise ValueError("unsupported dimension")
discr_tag_to_group_factory = {}
if use_quad:
qtag = dof_desc.DISCR_TAG_QUAD
else:
qtag = None
from meshmode.discretization.poly_element import (
QuadratureSimplexGroupFactory,
default_simplex_group_factory,
)
discr_tag_to_group_factory[dof_desc.DISCR_TAG_BASE] = \
default_simplex_group_factory(base_dim=dim-1, order=order)
if use_quad:
discr_tag_to_group_factory[qtag] = \
QuadratureSimplexGroupFactory(order=4*order)
from grudge.discretization import make_discretization_collection
dcoll = make_discretization_collection(
actx, mesh,
discr_tag_to_group_factory=discr_tag_to_group_factory
)
volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
logger.info("ndofs: %d", volume_discr.ndofs)
logger.info("nelements: %d", volume_discr.mesh.nelements)
# }}}
# {{{ Surface advection operator
# velocity field
x = actx.thaw(dcoll.nodes())
c = make_obj_array([-x[1], x[0], 0.0])[:dim]
def f_initial_condition(x):
return x[0]
from grudge.models.advection import SurfaceAdvectionOperator
adv_operator = SurfaceAdvectionOperator(
dcoll,
c,
flux_type=flux_type,
quad_tag=qtag
)
u0 = f_initial_condition(x)
def rhs(t, u):
return adv_operator.operator(t, u)
# check velocity is tangential
from grudge.geometry import normal
surf_normal = normal(actx, dcoll, dd=dof_desc.DD_VOLUME_ALL)
error = op.norm(dcoll, c.dot(surf_normal), 2)
logger.info("u_dot_n: %.5e", error)
# }}}
# {{{ time stepping
# FIXME: dt estimate is not necessarily valid for surfaces
dt = actx.to_numpy(
0.45 * adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u0))
nsteps = int(final_time // dt) + 1
logger.info("dt: %.5e", dt)
logger.info("nsteps: %d", nsteps)
from grudge.shortcuts import set_up_rk4
dt_stepper = set_up_rk4("u", float(dt), u0, rhs)
plot = Plotter(actx, dcoll, order, visualize=visualize)
norm_u = actx.to_numpy(op.norm(dcoll, u0, 2))
step = 0
event = dt_stepper.StateComputed(0.0, 0, 0, u0)
plot(event, f"fld-surface-{0:04d}")
if visualize and dim == 3:
from grudge.shortcuts import make_visualizer
vis = make_visualizer(dcoll)
vis.write_vtk_file(
"fld-surface-velocity.vtu",
[
("u", c),
("n", surf_normal)
],
overwrite=True
)
df = dof_desc.as_dofdesc(FACE_RESTR_INTERIOR)
face_discr = dcoll.discr_from_dd(df)
face_normal = geo.normal(actx, dcoll, dd=df)
from meshmode.discretization.visualization import make_visualizer
vis = make_visualizer(actx, face_discr)
vis.write_vtk_file("fld-surface-face-normals.vtu", [
("n", face_normal)
], overwrite=True)
for event in dt_stepper.run(t_end=final_time, max_steps=nsteps + 1):
if not isinstance(event, dt_stepper.StateComputed):
continue
step += 1
if step % 10 == 0:
norm_u = actx.to_numpy(op.norm(dcoll, event.state_component, 2))
plot(event, f"fld-surface-{step:04d}")
logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert norm_u < 3
# }}}
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--dim", choices=[2, 3], default=2, type=int)
parser.add_argument("--order", default=4, type=int)
parser.add_argument("--use-quad", action="store_true")
parser.add_argument("--visualize", action="store_true")
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
dim=args.dim,
order=args.order,
use_quad=args.use_quad,
visualize=args.visualize)
examples/advection/var-velocity.py
View file @ 49b9b7c7
# grudge - the Hybrid'n'Easy DG Environment
# Copyright (C) 2009 Andreas Stock
#
# 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 <http://www.gnu.org/licenses/>.
__copyright__ = """
Copyright (C) 2017 Bogdan Enache
Copyright (C) 2021 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 os
import numpy as np
import pyopencl as cl
import pyopencl.tools as cl_tools
from arraycontext import flatten
from meshmode.mesh import BTAG_ALL
from pytools.obj_array import flat_obj_array
import grudge.dof_desc as dof_desc
import grudge.op as op
from grudge.array_context import PyOpenCLArrayContext
logger = logging.getLogger(__name__)
# {{{ plotting (keep in sync with `weak.py`)
class Plotter:
def __init__(self, actx, dcoll, order, visualize=True, ylim=None):
self.actx = actx
self.dim = dcoll.ambient_dim
self.visualize = visualize
if not self.visualize:
return
if self.dim == 1:
import matplotlib.pyplot as pt
self.fig = pt.figure(figsize=(8, 8), dpi=300)
self.ylim = ylim
volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
self.x = actx.to_numpy(flatten(volume_discr.nodes()[0], self.actx))
else:
from grudge.shortcuts import make_visualizer
self.vis = make_visualizer(dcoll)
def __call__(self, evt, basename, overwrite=True):
if not self.visualize:
return
if self.dim == 1:
u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
filename = f"{basename}.png"
if not overwrite and os.path.exists(filename):
from meshmode import FileExistsError
raise FileExistsError(f"output file '{filename}' already exists")
ax = self.fig.gca()
ax.plot(self.x, u, "-")
ax.plot(self.x, u, "k.")
if self.ylim is not None:
ax.set_ylim(self.ylim)
ax.set_xlabel("$x$")
ax.set_ylabel("$u$")
ax.set_title(f"t = {evt.t:.2f}")
self.fig.savefig(filename)
self.fig.clf()
else:
self.vis.write_vtk_file(f"{basename}.vtu", [
("u", evt.state_component)
], overwrite=overwrite)
# }}}
def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False,
flux_type="upwind"):
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
actx = PyOpenCLArrayContext(queue, allocator=allocator)
# {{{ parameters
# domain [0, d]^dim
d = 1.0
# number of points in each dimension
npoints = 25
# final time
final_time = 1
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
if use_quad:
qtag = dof_desc.DISCR_TAG_QUAD
else:
qtag = None
# }}}
def main(write_output=True, flux_type_arg="central", use_quadrature=True,
final_time=20):
from math import sin, cos, pi, sqrt
# {{{ discretization
from grudge.backends import guess_run_context
rcon = guess_run_context()
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(0,)*dim, b=(d,)*dim,
npoints_per_axis=(npoints,)*dim,
order=order)
# 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)
from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
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}
if use_quad:
discr_tag_to_group_factory = {
qtag: QuadratureSimplexGroupFactory(order=4*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
discr_tag_to_group_factory = {}
from grudge.discretization import make_discretization_collection
logmgr = LogManager(log_file_name, "w", rcon.communicator)
add_run_info(logmgr)
add_general_quantities(logmgr)
add_simulation_quantities(logmgr)
discr.add_instrumentation(logmgr)
dcoll = make_discretization_collection(
actx, mesh, order=order,
discr_tag_to_group_factory=discr_tag_to_group_factory
)
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"))
# {{{ advection operator
logmgr.add_watches(["step.max", "t_sim.max", "l2_u", "t_step.max"])
# gaussian parameters
# Initialize v for data output:
v = op.advec_v.volume_interpolant(0, discr)
def f_halfcircle(x):
source_center = np.array([d/2, d/2, d/2])[:dim]
dist = x - source_center
return (
(0.5+0.5*actx.np.tanh(500*(-np.dot(dist, dist) + 0.4**2)))
* (0.5+0.5*actx.np.tanh(500*(dist[0]))))
# 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))
def zero_inflow_bc(dtag, t=0):
dd = dof_desc.as_dofdesc(dtag, qtag)
return dcoll.discr_from_dd(dd).zeros(actx)
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()
from grudge.models.advection import VariableCoefficientAdvectionOperator
x = actx.thaw(dcoll.nodes())
u = stepper(u, t, dt, rhs)
# velocity field
if dim == 1:
c = x
else:
# solid body rotation
c = flat_obj_array(
np.pi * (d/2 - x[1]),
np.pi * (x[0] - d/2),
0
)[:dim]
# 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)
adv_operator = VariableCoefficientAdvectionOperator(
dcoll,
c,
inflow_u=lambda t: zero_inflow_bc(BTAG_ALL, t),
quad_tag=qtag,
flux_type=flux_type
)
assert discr.norm(u) < 10
u = f_halfcircle(x)
finally:
if write_output:
vis.close()
def rhs(t, u):
return adv_operator.operator(t, u)
logmgr.close()
discr.close()
dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u))
logger.info("Timestep size: %g", dt)
# }}}
if __name__ == "__main__":
main()
# {{{ time stepping
from grudge.shortcuts import set_up_rk4
dt_stepper = set_up_rk4("u", float(dt), u, rhs)
plot = Plotter(actx, dcoll, order, visualize=visualize,
ylim=[-0.1, 1.1])
step = 0
for event in dt_stepper.run(t_end=final_time):
if not isinstance(event, dt_stepper.StateComputed):
continue
if step % 10 == 0:
norm_u = actx.to_numpy(op.norm(dcoll, event.state_component, 2))
plot(event, f"fld-var-velocity-{step:04d}")
# entry points for py.test ----------------------------------------------------
def test_var_velocity_advection():
from pytools.test import mark_test
mark_long = mark_test.long
logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert norm_u < 1
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"
step += 1
yield descr, mark_long(main), False, flux_type, use_quadrature, 1
# }}}
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--dim", default=2, type=int)
parser.add_argument("--order", default=4, type=int)
parser.add_argument("--use-quad", action="store_true")
parser.add_argument("--visualize", action="store_true")
parser.add_argument("--flux", default="upwind",
help="'central' or 'upwind'. Run with central to observe aliasing "
"instability. Add --use-quad to fix that instability.")
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
dim=args.dim,
order=args.order,
use_quad=args.use_quad,
visualize=args.visualize,
flux_type=args.flux)
examples/advection/weak.py 0 → 100644
View file @ 49b9b7c7
__copyright__ = """
Copyright (C) 2007 Andreas Kloeckner
Copyright (C) 2021 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 os
import numpy as np
import numpy.linalg as la
import pyopencl as cl
import pyopencl.tools as cl_tools
from arraycontext import flatten
from meshmode.mesh import BTAG_ALL
import grudge.dof_desc as dof_desc
import grudge.op as op
from grudge.array_context import PyOpenCLArrayContext
logger = logging.getLogger(__name__)
# {{{ plotting (keep in sync with `var-velocity.py`)
class Plotter:
def __init__(self, actx, dcoll, order, visualize=True, ylim=None):
self.actx = actx
self.dim = dcoll.ambient_dim
self.visualize = visualize
if not self.visualize:
return
if self.dim == 1:
import matplotlib.pyplot as pt
self.fig = pt.figure(figsize=(8, 8), dpi=300)
self.ylim = ylim
volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
self.x = actx.to_numpy(flatten(volume_discr.nodes()[0], self.actx))
else:
from grudge.shortcuts import make_visualizer
self.vis = make_visualizer(dcoll)
def __call__(self, evt, basename, overwrite=True):
if not self.visualize:
return
if self.dim == 1:
u = self.actx.to_numpy(flatten(evt.state_component, self.actx))
filename = f"{basename}.png"
if not overwrite and os.path.exists(filename):
from meshmode import FileExistsError
raise FileExistsError(f"output file '{filename}' already exists")
ax = self.fig.gca()
ax.plot(self.x, u, "-")
ax.plot(self.x, u, "k.")
if self.ylim is not None:
ax.set_ylim(self.ylim)
ax.set_xlabel("$x$")
ax.set_ylabel("$u$")
ax.set_title(f"t = {evt.t:.2f}")
self.fig.savefig(filename)
self.fig.clf()
else:
self.vis.write_vtk_file(f"{basename}.vtu", [
("u", evt.state_component)
], overwrite=overwrite)
# }}}
def main(ctx_factory, dim=2, order=4, visualize=False):
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
actx = PyOpenCLArrayContext(queue, allocator=allocator)
# {{{ parameters
# domain [-d/2, d/2]^dim
d = 1.0
# number of points in each dimension
npoints = 20
# final time
final_time = 1.0
# velocity field
c = np.array([0.5] * dim)
norm_c = la.norm(c)
# flux
flux_type = "central"
# }}}
# {{{ discretization
from meshmode.mesh.generation import generate_box_mesh
mesh = generate_box_mesh(
[np.linspace(-d/2, d/2, npoints) for _ in range(dim)],
order=order)
from grudge.discretization import make_discretization_collection
dcoll = make_discretization_collection(actx, mesh, order=order)
# }}}
# {{{ weak advection operator
def f(x):
return actx.np.sin(3 * x)
def u_analytic(x, t=0):
return f(-np.dot(c, x) / norm_c + t * norm_c)
from grudge.models.advection import WeakAdvectionOperator
adv_operator = WeakAdvectionOperator(
dcoll,
c,
inflow_u=lambda t: u_analytic(
actx.thaw(dcoll.nodes(dd=BTAG_ALL)),
t=t
),
flux_type=flux_type
)
nodes = actx.thaw(dcoll.nodes())
u = u_analytic(nodes, t=0)
def rhs(t, u):
return adv_operator.operator(t, u)
dt = actx.to_numpy(adv_operator.estimate_rk4_timestep(actx, dcoll, fields=u))
logger.info("Timestep size: %g", dt)
# }}}
# {{{ time stepping
from grudge.shortcuts import set_up_rk4
dt_stepper = set_up_rk4("u", float(dt), u, rhs)
plot = Plotter(actx, dcoll, order, visualize=visualize,
ylim=[-1.1, 1.1])
step = 0
norm_u = 0.0
for event in dt_stepper.run(t_end=final_time):
if not isinstance(event, dt_stepper.StateComputed):
continue
if step % 10 == 0:
norm_u = actx.to_numpy(op.norm(dcoll, event.state_component, 2))
plot(event, f"fld-weak-{step:04d}")
step += 1
logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u)
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert norm_u < 1
# }}}
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--dim", default=2, type=int)
parser.add_argument("--order", default=4, type=int)
parser.add_argument("--visualize", action="store_true")
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
dim=args.dim,
order=args.order,
visualize=args.visualize)
examples/burgers/burgers.py deleted 100644 → 0
View file @ 519e3e71
# 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 <http://www.gnu.org/licenses/>.
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
import numpy.linalg as la
from math import sin, cos, pi, sqrt
from pytools.test import mark_test
class ExactTestCase:
a = 0
b = 150
final_time = 1000
def u0(self, x):
return self.u_exact(x, 0)
def u_exact(self, x, t):
# CAUTION: This gets the shock speed wrong as soon as the pulse
# starts interacting with itself.
def f(x, shock_loc):
if x < (t-40)/4:
return 1/4
else:
if t < 40:
if x < (3*t)/4:
return (x+15)/(t+20)
elif x < (t+80)/4:
return (x-30)/(t-40)
else:
return 1/4
else:
if x < shock_loc:
return (x+15)/(t+20)
else:
return 1/4
from math import sqrt
shock_loc = 30*sqrt(2*t+40)/sqrt(120) + t/4 - 10
shock_win = (shock_loc + 20) // self.b
x += shock_win * 150
x -= 20
return max(f(x, shock_loc), f(x-self.b, shock_loc-self.b))
class OffCenterMigratingTestCase:
a = -pi
b = pi
final_time = 10
def u0(self, x):
return -0.4+sin(x+0.1)
class CenteredStationaryTestCase:
# does funny things to P-P
a = -pi
b = pi
final_time = 10
def u0(self, x):
return -sin(x)
class OffCenterStationaryTestCase:
# does funny things to P-P
a = -pi
b = pi
final_time = 10
def u0(self, x):
return -sin(x+0.3)
def main(write_output=True, flux_type_arg="upwind",
#case = CenteredStationaryTestCase(),
#case = OffCenterStationaryTestCase(),
#case = OffCenterMigratingTestCase(),
case = ExactTestCase(),
):
from grudge.backends import guess_run_context
rcon = guess_run_context()
order = 3
if rcon.is_head_rank:
if True:
from grudge.mesh.generator import make_uniform_1d_mesh
mesh = make_uniform_1d_mesh(case.a, case.b, 20, periodic=True)
else:
from grudge.mesh.generator import make_rect_mesh
print((pi*2)/(11*5*2))
mesh = make_rect_mesh((-pi, -1), (pi, 1),
periodicity=(True, True),
subdivisions=(11,5),
max_area=(pi*2)/(11*5*2)
)
if rcon.is_head_rank:
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
discr = rcon.make_discretization(mesh_data, order=order,
quad_min_degrees={"quad": 3*order})
if write_output:
from grudge.visualization import VtkVisualizer
vis = VtkVisualizer(discr, rcon, "fld")
# operator setup ----------------------------------------------------------
from grudge.second_order import IPDGSecondDerivative
from grudge.models.burgers import BurgersOperator
op = BurgersOperator(mesh.dimensions,
viscosity_scheme=IPDGSecondDerivative())
if rcon.is_head_rank:
print("%d elements" % len(discr.mesh.elements))
# exact solution ----------------------------------------------------------
import pymbolic
var = pymbolic.var
u = discr.interpolate_volume_function(lambda x, el: case.u0(x[0]))
# diagnostics setup -------------------------------------------------------
from pytools.log import LogManager, \
add_general_quantities, \
add_simulation_quantities, \
add_run_info
if write_output:
log_file_name = "burgers.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)
from grudge.log import LpNorm
u_getter = lambda: u
logmgr.add_quantity(LpNorm(u_getter, discr, p=1, name="l1_u"))
logmgr.add_watches(["step.max", "t_sim.max", "l1_u", "t_step.max"])
# timestep loop -----------------------------------------------------------
rhs = op.bind(discr)
from grudge.timestep.runge_kutta import ODE45TimeStepper, LSRK4TimeStepper
stepper = ODE45TimeStepper()
stepper.add_instrumentation(logmgr)
try:
from grudge.timestep import times_and_steps
# for visc=0.01
#stab_fac = 0.1 # RK4
#stab_fac = 1.6 # dumka3(3), central
#stab_fac = 3 # dumka3(4), central
#stab_fac = 0.01 # RK4
stab_fac = 0.2 # dumka3(3), central
#stab_fac = 3 # dumka3(4), central
dt = stab_fac*op.estimate_timestep(discr,
stepper=LSRK4TimeStepper(), t=0, fields=u)
step_it = times_and_steps(
final_time=case.final_time, logmgr=logmgr, max_dt_getter=lambda t: dt)
from grudge.symbolic import InverseVandermondeOperator
inv_vdm = InverseVandermondeOperator().bind(discr)
for step, t, dt in step_it:
if step % 3 == 0 and write_output:
if hasattr(case, "u_exact"):
extra_fields = [
("u_exact",
discr.interpolate_volume_function(
lambda x, el: case.u_exact(x[0], t)))]
else:
extra_fields = []
visf = vis.make_file("fld-%04d" % step)
vis.add_data(visf, [
("u", u),
] + extra_fields,
time=t,
step=step)
visf.close()
u = stepper(u, t, dt, rhs)
if isinstance(case, ExactTestCase):
assert discr.norm(u, 1) < 50
finally:
if write_output:
vis.close()
logmgr.save()
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
@mark_test.long
def test_stability():
main(write_output=False)
examples/conftest.py deleted 100644 → 0
View file @ 519e3e71
# This file tells py.test to scan for test_xxx() functions in
# any file below here, not just those named test_xxxx.py.
def pytest_collect_file(path, parent):
if "grudge/examples" in str(path.dirpath()) and path.ext == ".py":
return parent.Module(path, parent=parent)
examples/euler/acoustic_pulse.py 0 → 100644
View file @ 49b9b7c7
__copyright__ = """
Copyright (C) 2021 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 pyopencl.tools as cl_tools
from arraycontext import ArrayContext
from meshmode.mesh import BTAG_ALL
from pytools.obj_array import make_obj_array
import grudge.op as op
from grudge.array_context import PyOpenCLArrayContext, PytatoPyOpenCLArrayContext
from grudge.models.euler import ConservedEulerField, EulerOperator, InviscidWallBC
from grudge.shortcuts import rk4_step
logger = logging.getLogger(__name__)
def gaussian_profile(
actx: ArrayContext,
x_vec, t=0, rho0=1.0, rhoamp=1.0, p0=1.0, gamma=1.4,
center=None, velocity=None):
dim = len(x_vec)
if center is None:
center = np.zeros(shape=(dim,))
if velocity is None:
velocity = np.zeros(shape=(dim,))
lump_loc = center + t * velocity
# coordinates relative to lump center
rel_center = make_obj_array(
[x_vec[i] - lump_loc[i] for i in range(dim)]
)
r = actx.np.sqrt(np.dot(rel_center, rel_center))
expterm = rhoamp * actx.np.exp(1 - r ** 2)
mass = expterm + rho0
mom = velocity.astype(object) * mass
energy = (p0 / (gamma - 1.0)) + np.dot(mom, mom) / (2.0 * mass)
return ConservedEulerField(mass=mass, energy=energy, momentum=mom)
def make_pulse(amplitude, r0, w, r):
dim = len(r)
r_0 = np.zeros(dim)
r_0 = r_0 + r0
rel_center = make_obj_array(
[r[i] - r_0[i] for i in range(dim)]
)
actx = r[0].array_context
rms2 = w * w
r2 = np.dot(rel_center, rel_center) / rms2
return amplitude * actx.np.exp(-.5 * r2)
def acoustic_pulse_condition(actx: ArrayContext, x_vec, t=0):
dim = len(x_vec)
vel = np.zeros(shape=(dim,))
orig = np.zeros(shape=(dim,))
uniform_gaussian = gaussian_profile(
actx,
x_vec, t=t, center=orig, velocity=vel, rhoamp=0.0)
amplitude = 1.0
width = 0.1
pulse = make_pulse(amplitude, orig, width, x_vec)
return ConservedEulerField(
mass=uniform_gaussian.mass,
energy=uniform_gaussian.energy + pulse,
momentum=uniform_gaussian.momentum
)
def run_acoustic_pulse(actx,
order=3,
final_time=1,
resolution=16,
overintegration=False,
visualize=False):
# eos-related parameters
gamma = 1.4
# {{{ discretization
from meshmode.mesh.generation import generate_regular_rect_mesh
dim = 2
box_ll = -0.5
box_ur = 0.5
mesh = generate_regular_rect_mesh(
a=(box_ll,)*dim,
b=(box_ur,)*dim,
nelements_per_axis=(resolution,)*dim)
from meshmode.discretization.poly_element import (
QuadratureSimplexGroupFactory,
default_simplex_group_factory,
)
from grudge.discretization import make_discretization_collection
from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD
exp_name = f"fld-acoustic-pulse-N{order}-K{resolution}"
if overintegration:
exp_name += "-overintegrated"
quad_tag = DISCR_TAG_QUAD
else:
quad_tag = None
dcoll = make_discretization_collection(
actx, mesh,
discr_tag_to_group_factory={
DISCR_TAG_BASE: default_simplex_group_factory(
base_dim=mesh.dim, order=order),
DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order)
}
)
# }}}
# {{{ Euler operator
euler_operator = EulerOperator(
dcoll,
bdry_conditions={BTAG_ALL: InviscidWallBC()},
flux_type="lf",
gamma=gamma,
quadrature_tag=quad_tag
)
def rhs(t, q):
return euler_operator.operator(actx, t, q)
compiled_rhs = actx.compile(rhs)
from grudge.dt_utils import h_min_from_volume
cfl = 0.125
cn = 0.5*(order + 1)**2
dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn
fields = acoustic_pulse_condition(actx, actx.thaw(dcoll.nodes()))
logger.info("Timestep size: %g", dt)
# }}}
from grudge.shortcuts import make_visualizer
vis = make_visualizer(dcoll)
# {{{ time stepping
step = 0
t = 0.0
while t < final_time:
if step % 10 == 0:
norm_q = actx.to_numpy(op.norm(dcoll, fields, 2))
logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q)
if visualize:
vis.write_vtk_file(
f"{exp_name}-{step:04d}.vtu",
[
("rho", fields.mass),
("energy", fields.energy),
("momentum", fields.momentum)
]
)
assert norm_q < 5
fields = actx.thaw(actx.freeze(fields))
fields = rk4_step(fields, t, dt, compiled_rhs)
t += dt
step += 1
# }}}
def main(ctx_factory, order=3, final_time=1, resolution=16,
overintegration=False, visualize=False, lazy=False):
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
if lazy:
actx = PytatoPyOpenCLArrayContext(queue, allocator=allocator)
else:
actx = PyOpenCLArrayContext(queue, allocator=allocator)
run_acoustic_pulse(
actx,
order=order,
resolution=resolution,
overintegration=overintegration,
final_time=final_time,
visualize=visualize
)
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--order", default=3, type=int)
parser.add_argument("--tfinal", default=0.1, type=float)
parser.add_argument("--resolution", default=16, type=int)
parser.add_argument("--oi", action="store_true",
help="use overintegration")
parser.add_argument("--visualize", action="store_true",
help="write out vtk output")
parser.add_argument("--lazy", action="store_true",
help="switch to a lazy computation mode")
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
order=args.order,
final_time=args.tfinal,
resolution=args.resolution,
overintegration=args.oi,
visualize=args.visualize,
lazy=args.lazy)
examples/euler/vortex.py 0 → 100644
View file @ 49b9b7c7
__copyright__ = """
Copyright (C) 2021 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
import grudge.op as op
from grudge.array_context import PyOpenCLArrayContext, PytatoPyOpenCLArrayContext
from grudge.models.euler import EulerOperator, vortex_initial_condition
from grudge.shortcuts import rk4_step
logger = logging.getLogger(__name__)
def run_vortex(actx, order=3, resolution=8, final_time=5,
overintegration=False,
flux_type="central",
visualize=False):
logger.info(
"""
Isentropic vortex parameters:\n
order: %s\n
final_time: %s\n
resolution: %s\n
overintegration: %s\n
flux_type: %s\n
visualize: %s
""",
order, final_time, resolution,
overintegration, flux_type, visualize
)
# eos-related parameters
gamma = 1.4
# {{{ discretization
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(0, -5),
b=(20, 5),
nelements_per_axis=(2*resolution, resolution),
periodic=(True, True))
from meshmode.discretization.poly_element import (
QuadratureSimplexGroupFactory,
default_simplex_group_factory,
)
from grudge.discretization import make_discretization_collection
from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD
exp_name = f"fld-vortex-N{order}-K{resolution}-{flux_type}"
if overintegration:
exp_name += "-overintegrated"
quad_tag = DISCR_TAG_QUAD
else:
quad_tag = None
dcoll = make_discretization_collection(
actx, mesh,
discr_tag_to_group_factory={
DISCR_TAG_BASE: default_simplex_group_factory(
base_dim=mesh.dim, order=order),
DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order)
}
)
# }}}
# {{{ Euler operator
euler_operator = EulerOperator(
dcoll,
flux_type=flux_type,
gamma=gamma,
quadrature_tag=quad_tag
)
def rhs(t, q):
return euler_operator.operator(actx, t, q)
compiled_rhs = actx.compile(rhs)
fields = vortex_initial_condition(actx.thaw(dcoll.nodes()))
from grudge.dt_utils import h_min_from_volume
cfl = 0.01
cn = 0.5*(order + 1)**2
dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn
logger.info("Timestep size: %g", dt)
# }}}
from grudge.shortcuts import make_visualizer
vis = make_visualizer(dcoll)
# {{{ time stepping
step = 0
t = 0.0
while t < final_time:
if step % 10 == 0:
norm_q = actx.to_numpy(op.norm(dcoll, fields, 2))
logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q)
if visualize:
vis.write_vtk_file(
f"{exp_name}-{step:04d}.vtu",
[
("rho", fields.mass),
("energy", fields.energy),
("momentum", fields.momentum)
]
)
assert norm_q < 200
fields = actx.thaw(actx.freeze(fields))
fields = rk4_step(fields, t, dt, compiled_rhs)
t += dt
step += 1
# }}}
def main(ctx_factory, order=3, final_time=5, resolution=8,
overintegration=False,
lf_stabilization=False,
visualize=False,
lazy=False):
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
if lazy:
actx = PytatoPyOpenCLArrayContext(queue, allocator=allocator)
else:
actx = PyOpenCLArrayContext(queue, allocator=allocator)
if lf_stabilization:
flux_type = "lf"
else:
flux_type = "central"
run_vortex(
actx,
order=order,
resolution=resolution,
final_time=final_time,
overintegration=overintegration,
flux_type=flux_type,
visualize=visualize
)
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--order", default=3, type=int)
parser.add_argument("--tfinal", default=0.015, type=float)
parser.add_argument("--resolution", default=8, type=int)
parser.add_argument("--oi", action="store_true",
help="use overintegration")
parser.add_argument("--lf", action="store_true",
help="turn on lax-friedrichs dissipation")
parser.add_argument("--visualize", action="store_true",
help="write out vtk output")
parser.add_argument("--lazy", action="store_true",
help="switch to a lazy computation mode")
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
main(cl.create_some_context,
order=args.order,
final_time=args.tfinal,
resolution=args.resolution,
overintegration=args.oi,
lf_stabilization=args.lf,
visualize=args.visualize,
lazy=args.lazy)
examples/gas_dynamics/box-in-box.py deleted 100644 → 0
View file @ 519e3e71
# grudge - the Hybrid'n'Easy DG Environment
# Copyright (C) 2008 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 <http://www.gnu.org/licenses/>.
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
def make_boxmesh():
from meshpy.tet import MeshInfo, build
from meshpy.geometry import GeometryBuilder, Marker, make_box
geob = GeometryBuilder()
box_marker = Marker.FIRST_USER_MARKER
extent_small = 0.1*numpy.ones(3, dtype=numpy.float64)
geob.add_geometry(*make_box(-extent_small, extent_small))
# make small "separator box" for region attribute
geob.add_geometry(
*make_box(
-extent_small*4,
numpy.array([4, 0.4, 0.4], dtype=numpy.float64)))
geob.add_geometry(
*make_box(
numpy.array([-1, -1, -1], dtype=numpy.float64),
numpy.array([5, 1, 1], dtype=numpy.float64)))
mesh_info = MeshInfo()
geob.set(mesh_info)
mesh_info.set_holes([(0, 0, 0)])
# region attributes
mesh_info.regions.resize(1)
mesh_info.regions[0] = (
# point in region
list(extent_small*2) + [
# region number
1,
# max volume in region
#0.0001
0.005
])
mesh = build(mesh_info, max_volume=0.02,
volume_constraints=True, attributes=True)
print("%d elements" % len(mesh.elements))
#mesh.write_vtk("box-in-box.vtk")
#print "done writing"
fvi2fm = mesh.face_vertex_indices_to_face_marker
face_marker_to_tag = {
box_marker: "noslip",
Marker.MINUS_X: "inflow",
Marker.PLUS_X: "outflow",
Marker.MINUS_Y: "inflow",
Marker.PLUS_Y: "inflow",
Marker.PLUS_Z: "inflow",
Marker.MINUS_Z: "inflow"
}
def bdry_tagger(fvi, el, fn, all_v):
face_marker = fvi2fm[fvi]
return [face_marker_to_tag[face_marker]]
from grudge.mesh import make_conformal_mesh
return make_conformal_mesh(
mesh.points, mesh.elements, bdry_tagger)
def main():
from grudge.backends import guess_run_context
rcon = guess_run_context(["cuda"])
if rcon.is_head_rank:
mesh = make_boxmesh()
#from grudge.mesh import make_rect_mesh
#mesh = make_rect_mesh(
# boundary_tagger=lambda fvi, el, fn, all_v: ["inflow"])
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [3]:
from pytools import add_python_path_relative_to_script
add_python_path_relative_to_script("..")
from gas_dynamics_initials import UniformMachFlow
box = UniformMachFlow(angle_of_attack=0)
from grudge.models.gas_dynamics import GasDynamicsOperator
op = GasDynamicsOperator(dimensions=3,
gamma=box.gamma, mu=box.mu,
prandtl=box.prandtl, spec_gas_const=box.spec_gas_const,
bc_inflow=box, bc_outflow=box, bc_noslip=box,
inflow_tag="inflow", outflow_tag="outflow", noslip_tag="noslip")
discr = rcon.make_discretization(mesh_data, order=order,
debug=[
#"cuda_no_plan",
#"cuda_dump_kernels",
#"dump_dataflow_graph",
#"dump_optemplate_stages",
#"dump_dataflow_graph",
#"print_op_code",
"cuda_no_plan_el_local",
],
default_scalar_type=numpy.float32,
tune_for=op.sym_operator())
from grudge.visualization import SiloVisualizer, VtkVisualizer # noqa
#vis = VtkVisualizer(discr, rcon, "shearflow-%d" % order)
vis = SiloVisualizer(discr, rcon)
fields = box.volume_interpolant(0, discr)
navierstokes_ex = op.bind(discr)
max_eigval = [0]
def rhs(t, q):
ode_rhs, speed = navierstokes_ex(t, q)
max_eigval[0] = speed
return ode_rhs
rhs(0, fields)
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
from grudge.timestep import RK4TimeStepper
stepper = RK4TimeStepper()
# diagnostics setup ---------------------------------------------------
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
logmgr = LogManager("navierstokes-%d.dat" % order, "w", rcon.communicator)
add_run_info(logmgr)
add_general_quantities(logmgr)
add_simulation_quantities(logmgr)
discr.add_instrumentation(logmgr)
stepper.add_instrumentation(logmgr)
logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
from pytools.log import LogQuantity
class ChangeSinceLastStep(LogQuantity):
"""Records the change of a variable between a time step and the previous
one"""
def __init__(self, name="change"):
LogQuantity.__init__(self, name, "1", "Change since last time step")
self.old_fields = 0
def __call__(self):
result = discr.norm(fields - self.old_fields)
self.old_fields = fields
return result
logmgr.add_quantity(ChangeSinceLastStep())
# timestep loop -------------------------------------------------------
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=200,
#max_steps=500,
logmgr=logmgr,
max_dt_getter=lambda t: op.estimate_timestep(discr,
stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
for step, t, dt in step_it:
if step % 200 == 0:
#if False:
visf = vis.make_file("box-%d-%06d" % (order, step))
#rhs_fields = rhs(t, fields)
vis.add_data(visf,
[
("rho", discr.convert_volume(
op.rho(fields), kind="numpy")),
("e", discr.convert_volume(
op.e(fields), kind="numpy")),
("rho_u", discr.convert_volume(
op.rho_u(fields), kind="numpy")),
("u", discr.convert_volume(
op.u(fields), kind="numpy")),
# ("rhs_rho", discr.convert_volume(
# op.rho(rhs_fields), kind="numpy")),
# ("rhs_e", discr.convert_volume(
# op.e(rhs_fields), kind="numpy")),
# ("rhs_rho_u", discr.convert_volume(
# op.rho_u(rhs_fields), kind="numpy")),
],
expressions=[
("p", "(0.4)*(e- 0.5*(rho_u*u))"),
],
time=t, step=step
)
visf.close()
fields = stepper(fields, t, dt, rhs)
finally:
vis.close()
logmgr.save()
discr.close()
if __name__ == "__main__":
main()
examples/gas_dynamics/euler/sine-wave.py deleted 100644 → 0
View file @ 519e3e71
# grudge - the Hybrid'n'Easy DG Environment
# Copyright (C) 2008 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 <http://www.gnu.org/licenses/>.
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
import numpy.linalg as la
class SineWave:
def __init__(self):
self.gamma = 1.4
self.mu = 0
self.prandtl = 0.72
self.spec_gas_const = 287.1
def __call__(self, t, x_vec):
rho = 2 + numpy.sin(x_vec[0] + x_vec[1] + x_vec[2] - 2 * t)
velocity = numpy.array([1, 1, 0])
p = 1
e = p/(self.gamma-1) + rho/2 * numpy.dot(velocity, velocity)
rho_u = rho * velocity[0]
rho_v = rho * velocity[1]
rho_w = rho * velocity[2]
from grudge.tools import join_fields
return join_fields(rho, e, rho_u, rho_v, rho_w)
def properties(self):
return(self.gamma, self.mu, self.prandtl, self.spec_gas_const)
def volume_interpolant(self, t, discr):
return discr.convert_volume(
self(t, discr.nodes.T),
kind=discr.compute_kind)
def boundary_interpolant(self, t, discr, tag):
return discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T),
tag=tag, kind=discr.compute_kind)
def main(final_time=1, write_output=False):
from grudge.backends import guess_run_context
rcon = guess_run_context()
from grudge.tools import EOCRecorder, to_obj_array
eoc_rec = EOCRecorder()
if rcon.is_head_rank:
from grudge.mesh import make_box_mesh
mesh = make_box_mesh((0,0,0), (10,10,10), max_volume=0.5)
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [3, 4, 5]:
discr = rcon.make_discretization(mesh_data, order=order,
default_scalar_type=numpy.float64)
from grudge.visualization import SiloVisualizer, VtkVisualizer
vis = VtkVisualizer(discr, rcon, "sinewave-%d" % order)
#vis = SiloVisualizer(discr, rcon)
sinewave = SineWave()
fields = sinewave.volume_interpolant(0, discr)
gamma, mu, prandtl, spec_gas_const = sinewave.properties()
from grudge.mesh import BTAG_ALL
from grudge.models.gas_dynamics import GasDynamicsOperator
op = GasDynamicsOperator(dimensions=mesh.dimensions, gamma=gamma, mu=mu,
prandtl=prandtl, spec_gas_const=spec_gas_const,
bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave,
inflow_tag=BTAG_ALL, source=None)
euler_ex = op.bind(discr)
max_eigval = [0]
def rhs(t, q):
ode_rhs, speed = euler_ex(t, q)
max_eigval[0] = speed
return ode_rhs
rhs(0, fields)
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
from grudge.timestep import RK4TimeStepper
stepper = RK4TimeStepper()
# diagnostics setup ---------------------------------------------------
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
if write_output:
log_name = ("euler-sinewave-%(order)d-%(els)d.dat"
% {"order":order, "els":len(mesh.elements)})
else:
log_name = False
logmgr = LogManager(log_name, "w", rcon.communicator)
add_run_info(logmgr)
add_general_quantities(logmgr)
add_simulation_quantities(logmgr)
discr.add_instrumentation(logmgr)
stepper.add_instrumentation(logmgr)
logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
# timestep loop -------------------------------------------------------
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, max_eigenvalue=max_eigval[0]))
for step, t, dt in step_it:
#if step % 10 == 0:
if write_output:
visf = vis.make_file("sinewave-%d-%04d" % (order, step))
#from pyvisfile.silo import DB_VARTYPE_VECTOR
vis.add_data(visf,
[
("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
("e", discr.convert_volume(op.e(fields), kind="numpy")),
("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
("u", discr.convert_volume(op.u(fields), kind="numpy")),
#("true_rho", op.rho(true_fields)),
#("true_e", op.e(true_fields)),
#("true_rho_u", op.rho_u(true_fields)),
#("true_u", op.u(true_fields)),
#("rhs_rho", op.rho(rhs_fields)),
#("rhs_e", op.e(rhs_fields)),
#("rhs_rho_u", op.rho_u(rhs_fields)),
],
#expressions=[
#("diff_rho", "rho-true_rho"),
#("diff_e", "e-true_e"),
#("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
#("p", "0.4*(e- 0.5*(rho_u*u))"),
#],
time=t, step=step
)
visf.close()
fields = stepper(fields, t, dt, rhs)
finally:
vis.close()
logmgr.close()
discr.close()
true_fields = sinewave.volume_interpolant(t, discr)
eoc_rec.add_data_point(order, discr.norm(fields-true_fields))
print()
print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
from pytools.test import mark_test
@mark_test.long
def test_euler_sine_wave():
main(final_time=0.1, write_output=False)
examples/gas_dynamics/euler/sod-2d.py deleted 100644 → 0
View file @ 519e3e71
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
import numpy.linalg as la
class Sod:
def __init__(self, gamma):
self.gamma = gamma
self.prandtl = 0.72
def __call__(self, t, x_vec):
from grudge.tools import heaviside
from grudge.tools import heaviside_a
x_rel = x_vec[0]
y_rel = x_vec[1]
from math import pi
r = numpy.sqrt(x_rel**2+y_rel**2)
r_shift=r-3.0
u = 0.0
v = 0.0
from numpy import sign
rho = heaviside(-r_shift)+.125*heaviside_a(r_shift,1.0)
e = (1.0/(self.gamma-1.0))*(heaviside(-r_shift)+.1*heaviside_a(r_shift,1.0))
p = (self.gamma-1.0)*e
from grudge.tools import join_fields
return join_fields(rho, e, rho*u, rho*v)
def volume_interpolant(self, t, discr):
return discr.convert_volume(
self(t, discr.nodes.T),
kind=discr.compute_kind)
def boundary_interpolant(self, t, discr, tag):
return discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T),
tag=tag, kind=discr.compute_kind)
def main():
from grudge.backends import guess_run_context
rcon = guess_run_context()
from grudge.tools import to_obj_array
if rcon.is_head_rank:
from grudge.mesh.generator import make_rect_mesh
mesh = make_rect_mesh((-5,-5), (5,5), max_area=0.01)
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [1]:
discr = rcon.make_discretization(mesh_data, order=order,
default_scalar_type=numpy.float64)
from grudge.visualization import SiloVisualizer, VtkVisualizer
vis = VtkVisualizer(discr, rcon, "Sod2D-%d" % order)
#vis = SiloVisualizer(discr, rcon)
sod_field = Sod(gamma=1.4)
fields = sod_field.volume_interpolant(0, discr)
from grudge.models.gas_dynamics import GasDynamicsOperator
from grudge.mesh import BTAG_ALL
op = GasDynamicsOperator(dimensions=2, gamma=sod_field.gamma, mu=0.0,
prandtl=sod_field.prandtl,
bc_inflow=sod_field,
bc_outflow=sod_field,
bc_noslip=sod_field,
inflow_tag=BTAG_ALL,
source=None)
euler_ex = op.bind(discr)
max_eigval = [0]
def rhs(t, q):
ode_rhs, speed = euler_ex(t, q)
max_eigval[0] = speed
return ode_rhs
rhs(0, fields)
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
# limiter setup ------------------------------------------------------------
from grudge.models.gas_dynamics import SlopeLimiter1NEuler
limiter = SlopeLimiter1NEuler(discr, sod_field.gamma, 2, op)
# integrator setup---------------------------------------------------------
from grudge.timestep import SSPRK3TimeStepper, RK4TimeStepper
stepper = SSPRK3TimeStepper(limiter=limiter)
#stepper = SSPRK3TimeStepper()
#stepper = RK4TimeStepper()
# diagnostics setup ---------------------------------------------------
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
logmgr = LogManager("euler-%d.dat" % order, "w", rcon.communicator)
add_run_info(logmgr)
add_general_quantities(logmgr)
add_simulation_quantities(logmgr)
discr.add_instrumentation(logmgr)
stepper.add_instrumentation(logmgr)
logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
# filter setup-------------------------------------------------------------
from grudge.discretization import Filter, ExponentialFilterResponseFunction
mode_filter = Filter(discr,
ExponentialFilterResponseFunction(min_amplification=0.9,order=4))
# timestep loop -------------------------------------------------------
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=1.0, logmgr=logmgr,
max_dt_getter=lambda t: op.estimate_timestep(discr,
stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
for step, t, dt in step_it:
if step % 5 == 0:
#if False:
visf = vis.make_file("vortex-%d-%04d" % (order, step))
#true_fields = vortex.volume_interpolant(t, discr)
#from pyvisfile.silo import DB_VARTYPE_VECTOR
vis.add_data(visf,
[
("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
("e", discr.convert_volume(op.e(fields), kind="numpy")),
("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
("u", discr.convert_volume(op.u(fields), kind="numpy")),
#("true_rho", op.rho(true_fields)),
#("true_e", op.e(true_fields)),
#("true_rho_u", op.rho_u(true_fields)),
#("true_u", op.u(true_fields)),
#("rhs_rho", op.rho(rhs_fields)),
#("rhs_e", op.e(rhs_fields)),
#("rhs_rho_u", op.rho_u(rhs_fields)),
],
#expressions=[
#("diff_rho", "rho-true_rho"),
#("diff_e", "e-true_e"),
#("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
#("p", "0.4*(e- 0.5*(rho_u*u))"),
#],
time=t, step=step
)
visf.close()
fields = stepper(fields, t, dt, rhs)
# fields = limiter(fields)
# fields = mode_filter(fields)
assert not numpy.isnan(numpy.sum(fields[0]))
finally:
vis.close()
logmgr.close()
discr.close()
# not solution, just to check against when making code changes
true_fields = sod_field.volume_interpolant(t, discr)
print(discr.norm(fields-true_fields))
if __name__ == "__main__":
main()