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
Commits on Source (26)
Hide whitespace changes
Inline Side-by-side
Showing
with 801 additions and 1014 deletions
examples/maxwell/analytic_solutions.py
View file @ c2196857
# 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
from grudge.tools import \
cyl_bessel_j, \
cyl_bessel_j_prime
from math import sqrt, pi, sin, cos, atan2, acos
import numpy
from __future__ import division, absolute_import, print_function
__copyright__ = """
Copyright (C) 2007-17 Andreas Kloeckner
Copyright (C) 2017 Bogdan Enache
"""
__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.
"""
from six.moves import range, zip
import numpy as np
import numpy.linalg as la
import cmath
from six.moves import range
from six.moves import zip
from grudge import sym
# solution adapters -----------------------------------------------------------
class RealPartAdapter:
def __init__(self, adaptee):
self.adaptee = adaptee
@property
def shape(self):
return self.adaptee.shape
def __call__(self, x, el):
return [xi.real for xi in self.adaptee(x, el)]
class SplitComplexAdapter:
def __init__(self, adaptee):
self.adaptee = adaptee
@property
def shape(self):
(n,) = self.adaptee.shape
return (n*2,)
def __call__(self, x, el):
ad_x = self.adaptee(x, el)
return [xi.real for xi in ad_x] + [xi.imag for xi in ad_x]
# {{{ cylindrical field adapter
class CylindricalFieldAdapter:
"""Turns a field returned as (r, phi, z) components into
......@@ -90,8 +66,10 @@ class CylindricalFieldAdapter:
return result
# }}}
# {{{ spherical field adapter
class SphericalFieldAdapter:
"""Turns the adaptee's field C{(Er, Etheta, Ephi)} components into
......@@ -138,10 +116,11 @@ class SphericalFieldAdapter:
return result
# }}}
# {{{ cylindrical cavity mode
# actual solutions ------------------------------------------------------------
class CylindricalCavityMode:
"""A cylindrical TM cavity mode.
......@@ -150,7 +129,7 @@ class CylindricalCavityMode:
3rd edition, 2001.
ch 8.7, p. 368f.
"""
def __init__(self, m, n, p, radius, height, epsilon, mu):
try:
from bessel_zeros import bessel_zeros
......@@ -211,7 +190,7 @@ class CylindricalCavityMode:
# psi and derivatives ---------------------------------------------
psi = cyl_bessel_j(m, gamma * r) * phi_factor
psi_dr = gamma*cyl_bessel_j_prime(m, gamma*r) * phi_factor
psi_dphi = (cyl_bessel_j(m, gamma * r)
psi_dphi = (cyl_bessel_j(m, gamma * r)
* 1/r * phi_sign*1j*m * phi_factor)
# field components in polar coordinates ---------------------------
......@@ -237,91 +216,67 @@ class CylindricalCavityMode:
return [er, ephi, ez, hr, hphi, hz]
# }}}
# {{{ symmetric rectangular waveguide mode
class RectangularWaveguideMode:
"""A rectangular TM cavity mode."""
def __init__(self, epsilon, mu, mode_indices,
dimensions=(1,1,1), coefficients=(1,0,0),
forward_coeff=1, backward_coeff=0):
for n in mode_indices:
assert n >= 0 and n == int(n)
self.mode_indices = mode_indices
self.dimensions = dimensions
self.coefficients = coefficients
self.forward_coeff = forward_coeff
self.backward_coeff = backward_coeff
def get_rectangular_waveguide_mode(
ambient_dim, mode_indices,
dimensions=(1, 1, 1), coefficients=(1, 0, 0),
forward_coeff=1, backward_coeff=0,
):
"""A TM cavity mode.
self.epsilon = epsilon
self.mu = mu
self.t = 0
self.factors = [n*pi/a for n, a in zip(self.mode_indices, self.dimensions)]
c = 1/sqrt(mu*epsilon)
self.k = sqrt(sum(f**2 for f in self.factors))
self.omega = self.k*c
def set_time(self, t):
self.t = t
def __call__(self, discr):
f,g,h = self.factors
omega = self.omega
k = self.k
Returns an expression depending on *epsilon*, *mu*, and *t*.
"""
x = discr.nodes[:,0]
y = discr.nodes[:,1]
if discr.dimensions > 2:
z = discr.nodes[:,2]
else:
z = discr.volume_zeros()
f, g, h = factors = [n*np.pi/a for n, a in zip(mode_indices, dimensions)]
k = np.sqrt(sum(f**2 for f in factors))
sx = numpy.sin(f*x)
sy = numpy.sin(g*y)
cx = numpy.cos(f*x)
cy = numpy.cos(g*y)
epsilon = sym.ScalarVariable("epsilon")
mu = sym.ScalarVariable("mu")
zdep_add = numpy.exp(1j*h*z)+numpy.exp(-1j*h*z)
zdep_sub = numpy.exp(1j*h*z)-numpy.exp(-1j*h*z)
c = sym.cse(1/sym.sqrt(mu*epsilon), "c")
omega = k*c
tdep = numpy.exp(-1j * omega * self.t)
nodes = sym.nodes(ambient_dim)
x = nodes[0]
y = nodes[1]
if nodes.shape[0] > 2:
z = nodes[2]
else:
z = 0
C = 1j/(f**2+g**2)
sx = sym.sin(f*x)
sy = sym.sin(g*y)
cx = sym.cos(f*x)
cy = sym.cos(g*y)
result = numpy.zeros((6,len(x)), dtype=zdep_add.dtype)
zdep_add = sym.exp(1j*h*z)+sym.exp(-1j*h*z)
zdep_sub = sym.exp(1j*h*z)-sym.exp(-1j*h*z)
result[0] = C*f*h*cx*sy*zdep_sub*tdep # ex
result[1] = C*g*h*sx*cy*zdep_sub*tdep # ey
result[2] = sx*sy*zdep_add*tdep # ez
result[3] = -C*g*self.epsilon*omega*sx*cy*zdep_add*tdep # hx
result[4] = C*f*self.epsilon*omega*cx*sy*zdep_add*tdep
tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
return result
const = 1j/(f**2+g**2)
result = sym.join_fields(
const*f*h*cx*sy*zdep_sub*tdep, # ex
const*g*h*sx*cy*zdep_sub*tdep, # ey
sx*sy*zdep_add*tdep, # ez
-const*g*epsilon*omega*sx*cy*zdep_add*tdep, # hx
const*f*epsilon*omega*cx*sy*zdep_add*tdep,
const*f*epsilon*omega*cx*sy*zdep_add*tdep,
) # THIS IS BAD, last line should be all zeros
return result
class RectangularCavityMode(RectangularWaveguideMode):
"""A rectangular TM cavity mode."""
def __init__(self, *args, **kwargs):
if "scale" in kwargs:
kwargs["forward_coeff"] = scale
kwargs["backward_coeff"] = scale
else:
kwargs["forward_coeff"] = 1
kwargs["backward_coeff"] = 1
RectangularWaveguideMode.__init__(self, *args, **kwargs)
# }}}
# {{{ dipole solution
# dipole solution -------------------------------------------------------------
class DipoleFarField:
"""Dipole EM field.
......@@ -384,16 +339,14 @@ class DipoleFarField:
j0 = self.q*self.d*self.omega/self.epsilon
return j0*cos(self.omega*t)
# }}}
# {{{ check_time_harmonic_solution
# analytic solution tools -----------------------------------------------------
def check_time_harmonic_solution(discr, mode, c_sol):
from grudge.discretization import bind_nabla, bind_mass_matrix
from grudge.visualization import SiloVisualizer
from grudge.silo import SiloFile
from grudge.tools import dot, cross
from grudge.silo import DB_VARTYPE_VECTOR
def curl(field):
return cross(nabla, field)
......@@ -432,15 +385,15 @@ def check_time_harmonic_solution(discr, mode, c_sol):
silo = SiloFile("em-complex-%04d.silo" % step)
vis.add_to_silo(silo,
vectors=[
("curl_er", curl(er)),
("om_hi", -mode.mu*mode.omega*hi),
("curl_hr", curl(hr)),
("om_ei", mode.epsilon*mode.omega*hi),
("curl_er", curl(er)),
("om_hi", -mode.mu*mode.omega*hi),
("curl_hr", curl(hr)),
("om_ei", mode.epsilon*mode.omega*hi),
],
expressions=[
("diff_er", "curl_er-om_hi", DB_VARTYPE_VECTOR),
("diff_hr", "curl_hr-om_ei", DB_VARTYPE_VECTOR),
],
("diff_er", "curl_er-om_hi", DB_VARTYPE_VECTOR),
("diff_hr", "curl_hr-om_ei", DB_VARTYPE_VECTOR),
],
write_coarse_mesh=True,
time=t, step=step
)
......@@ -458,3 +411,6 @@ def check_time_harmonic_solution(discr, mode, c_sol):
rel_l2_error(hi_res, hi),
))
# }}}
# vim: foldmethod=marker
examples/maxwell/cavities.py
View file @ c2196857
# 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 as np
import logging
logger = logging.getLogger(__name__)
def main(write_output=True, allow_features=None, flux_type_arg=1,
bdry_flux_type_arg=None, extra_discr_args={}):
from grudge.mesh.generator import make_cylinder_mesh, make_box_mesh
from grudge.tools import EOCRecorder, to_obj_array
from math import sqrt, pi # noqa
from analytic_solutions import ( # noqa
RealPartAdapter,
SplitComplexAdapter,
CylindricalFieldAdapter,
CylindricalCavityMode,
RectangularWaveguideMode,
RectangularCavityMode)
from grudge.models.em import MaxwellOperator
"""Minimal example of a grudge driver."""
logging.basicConfig(level=logging.DEBUG)
from __future__ import division, print_function
from grudge.backends import guess_run_context
rcon = guess_run_context(allow_features)
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
mu0 = 4*pi*1e-7 # N/A**2.
epsilon = 1*epsilon0
mu = 1*mu0
__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:
eoc_rec = EOCRecorder()
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
cylindrical = False
periodic = False
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.
"""
if cylindrical:
R = 1
d = 2
mode = CylindricalCavityMode(m=1, n=1, p=1,
radius=R, height=d,
epsilon=epsilon, mu=mu)
# r_sol = CylindricalFieldAdapter(RealPartAdapter(mode))
# c_sol = SplitComplexAdapter(CylindricalFieldAdapter(mode))
if rcon.is_head_rank:
mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
else:
if periodic:
mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1))
periodicity = (False, False, True)
else:
periodicity = None
mode = RectangularCavityMode(epsilon, mu, (1, 2, 2))
if rcon.is_head_rank:
mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
if rcon.is_head_rank:
mesh_data = rcon.distribute_mesh(mesh)
import numpy as np
import pyopencl as cl
from grudge.shortcuts import set_up_rk4, set_up_forward_euler # noqa: F401
from grudge import sym, bind, Discretization
from analytic_solutions import ( # noqa
get_rectangular_waveguide_mode,
#RectangularCavityMode,
#RectangularWaveguideMode,
#RectangularCavityMode
)
def main(write_output=True, order=4):
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
dims = 2
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dims,
b=(0.5,)*dims,
n=(20,)*dims)
discr = Discretization(cl_ctx, mesh, order=order)
if 0:
epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
mu0 = 4*np.pi*1e-7 # N/A**2.
epsilon = 1*epsilon0
mu = 1*mu0
else:
mesh_data = rcon.receive_mesh()
for order in [4, 5, 6]:
#for order in [1,2,3,4,5,6]:
extra_discr_args.setdefault("debug", []).extend([
"cuda_no_plan", "cuda_dump_kernels"])
op = MaxwellOperator(epsilon, mu,
flux_type=flux_type_arg,
bdry_flux_type=bdry_flux_type_arg)
discr = rcon.make_discretization(mesh_data, order=order,
tune_for=op.sym_operator(),
**extra_discr_args)
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "em-%d" % order)
mode.set_time(0)
def get_true_field():
return discr.convert_volume(
to_obj_array(mode(discr)
.real.astype(discr.default_scalar_type).copy()),
kind=discr.compute_kind)
fields = get_true_field()
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper(dtype=discr.default_scalar_type, rcon=rcon)
#from grudge.timestep.dumka3 import Dumka3TimeStepper
#stepper = Dumka3TimeStepper(3, dtype=discr.default_scalar_type, rcon=rcon)
# {{{ diagnostics setup
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
if write_output:
log_file_name = "maxwell-%d.dat" % order
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 pytools.log import IntervalTimer
vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
logmgr.add_quantity(vis_timer)
epsilon = 1
mu = 1
from grudge.log import EMFieldGetter, add_em_quantities
field_getter = EMFieldGetter(discr, op, lambda: fields)
add_em_quantities(logmgr, op, field_getter)
logmgr.add_watches(
["step.max", "t_sim.max",
("W_field", "W_el+W_mag"),
"t_step.max"]
)
# }}}
# {{{ timestep loop
rhs = op.bind(discr)
final_time = 0.5e-9
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=fields))
for step, t, dt in step_it:
if step % 50 == 0 and write_output:
sub_timer = vis_timer.start_sub_timer()
e, h = op.split_eh(fields)
visf = vis.make_file("em-%d-%04d" % (order, step))
vis.add_data(
visf,
[
("e",
discr.convert_volume(e, kind="numpy")),
("h",
discr.convert_volume(h, kind="numpy")),
],
time=t, step=step)
visf.close()
sub_timer.stop().submit()
fields = stepper(fields, t, dt, rhs)
mode.set_time(final_time)
eoc_rec.add_data_point(order,
discr.norm(fields-get_true_field()))
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
from grudge.models.em import MaxwellOperator
op = MaxwellOperator(epsilon, mu, flux_type="lf", dimensions=dims)
if rcon.is_head_rank:
print()
print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
queue = cl.CommandQueue(discr.cl_context)
# }}}
sym_mode = get_rectangular_waveguide_mode(mesh.ambient_dim, (1, 2, 2))
fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
#fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr)
assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
#fields = (discr.volume_discr)
#fields = join_fields(discr.zeros(queue),
#[discr.zeros(queue) for i in range(discr.dim*2 - 1)])
# FIXME
#dt = op.estimate_rk4_timestep(discr, fields=fields)
# {{{ entry points for py.test
op.check_bc_coverage(mesh)
from pytools.test import mark_test
# print(sym.pretty(op.sym_operator()))
bound_op = bind(discr, op.sym_operator())
print(bound_op)
# 1/0
def rhs(t, w):
return bound_op(queue, t=t, w=w)
@mark_test.long
def test_maxwell_cavities():
main(write_output=False)
if mesh.dim == 2:
dt = 0.0004
elif mesh.dim == 3:
dt = 0.0002
dt_stepper = set_up_forward_euler("w", dt, fields, rhs)
@mark_test.long
def test_maxwell_cavities_lf():
main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1)
final_t = dt * 300
nsteps = int(final_t/dt)
print("dt=%g nsteps=%d" % (dt, nsteps))
@mark_test.mpi
@mark_test.long
def test_maxwell_cavities_mpi():
from pytools.mpi import run_with_mpi_ranks
run_with_mpi_ranks(__file__, 2, main,
write_output=False, allow_features=["mpi"])
from grudge.shortcuts import make_visualizer
vis = make_visualizer(discr, vis_order=order)
step = 0
def test_cuda():
try:
from pycuda.tools import mark_cuda_test
except ImportError:
pass
norm = bind(discr, sym.norm(2, sym.var("u")))
yield "SP CUDA Maxwell", mark_cuda_test(
do_test_maxwell_cavities_cuda), np.float32
yield "DP CUDA Maxwell", mark_cuda_test(
do_test_maxwell_cavities_cuda), np.float64
from time import time
t_last_step = time()
e, h = op.split_eh(fields)
def do_test_maxwell_cavities_cuda(dtype):
import pytest # noqa
if 1:
vis.write_vtk_file("fld-%04d.vtu" % step,
[
("e", e),
("h", h),
])
main(write_output=False, allow_features=["cuda"],
extra_discr_args=dict(
debug=["cuda_no_plan"],
default_scalar_type=dtype,
))
for event in dt_stepper.run(t_end=final_t):
if isinstance(event, dt_stepper.StateComputed):
assert event.component_id == "w"
# }}}
step += 1
print(step, event.t, norm(queue, u=event.state_component[0]),
time()-t_last_step)
if step % 10 == 0:
e, h = op.split_eh(event.state_component)
vis.write_vtk_file("fld-%04d.vtu" % step,
[
("e", e),
("h", h),
])
t_last_step = time()
# entry point -----------------------------------------------------------------
if __name__ == "__main__":
from pytools.mpi import in_mpi_relaunch
if in_mpi_relaunch():
test_maxwell_cavities_mpi()
else:
do_test_maxwell_cavities_cuda(np.float32)
main()
examples/maxwell/testing.py 0 → 100644
View file @ c2196857
"""Minimal example of a grudge driver."""
from __future__ import division, print_function
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import numpy as np
import pyopencl as cl
from grudge.shortcuts import set_up_rk4
from grudge import sym, bind, Discretization
from analytic_solutions import ( # noqa
RectangularWaveguideMode,
RectangularCavityMode,
SymRectangularWaveguideMode,
SymRectangularCavityMode)
def main(write_output=True, order=4):
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
from grudge.models.em import MaxwellOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
queue = cl.CommandQueue(cl_ctx)
from pytools.obj_array import join_fields
from pymbolic.primitives import Variable as var
import pymbolic.primitives as p
import loopy as lp
nom = "abs"
i = var("i")
func = var(nom)
if("fabs" == nom):
func = var("abs")
def knl():
knl = lp.make_kernel(
"{[i]: 0<=i<n}",
[
lp.Assignment(var("out")[i],
func(var("a")[i]))
])
return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
numpyA = np.array([-0.5j,-0.1,0,0.1j,0.2 + 3j,0.3,0.4,-6+3j])
a = cl.array.Array(queue,numpyA.shape,dtype= numpyA.dtype)
a.set(numpyA)
evt, (out,) = knl()(queue, a=a)
print(out)
if __name__ == "__main__":
main()
examples/wave/var2.py 0 → 100644
View file @ c2196857
"""Minimal example of a grudge driver."""
from __future__ import division, print_function
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import numpy as np
import pyopencl as cl
from grudge.shortcuts import set_up_rk4
from grudge import sym, bind, Discretization
def main(write_output=True, order=4):
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
dims = 2
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dims,
b=(0.5,)*dims,
n=(8,)*dims)
discr = Discretization(cl_ctx, mesh, order=order)
source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
source_width = 0.05
source_omega = 3
sym_x = sym.nodes(mesh.dim)
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
c = sym.If(sym.Comparison(
np.dot(sym_x, sym_x), "<", 0.2),
-0.1, -0.2)
#c = np.dot(sym_x,sym_x)
from grudge.models.wave import VariableCoefficientWeakWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
op = VariableCoefficientWeakWaveOperator(c,
discr.dim,
source_f=(
sym.sin(source_omega*sym_t)
* sym.exp(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2)),
dirichlet_tag=BTAG_NONE,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_ALL,
flux_type="upwind")
queue = cl.CommandQueue(discr.cl_context)
from pytools.obj_array import join_fields
fields = join_fields(discr.zeros(queue),
[discr.zeros(queue) for i in range(discr.dim)])
# FIXME
#dt = op.estimate_rk4_timestep(discr, fields=fields)
op.check_bc_coverage(mesh)
# print(sym.pretty(op.sym_operator()))
bound_op = bind(discr, op.sym_operator())
print(bound_op)
# 1/0
def rhs(t, w):
return bound_op(queue, t=t, w=w)
if mesh.dim == 2:
dt = 0.04
elif mesh.dim == 3:
dt = 0.02
dt_stepper = set_up_rk4("w", dt, fields, rhs)
final_t = 10
nsteps = int(final_t/dt)
print("dt=%g nsteps=%d" % (dt, nsteps))
from grudge.shortcuts import make_visualizer
vis = make_visualizer(discr, vis_order=order)
step = 0
norm = bind(discr, sym.norm(2, sym.var("u")))
from time import time
t_last_step = time()
for event in dt_stepper.run(t_end=final_t):
if isinstance(event, dt_stepper.StateComputed):
assert event.component_id == "w"
step += 1
print(step, event.t, norm(queue, u=event.state_component[0]),
time()-t_last_step)
if step % 10 == 0:
vis.write_vtk_file("fld-%04d.vtu" % step,
[
("u", event.state_component[0]),
("v", event.state_component[1:]),
])
t_last_step = time()
if __name__ == "__main__":
main()
examples/wave/weak.py 0 → 100644
View file @ c2196857
"""Minimal example of a grudge driver."""
from __future__ import division, print_function
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import numpy as np
import pyopencl as cl
from grudge.shortcuts import set_up_rk4
from grudge import sym, bind, Discretization
def main(write_output=True, order=4):
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
dims = 2
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dims,
b=(0.5,)*dims,
n=(8,)*dims)
discr = Discretization(cl_ctx, mesh, order=order)
source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
source_width = 0.05
source_omega = 3
sym_x = sym.nodes(mesh.dim)
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
from grudge.models.wave import WeakWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
op = WeakWaveOperator(-0.1, discr.dim,
source_f=(
sym.sin(source_omega*sym_t)
* sym.exp(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2)),
dirichlet_tag=BTAG_NONE,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_ALL,
flux_type="upwind")
queue = cl.CommandQueue(discr.cl_context)
from pytools.obj_array import join_fields
fields = join_fields(discr.zeros(queue),
[discr.zeros(queue) for i in range(discr.dim)])
# FIXME
#dt = op.estimate_rk4_timestep(discr, fields=fields)
op.check_bc_coverage(mesh)
# print(sym.pretty(op.sym_operator()))
bound_op = bind(discr, op.sym_operator())
print(bound_op)
# 1/0
def rhs(t, w):
return bound_op(queue, t=t, w=w)
if mesh.dim == 2:
dt = 0.04
elif mesh.dim == 3:
dt = 0.02
dt_stepper = set_up_rk4("w", dt, fields, rhs)
final_t = 10
nsteps = int(final_t/dt)
print("dt=%g nsteps=%d" % (dt, nsteps))
from grudge.shortcuts import make_visualizer
vis = make_visualizer(discr, vis_order=order)
step = 0
norm = bind(discr, sym.norm(2, sym.var("u")))
from time import time
t_last_step = time()
for event in dt_stepper.run(t_end=final_t):
if isinstance(event, dt_stepper.StateComputed):
assert event.component_id == "w"
step += 1
print(step, event.t, norm(queue, u=event.state_component[0]),
time()-t_last_step)
if step % 10 == 0:
vis.write_vtk_file("fld-%04d.vtu" % step,
[
("u", event.state_component[0]),
("v", event.state_component[1:]),
])
t_last_step = time()
if __name__ == "__main__":
main()
grudge/execution.py
View file @ c2196857
......@@ -90,14 +90,40 @@ class ExecutionMapper(mappers.Evaluator,
# FIXME: Make a way to register functions
pars = [self.rec(p) for p in expr.parameters]
if any(isinstance(par, cl.array.Array) for par in pars):
import pyopencl.clmath as clmath
func = getattr(clmath, expr.function.name)
else:
args = [self.rec(p) for p in expr.parameters]
from numbers import Number
representative_arg = args[0]
if (
isinstance(representative_arg, Number)
or (isinstance(representative_arg, np.ndarray)
and representative_arg.shape == ())):
func = getattr(np, expr.function.name)
return func(*args)
cached_name = "map_call_knl_"
i = Variable("i")
func = Variable(expr.function.name)
if expr.function.name == "fabs": # FIXME
func = Variable("abs")
cached_name += "abs"
else:
cached_name += expr.function.name
@memoize_in(self.bound_op, cached_name)
def knl():
knl = lp.make_kernel(
"{[i]: 0<=i<n}",
[
lp.Assignment(Variable("out")[i],
func(Variable("a")[i]))
])
return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
return func(*pars)
assert len(args) == 1
evt, (out,) = knl()(self.queue, a=args[0])
return out
def map_nodal_sum(self, op, field_expr):
# FIXME: Could allow array scalars
......@@ -330,6 +356,23 @@ class ExecutionMapper(mappers.Evaluator,
# }}}
# {{{ code execution functions
# DEBUG stuff
def get_max(self, arr):
from pymbolic.primitives import Variable
i = Variable("i")
def knl2():
knl = lp.make_kernel(
"{[i]: 0<=i<n}",
[
lp.Assignment(Variable("out")[i],
Variable("abs")(Variable("a")[i]))
])
return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
evt, (out,) = knl2()(self.queue, a=arr)
return pyopencl.array.max(out)
def map_insn_loopy_kernel(self, insn):
kwargs = {}
......@@ -485,27 +528,41 @@ class BoundOperator(object):
# {{{ process_sym_operator function
def process_sym_operator(sym_operator, post_bind_mapper=None,
dumper=lambda name, sym_operator: None, mesh=None):
def process_sym_operator(discr, sym_operator, post_bind_mapper=None,
dumper=lambda name, sym_operator: None):
import grudge.symbolic.mappers as mappers
dumper("before-bind", sym_operator)
sym_operator = mappers.OperatorBinder()(sym_operator)
mappers.ErrorChecker(mesh)(sym_operator)
mappers.ErrorChecker(discr.mesh)(sym_operator)
if post_bind_mapper is not None:
dumper("before-postbind", sym_operator)
sym_operator = post_bind_mapper(sym_operator)
if mesh is not None:
dumper("before-empty-flux-killer", sym_operator)
sym_operator = mappers.EmptyFluxKiller(mesh)(sym_operator)
dumper("before-empty-flux-killer", sym_operator)
sym_operator = mappers.EmptyFluxKiller(discr.volume_discr.mesh)(sym_operator)
dumper("before-cfold", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
# Work around https://github.com/numpy/numpy/issues/9438
#
# The idea is that we need 1j as an expression to survive
# until code generation time. If it is evaluated and combined
# with other constants, we will need to determine its size
# (as np.complex64/128) within the expression. But because
# of the above numpy bug, sized numbers are not likely to survive
# expression building--so that's why we step in here to fix that.
dumper("before-csize", sym_operator)
sym_operator = mappers.ConstantToNumpyConversionMapper(
real_type=discr.volume_discr.real_dtype.type,
complex_type=discr.volume_discr.complex_dtype.type,
)(sym_operator)
# Ordering restriction:
#
# - Must run constant fold before first type inference pass, because zeros,
......@@ -526,9 +583,8 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
# because otherwise it won't differentiate the type of grids (node or quadrature
# grids) that the operators will apply on.
assert mesh is not None
dumper("before-global-to-reference", sym_operator)
sym_operator = mappers.GlobalToReferenceMapper(mesh.ambient_dim)(sym_operator)
sym_operator = mappers.GlobalToReferenceMapper(discr.ambient_dim)(sym_operator)
# Ordering restriction:
#
......@@ -569,10 +625,11 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
pretty(sym_operator))
stage[0] += 1
sym_operator = process_sym_operator(sym_operator,
sym_operator = process_sym_operator(
discr,
sym_operator,
post_bind_mapper=post_bind_mapper,
dumper=dump_optemplate,
mesh=discr.mesh)
dumper=dump_optemplate)
from grudge.symbolic.compiler import OperatorCompiler
discr_code, eval_code = OperatorCompiler(discr)(sym_operator)
......
grudge/models/em.py
View file @ c2196857
# -*- coding: utf8 -*-
"""grudge operators modelling electromagnetic phenomena."""
from __future__ import division
from __future__ import absolute_import
from six.moves import range
from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2007, 2010 Andreas Kloeckner, David Powell"
__copyright__ = """
Copyright (C) 2007-2017 Andreas Kloeckner
Copyright (C) 2010 David Powell
Copyright (C) 2017 Bogdan Enache
"""
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
......@@ -27,11 +29,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
from six.moves import range
from pytools import memoize_method
from grudge.models import HyperbolicOperator
from grudge.symbolic.primitives import make_common_subexpression as cse
from meshmode.mesh import BTAG_ALL, BTAG_NONE
from grudge import sym
from pytools.obj_array import join_fields, make_obj_array
# TODO: Check PML
......@@ -106,17 +110,8 @@ class MaxwellOperator(HyperbolicOperator):
self.current = current
self.incident_bc_data = incident_bc
@property
def c(self):
from warnings import warn
warn("MaxwellOperator.c is deprecated", DeprecationWarning)
if not self.fixed_material:
raise RuntimeError("Cannot compute speed of light "
"for non-constant material")
return 1/(self.mu*self.epsilon)**0.5
def flux(self, flux_type):
def flux(self, w):
"""The template for the numerical flux for variable coefficients.
:param flux_type: can be in [0,1] for anything between central and upwind,
......@@ -124,88 +119,72 @@ class MaxwellOperator(HyperbolicOperator):
As per Hesthaven and Warburton page 433.
"""
from grudge.flux import (make_normal, FluxVectorPlaceholder,
FluxConstantPlaceholder)
normal = make_normal(self.dimensions)
normal = sym.normal(w.dd, self.dimensions)
if self.fixed_material:
from grudge.tools import count_subset
w = FluxVectorPlaceholder(count_subset(self.get_eh_subset()))
e, h = self.split_eh(w)
epsilon = FluxConstantPlaceholder(self.epsilon)
mu = FluxConstantPlaceholder(self.mu)
else:
from grudge.tools import count_subset
w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())+2)
epsilon = self.epsilon
mu = self.mu
epsilon, mu, e, h = self.split_eps_mu_eh(w)
Z_int = (mu.int/epsilon.int)**0.5
Z_int = (mu/epsilon)**0.5
Y_int = 1/Z_int
Z_ext = (mu.ext/epsilon.ext)**0.5
Z_ext = (mu/epsilon)**0.5
Y_ext = 1/Z_ext
if flux_type == "lf":
if self.flux_type == "lf":
if self.fixed_material:
max_c = (self.epsilon*self.mu)**(-0.5)
else:
from grudge.flux import Max
c_int = (epsilon.int*mu.int)**(-0.5)
c_ext = (epsilon.ext*mu.ext)**(-0.5)
max_c = Max(c_int, c_ext) # noqa
return join_fields(
# flux e,
1/2*(
-self.space_cross_h(normal, h.int-h.ext)
# multiplication by epsilon undoes material divisor below
#-max_c*(epsilon.int*e.int - epsilon.ext*e.ext)
#-max_c*(epsilon*e.int - epsilon*e.ext)
),
# flux h
1/2*(
self.space_cross_e(normal, e.int-e.ext)
# multiplication by mu undoes material divisor below
#-max_c*(mu.int*h.int - mu.ext*h.ext)
#-max_c*(mu*h.int - mu*h.ext)
))
elif isinstance(flux_type, (int, float)):
elif isinstance(self.flux_type, (int, float)):
# see doc/maxima/maxwell.mac
return join_fields(
# flux e,
(
-1/(Z_int+Z_ext)*self.space_cross_h(normal,
Z_ext*(h.int-h.ext)
- flux_type*self.space_cross_e(normal, e.int-e.ext))
- self.flux_type*self.space_cross_e(normal, e.int-e.ext))
),
# flux h
(
1/(Y_int + Y_ext)*self.space_cross_e(normal,
Y_ext*(e.int-e.ext)
+ flux_type*self.space_cross_h(normal, h.int-h.ext))
+ self.flux_type*self.space_cross_h(normal, h.int-h.ext))
),
)
else:
raise ValueError("maxwell: invalid flux_type (%s)"
% self.flux_type)
def local_derivatives(self, w=None):
def local_derivatives(self, w):
"""Template for the spatial derivatives of the relevant components of
:math:`E` and :math:`H`
"""
e, h = self.split_eh(self.field_placeholder(w))
e, h = self.split_eh(w)
nabla = sym.nabla(self.dimensions)
def e_curl(field):
return self.space_cross_e(nabla, field)
def h_curl(field):
return self.space_cross_h(nabla, field)
from grudge.symbolic import make_nabla
nabla = make_nabla(self.dimensions)
# in conservation form: u_t + A u_x = 0
return join_fields(
......@@ -213,62 +192,44 @@ class MaxwellOperator(HyperbolicOperator):
e_curl(e)
)
def field_placeholder(self, w=None):
"A placeholder for E and H."
from grudge.tools import count_subset
fld_cnt = count_subset(self.get_eh_subset())
if w is None:
from grudge.symbolic import make_sym_vector
w = make_sym_vector("w", fld_cnt)
return w
def pec_bc(self, w=None):
def pec_bc(self, w):
"Construct part of the flux operator template for PEC boundary conditions"
e, h = self.split_eh(self.field_placeholder(w))
e, h = self.split_eh(w)
pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e))
pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h))
from grudge.symbolic import RestrictToBoundary
pec_e = RestrictToBoundary(self.pec_tag)(e)
pec_h = RestrictToBoundary(self.pec_tag)(h)
return join_fields(-pec_e, pec_h)
def pmc_bc(self, w=None):
def pmc_bc(self, w):
"Construct part of the flux operator template for PMC boundary conditions"
e, h = self.split_eh(self.field_placeholder(w))
e, h = self.split_eh(w)
from grudge.symbolic import RestrictToBoundary
pmc_e = RestrictToBoundary(self.pmc_tag)(e)
pmc_h = RestrictToBoundary(self.pmc_tag)(h)
pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e))
pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h))
return join_fields(pmc_e, -pmc_h)
def absorbing_bc(self, w=None):
def absorbing_bc(self, w):
"""Construct part of the flux operator template for 1st order
absorbing boundary conditions.
"""
from grudge.symbolic import normal
absorb_normal = normal(self.absorb_tag, self.dimensions)
absorb_normal = sym.normal(self.absorb_tag, self.dimensions)
from grudge.symbolic import RestrictToBoundary, Field
e, h = self.split_eh(self.field_placeholder(w))
e, h = self.split_eh(w)
if self.fixed_material:
epsilon = self.epsilon
mu = self.mu
else:
epsilon = cse(
RestrictToBoundary(self.absorb_tag)(Field("epsilon")))
mu = cse(
RestrictToBoundary(self.absorb_tag)(Field("mu")))
absorb_Z = (mu/epsilon)**0.5
absorb_Y = 1/absorb_Z
absorb_e = RestrictToBoundary(self.absorb_tag)(e)
absorb_h = RestrictToBoundary(self.absorb_tag)(h)
absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e))
absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h))
bc = join_fields(
absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e(
......@@ -281,17 +242,12 @@ class MaxwellOperator(HyperbolicOperator):
return bc
def incident_bc(self, w=None):
def incident_bc(self, w):
"Flux terms for incident boundary conditions"
# NOTE: Untested for inhomogeneous materials, but would usually be
# physically meaningless anyway (are there exceptions to this?)
e, h = self.split_eh(self.field_placeholder(w))
if not self.fixed_material:
from warnings import warn
if self.incident_tag != BTAG_NONE:
warn("Incident boundary conditions assume homogeneous"
" background material, results may be unphysical")
e, h = self.split_eh(w)
from grudge.tools import count_subset
fld_cnt = count_subset(self.get_eh_subset())
......@@ -301,7 +257,7 @@ class MaxwellOperator(HyperbolicOperator):
if is_zero(incident_bc_data):
return make_obj_array([0]*fld_cnt)
else:
return cse(-incident_bc_data)
return sym.cse(-incident_bc_data)
def sym_operator(self, w=None):
"""The full operator template - the high level description of
......@@ -310,23 +266,9 @@ class MaxwellOperator(HyperbolicOperator):
Combines the relevant operator templates for spatial
derivatives, flux, boundary conditions etc.
"""
w = self.field_placeholder(w)
if self.fixed_material:
flux_w = w
else:
epsilon = self.epsilon
mu = self.mu
flux_w = join_fields(epsilon, mu, w)
from grudge.symbolic import BoundaryPair, \
InverseMassOperator, get_flux_operator
flux_op = get_flux_operator(self.flux(self.flux_type))
bdry_flux_op = get_flux_operator(self.flux(self.bdry_flux_type))
from grudge.tools import count_subset
w = sym.make_sym_array("w", count_subset(self.get_eh_subset()))
from grudge.tools.indexing import count_subset
elec_components = count_subset(self.get_eh_subset()[0:3])
mag_components = count_subset(self.get_eh_subset()[3:6])
......@@ -334,10 +276,6 @@ class MaxwellOperator(HyperbolicOperator):
# need to check this
material_divisor = (
[self.epsilon]*elec_components+[self.mu]*mag_components)
else:
material_divisor = join_fields(
[epsilon]*elec_components,
[mu]*mag_components)
tags_and_bcs = [
(self.pec_tag, self.pec_bc(w)),
......@@ -346,67 +284,20 @@ class MaxwellOperator(HyperbolicOperator):
(self.incident_tag, self.incident_bc(w)),
]
def make_flux_bc_vector(tag, bc):
if self.fixed_material:
return bc
else:
from grudge.symbolic import RestrictToBoundary
return join_fields(
cse(RestrictToBoundary(tag)(epsilon)),
cse(RestrictToBoundary(tag)(mu)),
bc)
def flux(pair):
return sym.interp(pair.dd, "all_faces")(self.flux(pair))
return (
- self.local_derivatives(w)
+ InverseMassOperator()(
flux_op(flux_w)
- sym.InverseMassOperator()(sym.FaceMassOperator()(
flux(sym.int_tpair(w))
+ sum(
bdry_flux_op(BoundaryPair(
flux_w, make_flux_bc_vector(tag, bc), tag))
for tag, bc in tags_and_bcs))
) / material_divisor
flux(sym.bv_tpair(tag, w, bc))
for tag, bc in tags_and_bcs)
) )) / material_divisor
def bind(self, discr):
"Convert the abstract operator template into compiled code."
from grudge.mesh import check_bc_coverage
check_bc_coverage(discr.mesh, [
self.pec_tag, self.absorb_tag, self.incident_tag])
compiled_sym_operator = discr.compile(self.op_template())
def rhs(t, w, **extra_context):
kwargs = {}
kwargs.update(extra_context)
return compiled_sym_operator(w=w, t=t, **kwargs)
return rhs
def assemble_eh(self, e=None, h=None, discr=None):
"Combines separate E and H vectors into a single array."
if discr is None:
def zero():
return 0
else:
def zero():
return discr.volume_zeros()
from grudge.tools import count_subset
e_components = count_subset(self.get_eh_subset()[0:3])
h_components = count_subset(self.get_eh_subset()[3:6])
def default_fld(fld, comp):
if fld is None:
return [zero() for i in range(comp)]
else:
return fld
e = default_fld(e, e_components)
h = default_fld(h, h_components)
return join_fields(e, h)
assemble_fields = assemble_eh
@memoize_method
def partial_to_eh_subsets(self):
......@@ -445,11 +336,8 @@ class MaxwellOperator(HyperbolicOperator):
e_idx, h_idx = self.partial_to_eh_subsets()
e, h = w[e_idx], w[h_idx]
from grudge.flux import FluxVectorPlaceholder as FVP
if isinstance(w, FVP):
return FVP(scalars=e), FVP(scalars=h)
else:
return make_obj_array(e), make_obj_array(h)
return e,h
#return make_obj_array(e), make_obj_array(h)
def get_eh_subset(self):
"""Return a 6-tuple of :class:`bool` objects indicating whether field
......@@ -476,6 +364,14 @@ class MaxwellOperator(HyperbolicOperator):
raise ValueError("max_eigenvalue is no longer supported for "
"variable-coefficient problems--use max_eigenvalue_expr")
def check_bc_coverage(self, mesh):
from meshmode.mesh import check_bc_coverage
check_bc_coverage(mesh, [
self.pec_tag,
self.pmc_tag,
self.absorb_tag,
self.incident_tag])
class TMMaxwellOperator(MaxwellOperator):
"""A 2D TM Maxwell operator with PEC boundaries.
......
grudge/models/wave.py deleted 100644 → 0
View file @ c5e574c0
# -*- coding: utf8 -*-
"""Wave equation operators."""
from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import numpy as np
from grudge.models import HyperbolicOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
from grudge import sym
from pytools.obj_array import join_fields
# {{{ constant-velocity
class StrongWaveOperator(HyperbolicOperator):
"""This operator discretizes the wave equation
:math:`\\partial_t^2 u = c^2 \\Delta u`.
To be precise, we discretize the hyperbolic system
.. math::
\partial_t u - c \\nabla \\cdot v = 0
\partial_t v - c \\nabla u = 0
The sign of :math:`v` determines whether we discretize the forward or the
backward wave equation.
:math:`c` is assumed to be constant across all space.
"""
def __init__(self, c, ambient_dim, source_f=0,
flux_type="upwind",
dirichlet_tag=BTAG_ALL,
dirichlet_bc_f=0,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_NONE):
assert isinstance(ambient_dim, int)
self.c = c
self.ambient_dim = ambient_dim
self.source_f = source_f
if self.c > 0:
self.sign = 1
else:
self.sign = -1
self.dirichlet_tag = dirichlet_tag
self.neumann_tag = neumann_tag
self.radiation_tag = radiation_tag
self.dirichlet_bc_f = dirichlet_bc_f
self.flux_type = flux_type
def flux(self, w):
u = w[0]
v = w[1:]
normal = sym.normal(w.dd, self.ambient_dim)
flux_weak = join_fields(
np.dot(v.avg, normal),
u.avg * normal)
if self.flux_type == "central":
pass
elif self.flux_type == "upwind":
# see doc/notes/grudge-notes.tm
flux_weak -= self.sign*join_fields(
0.5*(u.int-u.ext),
0.5*(normal * np.dot(normal, v.int-v.ext)))
else:
raise ValueError("invalid flux type '%s'" % self.flux_type)
flux_strong = join_fields(
np.dot(v.int, normal),
u.int * normal) - flux_weak
return self.c*flux_strong
def sym_operator(self):
d = self.ambient_dim
w = sym.make_sym_array("w", d+1)
u = w[0]
v = w[1:]
# boundary conditions -------------------------------------------------
# dirichlet BCs -------------------------------------------------------
dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
if self.dirichlet_bc_f:
# FIXME
from warnings import warn
warn("Inhomogeneous Dirichlet conditions on the wave equation "
"are still having issues.")
dir_g = sym.Field("dir_bc_u")
dir_bc = join_fields(2*dir_g - dir_u, dir_v)
else:
dir_bc = join_fields(-dir_u, dir_v)
dir_bc = sym.cse(dir_bc, "dir_bc")
# neumann BCs ---------------------------------------------------------
neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
# radiation BCs -------------------------------------------------------
rad_normal = sym.normal(self.radiation_tag, d)
rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
rad_bc = sym.cse(join_fields(
0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
), "rad_bc")
# entire operator -----------------------------------------------------
nabla = sym.nabla(d)
def flux(pair):
return sym.interp(pair.dd, "all_faces")(
self.flux(pair))
result = (
- join_fields(
-self.c*np.dot(nabla, v),
-self.c*(nabla*u)
)
+
sym.InverseMassOperator()(
sym.FaceMassOperator()(
flux(sym.int_tpair(w))
+ flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
+ flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
+ flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
)))
result[0] += self.source_f
return result
def check_bc_coverage(self, mesh):
from meshmode.mesh import check_bc_coverage
check_bc_coverage(mesh, [
self.dirichlet_tag,
self.neumann_tag,
self.radiation_tag])
def max_eigenvalue(self, t, fields=None, discr=None):
return abs(self.c)
class WeakWaveOperator(HyperbolicOperator):
"""This operator discretizes the wave equation
:math:`\\partial_t^2 u = c^2 \\Delta u`.
To be precise, we discretize the hyperbolic system
.. math::
\partial_t u - c \\nabla \\cdot v = 0
\partial_t v - c \\nabla u = 0
The sign of :math:`v` determines whether we discretize the forward or the
backward wave equation.
:math:`c` is assumed to be constant across all space.
"""
def __init__(self, c, ambient_dim, source_f=0,
flux_type="upwind",
dirichlet_tag=BTAG_ALL,
dirichlet_bc_f=0,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_NONE):
assert isinstance(ambient_dim, int)
self.c = c
self.ambient_dim = ambient_dim
self.source_f = source_f
if self.c > 0:
self.sign = 1
else:
self.sign = -1
self.dirichlet_tag = dirichlet_tag
self.neumann_tag = neumann_tag
self.radiation_tag = radiation_tag
self.dirichlet_bc_f = dirichlet_bc_f
self.flux_type = flux_type
def flux(self, w):
u = w[0]
v = w[1:]
normal = sym.normal(w.dd, self.ambient_dim)
flux_weak = join_fields(
np.dot(v.avg, normal),
u.avg * normal)
if self.flux_type == "central":
return -self.c*flux_weak
elif self.flux_type == "upwind":
return -self.c*(flux_weak + self.sign*join_fields(
0.5*(u.int-u.ext),
0.5*(normal * np.dot(normal, v.int-v.ext))))
else:
raise ValueError("invalid flux type '%s'" % self.flux_type)
def sym_operator(self):
d = self.ambient_dim
w = sym.make_sym_array("w", d+1)
u = w[0]
v = w[1:]
# boundary conditions -------------------------------------------------
# dirichlet BCs -------------------------------------------------------
dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
if self.dirichlet_bc_f:
# FIXME
from warnings import warn
warn("Inhomogeneous Dirichlet conditions on the wave equation "
"are still having issues.")
dir_g = sym.Field("dir_bc_u")
dir_bc = join_fields(2*dir_g - dir_u, dir_v)
else:
dir_bc = join_fields(-dir_u, dir_v)
dir_bc = sym.cse(dir_bc, "dir_bc")
# neumann BCs ---------------------------------------------------------
neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
# radiation BCs -------------------------------------------------------
rad_normal = sym.normal(self.radiation_tag, d)
rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
rad_bc = sym.cse(join_fields(
0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
), "rad_bc")
# entire operator -----------------------------------------------------
def flux(pair):
return sym.interp(pair.dd, "all_faces")(self.flux(pair))
result = sym.InverseMassOperator()(
join_fields(
-self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
-self.c*(sym.stiffness_t(self.ambient_dim)*u)
)
- sym.FaceMassOperator()(flux(sym.int_tpair(w))
+ flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
+ flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
+ flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
))
result[0] += self.source_f
return result
def check_bc_coverage(self, mesh):
from meshmode.mesh import check_bc_coverage
check_bc_coverage(mesh, [
self.dirichlet_tag,
self.neumann_tag,
self.radiation_tag])
def max_eigenvalue(self, t, fields=None, discr=None):
return abs(self.c)
# }}}
# {{{ variable-velocity
class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
"""This operator discretizes the wave equation
:math:`\\partial_t^2 u = c^2 \\Delta u`.
To be precise, we discretize the hyperbolic system
.. math::
\partial_t u - c \\nabla \\cdot v = 0
\partial_t v - c \\nabla u = 0
The sign of :math:`v` determines whether we discretize the forward or the
backward wave equation.
:math:`c` is assumed to be constant across all space.
"""
def __init__(self, c, ambient_dim, source_f=0,
flux_type="upwind",
dirichlet_tag=BTAG_ALL,
dirichlet_bc_f=0,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_NONE):
assert isinstance(ambient_dim, int)
self.c = c
self.ambient_dim = ambient_dim
self.source_f = source_f
self.sign = sym.If(sym.Comparison(
self.c, ">", 0),
np.int32(1), np.int32(-1))
self.dirichlet_tag = dirichlet_tag
self.neumann_tag = neumann_tag
self.radiation_tag = radiation_tag
self.dirichlet_bc_f = dirichlet_bc_f
self.flux_type = flux_type
def flux(self, w):
c = w[0]
u = w[1]
v = w[2:]
normal = sym.normal(w.dd, self.ambient_dim)
if self.flux_type == "central":
return -0.5 * join_fields(
np.dot(v.int*c.int + v.ext*c.ext, normal),
(u.int * c.int + u.ext*c.ext) * normal)
elif self.flux_type == "upwind":
return -0.5 * join_fields(np.dot(normal, c.ext * v.ext + c.int * v.int)
+ c.ext*u.ext - c.int * u.int,
normal * (np.dot(normal, c.ext * v.ext - c.int * v.int)
+ c.ext*u.ext + c.int * u.int))
else:
raise ValueError("invalid flux type '%s'" % self.flux_type)
def sym_operator(self):
d = self.ambient_dim
w = sym.make_sym_array("w", d+1)
u = w[0]
v = w[1:]
flux_w = join_fields(self.c, w)
# boundary conditions -------------------------------------------------
# dirichlet BCs -------------------------------------------------------
dir_c = sym.cse(sym.interp("vol", self.dirichlet_tag)(self.c))
dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
if self.dirichlet_bc_f:
# FIXME
from warnings import warn
warn("Inhomogeneous Dirichlet conditions on the wave equation "
"are still having issues.")
dir_g = sym.Field("dir_bc_u")
dir_bc = join_fields(dir_c, 2*dir_g - dir_u, dir_v)
else:
dir_bc = join_fields(dir_c, -dir_u, dir_v)
dir_bc = sym.cse(dir_bc, "dir_bc")
# neumann BCs ---------------------------------------------------------
neu_c = sym.cse(sym.interp("vol", self.neumann_tag)(self.c))
neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
neu_bc = sym.cse(join_fields(neu_c, neu_u, -neu_v), "neu_bc")
# radiation BCs -------------------------------------------------------
rad_normal = sym.normal(self.radiation_tag, d)
rad_c = sym.cse(sym.interp("vol", self.radiation_tag)(self.c))
rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
rad_bc = sym.cse(join_fields(rad_c,
0.5*(rad_u - sym.interp("vol", self.radiation_tag)(self.sign)
* np.dot(rad_normal, rad_v)),
0.5*rad_normal*(np.dot(rad_normal, rad_v)
- sym.interp("vol", self.radiation_tag)(self.sign)*rad_u)
), "rad_bc")
# entire operator -----------------------------------------------------
def flux(pair):
return sym.interp(pair.dd, "all_faces")(self.flux(pair))
result = sym.InverseMassOperator()(
join_fields(
-self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
-self.c*(sym.stiffness_t(self.ambient_dim)*u)
)
- sym.FaceMassOperator()(flux(sym.int_tpair(flux_w))
+ flux(sym.bv_tpair(self.dirichlet_tag, flux_w, dir_bc))
+ flux(sym.bv_tpair(self.neumann_tag, flux_w, neu_bc))
+ flux(sym.bv_tpair(self.radiation_tag, flux_w, rad_bc))
))
result[0] += self.source_f
return result
def check_bc_coverage(self, mesh):
from meshmode.mesh import check_bc_coverage
check_bc_coverage(mesh, [
self.dirichlet_tag,
self.neumann_tag,
self.radiation_tag])
def max_eigenvalue(self, t, fields=None, discr=None):
return sym.NodalMax()(sym.CFunction("fabs")(self.c))
# }}}
# vim: foldmethod=marker
grudge/shortcuts.py
View file @ c2196857
......@@ -43,6 +43,21 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
return dt_stepper
def set_up_forward_euler(field_var_name, dt, fields, rhs, t_start=0):
from leap.rk import ForwardEulerMethod
from dagrt.codegen import PythonCodeGenerator
dt_method = ForwardEulerMethod(component_id=field_var_name)
dt_code = dt_method.generate()
dt_stepper_class = PythonCodeGenerator("TimeStep").get_class(dt_code)
dt_stepper = dt_stepper_class({"<func>"+dt_method.component_id: rhs})
dt_stepper.set_up(
t_start=t_start, dt_start=dt,
context={dt_method.component_id: fields})
return dt_stepper
def make_visualizer(discr, vis_order):
from meshmode.discretization.visualization import make_visualizer
......
grudge/symbolic/compiler.py
View file @ c2196857
......@@ -852,7 +852,11 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper):
def map_call(self, expr):
if isinstance(expr.function, sym.CFunction):
from pymbolic import var
return var(expr.function.name)(
func_name = expr.function.name
if func_name == "fabs":
func_name = "abs"
return var(func_name)(
*[self.rec(par) for par in expr.parameters])
else:
raise NotImplementedError(
......
grudge/symbolic/mappers/__init__.py
View file @ c2196857
......@@ -34,6 +34,7 @@ import pymbolic.mapper.evaluator
import pymbolic.mapper.dependency
import pymbolic.mapper.substitutor
import pymbolic.mapper.constant_folder
import pymbolic.mapper.constant_converter
import pymbolic.mapper.flop_counter
from pymbolic.mapper import CSECachingMapperMixin
......@@ -812,6 +813,12 @@ class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper):
# {{{ simplification / optimization
class ConstantToNumpyConversionMapper(
pymbolic.mapper.constant_converter.ConstantToNumpyConversionMapper,
IdentityMapperMixin):
pass
class CommutativeConstantFoldingMapper(
pymbolic.mapper.constant_folder.CommutativeConstantFoldingMapper,
IdentityMapperMixin):
......
grudge/tools.py
View file @ c2196857
......@@ -36,3 +36,105 @@ def is_zero(x):
return x == 0
else:
return False
def levi_civita(tuple):
"""Compute an entry of the Levi-Civita tensor for the indices *tuple*.
Only three-tuples are supported for now.
"""
if len(tuple) == 3:
if tuple in [(0,1,2), (2,0,1), (1,2,0)]:
return 1
elif tuple in [(2,1,0), (0,2,1), (1,0,2)]:
return -1
else:
return 0
else:
raise NotImplementedError
class SubsettableCrossProduct:
"""A cross product that can operate on an arbitrary subsets of its
two operands and return an arbitrary subset of its result.
"""
full_subset = (True, True, True)
def __init__(self, op1_subset=full_subset, op2_subset=full_subset, result_subset=full_subset):
"""Construct a subset-able cross product.
:param op1_subset: The subset of indices of operand 1 to be taken into account.
Given as a 3-sequence of bools.
:param op2_subset: The subset of indices of operand 2 to be taken into account.
Given as a 3-sequence of bools.
:param result_subset: The subset of indices of the result that are calculated.
Given as a 3-sequence of bools.
"""
def subset_indices(subset):
return [i for i, use_component in enumerate(subset)
if use_component]
self.op1_subset = op1_subset
self.op2_subset = op2_subset
self.result_subset = result_subset
import pymbolic
op1 = pymbolic.var("x")
op2 = pymbolic.var("y")
self.functions = []
self.component_lcjk = []
for i, use_component in enumerate(result_subset):
if use_component:
this_expr = 0
this_component = []
for j, j_real in enumerate(subset_indices(op1_subset)):
for k, k_real in enumerate(subset_indices(op2_subset)):
lc = levi_civita((i, j_real, k_real))
if lc != 0:
this_expr += lc*op1.index(j)*op2.index(k)
this_component.append((lc, j, k))
self.functions.append(pymbolic.compile(this_expr,
variables=[op1, op2]))
self.component_lcjk.append(this_component)
def __call__(self, x, y, three_mult=None):
"""Compute the subsetted cross product on the indexables *x* and *y*.
:param three_mult: a function of three arguments *sign, xj, yk*
used in place of the product *sign*xj*yk*. Defaults to just this
product if not given.
"""
from pytools.obj_array import join_fields
if three_mult is None:
return join_fields(*[f(x, y) for f in self.functions])
else:
return join_fields(
*[sum(three_mult(lc, x[j], y[k]) for lc, j, k in lcjk)
for lcjk in self.component_lcjk])
cross = SubsettableCrossProduct()
def count_subset(subset):
from pytools import len_iterable
return len_iterable(uc for uc in subset if uc)
def partial_to_all_subset_indices(subsets, base=0):
"""Takes a sequence of bools and generates it into an array of indices
to be used to insert the subset into the full set.
Example:
>>> list(partial_to_all_subset_indices([[False, True, True], [True,False,True]]))
[array([0 1]), array([2 3]
"""
idx = base
for subset in subsets:
result = []
for is_in in subset:
if is_in:
result.append(idx)
idx += 1
yield np.array(result, dtype=np.intp)