# 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/>.
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import numpy
def main(write_output=True, flux_type_arg="central", use_quadrature=True,
final_time=20):
from math import sin, cos, pi, sqrt
from grudge.backends import guess_run_context
rcon = guess_run_context()
# mesh setup --------------------------------------------------------------
if rcon.is_head_rank:
#from grudge.mesh.generator import make_disk_mesh
#mesh = make_disk_mesh()
from grudge.mesh.generator import make_rect_mesh
mesh = make_rect_mesh(a=(-1,-1),b=(1,1),max_area=0.008)
if rcon.is_head_rank:
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
# space-time-dependent-velocity-field -------------------------------------
# simple vortex
class TimeDependentVField:
""" `TimeDependentVField` is a callable expecting `(x, t)` representing space and time
`x` is of the length of the spatial dimension and `t` is the time."""
shape = (2,)
def __call__(self, pt, el, t):
x, y = pt
# Correction-Factor to make the speed zero on the on the boundary
#fac = (1-x**2)*(1-y**2)
fac = 1.
return numpy.array([-y*fac, x*fac]) * cos(pi*t)
class VField:
""" `VField` is a callable expecting `(x)` representing space
`x` is of the length of the spatial dimension."""
shape = (2,)
def __call__(self, pt, el):
x, y = pt
# Correction-Factor to make the speed zero on the on the boundary
#fac = (1-x**2)*(1-y**2)
fac = 1.
return numpy.array([-y*fac, x*fac])
# space-time-dependent State BC (optional)-----------------------------------
class TimeDependentBc_u:
""" space and time dependent BC for state u"""
def __call__(self, pt, el, t):
x, y = pt
if t <= 0.5:
if x > 0:
return 1
else:
return 0
else:
return 0
class Bc_u:
""" Only space dependent BC for state u"""
def __call__(seld, pt, el):
x, y = pt
if x > 0:
return 1
else:
return 0
# operator setup ----------------------------------------------------------
# In the operator setup it is possible to switch between a only space
# dependent velocity field `VField` or a time and space dependent
# `TimeDependentVField`.
# For `TimeDependentVField`: advec_v=TimeDependentGivenFunction(VField())
# For `VField`: advec_v=TimeConstantGivenFunction(GivenFunction(VField()))
# Same for the Bc_u Function! If you don't define Bc_u then the BC for u = 0.
from grudge.data import \
ConstantGivenFunction, \
TimeConstantGivenFunction, \
TimeDependentGivenFunction, \
GivenFunction
from grudge.models.advection import VariableCoefficientAdvectionOperator
op = VariableCoefficientAdvectionOperator(mesh.dimensions,
#advec_v=TimeDependentGivenFunction(
# TimeDependentVField()),
advec_v=TimeConstantGivenFunction(
GivenFunction(VField())),
#bc_u_f=TimeDependentGivenFunction(
# TimeDependentBc_u()),
bc_u_f=TimeConstantGivenFunction(
GivenFunction(Bc_u())),
flux_type=flux_type_arg)
# discretization setup ----------------------------------------------------
order = 5
if use_quadrature:
quad_min_degrees = {"quad": 3*order}
else:
quad_min_degrees = {}
discr = rcon.make_discretization(mesh_data, order=order,
default_scalar_type=numpy.float64,
debug=["cuda_no_plan"],
quad_min_degrees=quad_min_degrees,
tune_for=op.sym_operator(),
)
vis_discr = discr
# visualization setup -----------------------------------------------------
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(vis_discr, rcon, "fld")
# initial condition -------------------------------------------------------
if True:
def initial(pt, el):
# Gauss pulse
from math import exp
x = (pt-numpy.array([0.3, 0.5]))*8
return exp(-numpy.dot(x, x))
else:
def initial(pt, el):
# Rectangle
x, y = pt
if abs(x) < 0.5 and abs(y) < 0.2:
return 2
else:
return 1
u = discr.interpolate_volume_function(initial)
# timestep setup ----------------------------------------------------------
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper(
vector_primitive_factory=discr.get_vector_primitive_factory())
if rcon.is_head_rank:
print("%d elements" % len(discr.mesh.elements))
# filter setup-------------------------------------------------------------
from grudge.discretization import ExponentialFilterResponseFunction
from grudge.symbolic.operators import FilterOperator
mode_filter = FilterOperator(
ExponentialFilterResponseFunction(min_amplification=0.9,order=4))\
.bind(discr)
# diagnostics setup -------------------------------------------------------
from pytools.log import LogManager, \
add_general_quantities, \
add_simulation_quantities, \
add_run_info
if write_output:
log_file_name = "space-dep.dat"
else:
log_file_name = None
logmgr = LogManager(log_file_name, "w", rcon.communicator)
add_run_info(logmgr)
add_general_quantities(logmgr)
add_simulation_quantities(logmgr)
discr.add_instrumentation(logmgr)
stepper.add_instrumentation(logmgr)
from grudge.log import Integral, LpNorm
u_getter = lambda: u
logmgr.add_quantity(Integral(u_getter, discr, name="int_u"))
logmgr.add_quantity(LpNorm(u_getter, discr, p=1, name="l1_u"))
logmgr.add_quantity(LpNorm(u_getter, discr, name="l2_u"))
logmgr.add_watches(["step.max", "t_sim.max", "l2_u", "t_step.max"])
# Initialize v for data output:
v = op.advec_v.volume_interpolant(0, discr)
# timestep loop -----------------------------------------------------------
rhs = op.bind(discr)
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=final_time, logmgr=logmgr,
max_dt_getter=lambda t: op.estimate_timestep(discr,
stepper=stepper, t=t, fields=u))
for step, t, dt in step_it:
if step % 10 == 0 and write_output:
visf = vis.make_file("fld-%04d" % step)
vis.add_data(visf, [
("u", discr.convert_volume(u, kind="numpy")),
("v", discr.convert_volume(v, kind="numpy"))
], time=t, step=step)
visf.close()
u = stepper(u, t, dt, rhs)
# We're feeding in a discontinuity through the BCs.
# Quadrature does not help with shock capturing--
# therefore we do need to filter here, regardless
# of whether quadrature is enabled.
u = mode_filter(u)
assert discr.norm(u) < 10
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
def test_var_velocity_advection():
from pytools.test import mark_test
mark_long = mark_test.long
for flux_type in ["upwind", "central", "lf"]:
for use_quadrature in [False, True]:
descr = "variable-velocity-advection with %s flux" % flux_type
if use_quadrature:
descr += " and quadrature"
yield descr, mark_long(main), False, flux_type, use_quadrature, 1