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 4191 additions and 2265 deletions
examples/maxwell/maxwell-pml.py deleted 100644 → 0
View file @ 90acbf2e
"""grudge is the Hybrid'n'Easy Discontinuous Galerkin Environment."""
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
__copyright__ = "Copyright (C) 2007 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
def make_mesh(a, b, pml_width=0.25, **kwargs):
from meshpy.geometry import GeometryBuilder, make_circle
geob = GeometryBuilder()
circle_centers = [(-1.5, 0), (1.5, 0)]
for cent in circle_centers:
geob.add_geometry(*make_circle(1, cent))
geob.wrap_in_box(1)
geob.wrap_in_box(pml_width)
mesh_mod = geob.mesher_module()
mi = mesh_mod.MeshInfo()
geob.set(mi)
mi.set_holes(circle_centers)
built_mi = mesh_mod.build(mi, **kwargs)
def boundary_tagger(fvi, el, fn, points):
return []
from grudge.mesh import make_conformal_mesh_ext
from grudge.mesh.element import Triangle
pts = np.asarray(built_mi.points, dtype=np.float64)
return make_conformal_mesh_ext(
pts,
[Triangle(i, el, pts)
for i, el in enumerate(built_mi.elements)],
boundary_tagger)
def main(write_output=True):
from grudge.timestep.runge_kutta import LSRK4TimeStepper
from math import sqrt, pi, exp
from grudge.backends import guess_run_context
rcon = guess_run_context()
epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
mu0 = 4*pi*1e-7 # N/A**2.
epsilon = 1*epsilon0
mu = 1*mu0
c = 1/sqrt(mu*epsilon)
pml_width = 0.5
#mesh = make_mesh(a=np.array((-1,-1,-1)), b=np.array((1,1,1)),
#mesh = make_mesh(a=np.array((-3,-3)), b=np.array((3,3)),
mesh = make_mesh(a=np.array((-1,-1)), b=np.array((1,1)),
#mesh = make_mesh(a=np.array((-2,-2)), b=np.array((2,2)),
pml_width=pml_width, max_volume=0.01)
if rcon.is_head_rank:
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
class Current:
def volume_interpolant(self, t, discr):
from grudge.tools import make_obj_array
result = discr.volume_zeros(kind="numpy", dtype=np.float64)
omega = 6*c
if omega*t > 2*pi:
return make_obj_array([result, result, result])
x = make_obj_array(discr.nodes.T)
r = np.sqrt(np.dot(x, x))
idx = r<0.3
result[idx] = (1+np.cos(pi*r/0.3))[idx] \
*np.sin(omega*t)**3
result = discr.convert_volume(result, kind=discr.compute_kind,
dtype=discr.default_scalar_type)
return make_obj_array([-result, result, result])
order = 3
discr = rcon.make_discretization(mesh_data, order=order,
debug=["cuda_no_plan"])
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "em-%d" % order)
from grudge.mesh import BTAG_ALL, BTAG_NONE
from grudge.data import GivenFunction, TimeHarmonicGivenFunction, TimeIntervalGivenFunction
from grudge.models.em import MaxwellOperator
from grudge.models.pml import \
AbarbanelGottliebPMLMaxwellOperator, \
AbarbanelGottliebPMLTMMaxwellOperator, \
AbarbanelGottliebPMLTEMaxwellOperator
op = AbarbanelGottliebPMLTEMaxwellOperator(epsilon, mu, flux_type=1,
current=Current(),
pec_tag=BTAG_ALL,
absorb_tag=BTAG_NONE,
add_decay=True
)
fields = op.assemble_ehpq(discr=discr)
stepper = LSRK4TimeStepper()
if rcon.is_head_rank:
print("order %d" % order)
print("#elements=", len(mesh.elements))
# 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)
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"])
from grudge.log import LpNorm
class FieldIdxGetter:
def __init__(self, whole_getter, idx):
self.whole_getter = whole_getter
self.idx = idx
def __call__(self):
return self.whole_getter()[self.idx]
# timestep loop -------------------------------------------------------
t = 0
pml_coeff = op.coefficients_from_width(discr, width=pml_width)
rhs = op.bind(discr, pml_coeff)
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=4/c, 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 % 10 == 0 and write_output:
e, h, p, q = op.split_ehpq(fields)
visf = vis.make_file("em-%d-%04d" % (order, step))
#pml_rhs_e, pml_rhs_h, pml_rhs_p, pml_rhs_q = \
#op.split_ehpq(rhs(t, fields))
j = Current().volume_interpolant(t, discr)
vis.add_data(visf, [
("e", discr.convert_volume(e, "numpy")),
("h", discr.convert_volume(h, "numpy")),
("p", discr.convert_volume(p, "numpy")),
("q", discr.convert_volume(q, "numpy")),
("j", discr.convert_volume(j, "numpy")),
#("pml_rhs_e", pml_rhs_e),
#("pml_rhs_h", pml_rhs_h),
#("pml_rhs_p", pml_rhs_p),
#("pml_rhs_q", pml_rhs_q),
#("max_rhs_e", max_rhs_e),
#("max_rhs_h", max_rhs_h),
#("max_rhs_p", max_rhs_p),
#("max_rhs_q", max_rhs_q),
],
time=t, step=step)
visf.close()
fields = stepper(fields, t, dt, rhs)
_, _, energies_data = logmgr.get_expr_dataset("W_el+W_mag")
energies = [value for tick_nbr, value in energies_data]
assert energies[-1] < max(energies) * 1e-2
finally:
logmgr.close()
if write_output:
vis.close()
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
from pytools.test import mark_test
@mark_test.long
def test_maxwell_pml():
main(write_output=False)
examples/maxwell/notes.tm deleted 100644 → 0
View file @ 90acbf2e
<TeXmacs|1.0.6>
<style|<tuple|generic|maxima|axiom>>
<\body>
<section|Cylindrical TM Maxwell Cavity Mode>
<with|prog-language|axiom|prog-session|default|<\session>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)clear all
</input>
<\output>
\ \ \ All user variables and function definitions have been cleared.
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)library )dir "/home/andreas/axiom"
</input>
<\output>
\ \ \ TexFormat is already explicitly exposed in frame initial\
\ \ \ TexFormat will be automatically loaded when needed from\
\ \ \ \ \ \ /home/andreas/axiom/TEX.NRLIB/code
\ \ \ TexFormat1 is already explicitly exposed in frame initial\
\ \ \ TexFormat1 will be automatically loaded when needed from\
\ \ \ \ \ \ /home/andreas/axiom/TEX1.NRLIB/code
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
J:=operator 'J
</input>
<\output>
<with|mode|math|math-display|true|J<leqno>(1)>
<axiomtype|BasicOperator >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
psi(rho,phi) == J(gamma*rho)*exp(PP*%i*m*phi)
</input>
<\output>
<axiomtype|Void >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
psiglob:=psi(sqrt(x^2+y^2),atan(y/x))
</input>
<\output>
\ \ \ Compiling function psi with type (Expression Integer,Expression\
\ \ \ \ \ \ Integer) -\<gtr\> Expression Complex Integer\
<with|mode|math|math-display|true|J<left|(>\<gamma\><sqrt|y<rsup|2>+x<rsup|2>><right|)>e<rsup|<left|(>i*P*P*m*arctan
<left|(><frac|y|x><right|)><right|)>><leqno>(3)>
<axiomtype|Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
D(psi(rho,phi),rho)
</input>
<\output>
\ \ \ Compiling function psi with type (Variable rho,Variable phi)
-\<gtr\>\
\ \ \ \ \ \ Expression Complex Integer\
<with|mode|math|math-display|true|\<gamma\>e<rsup|<left|(>i*P*P*m\<phi\><right|)>>J<rsub|
><rsup|,><left|(>\<gamma\>\<rho\><right|)><leqno>(5)>
<axiomtype|Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
cross(vector [0,0,1], vector [x,y,0])
</input>
<\output>
<with|mode|math|math-display|true|<left|[>-y,<space|0.5spc>x,<space|0.5spc>0<right|]><leqno>(7)>
<axiomtype|Vector Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
\;
</input>
</session>>
<section|Rectangular Cavity Mode>
According to Jackson, p. 357, (8.17), we need to solve the Helmholtz
equation
<\eqnarray*>
<tformat|<cwith|1|1|3|3|cell-halign|l>|<table|<row|<cell|(\<nabla\><rsup|2>+\<mu\>\<varepsilon\>\<omega\><rsup|2>)<matrix|<tformat|<table|<row|<cell|\<b-E\>>>|<row|<cell|\<b-B\>>>>>>>|<cell|=>|<cell|\<b-0\>,>>>>
</eqnarray*>
subject to <with|mode|math|n\<times\>\<b-E\>=0> and
<with|mode|math|n\<cdot\>\<b-B\>=0>. The ansatz is
<\equation*>
\<b-E\>=<matrix|<tformat|<table|<row|<cell|E<rsub|x,x>(x)E<rsub|x,y>(y)E<rsub|x,z>(z)>>|<row|<cell|E<rsub|y,x>(x)E<rsub|y,y>(y)E<rsub|y,z>(z)>>|<row|<cell|E<rsub|z,x>(x)E<rsub|z,y>(y)E<rsub|z,z>(z)>>>>>
</equation*>
and likewise for <with|mode|math|\<b-B\>>. The boundary conditions are
<\eqnarray*>
<tformat|<table|<row|<cell|E<rsub|x>(x,<with|math-level|1|<tabular|<tformat|<table|<row|<cell|0>>|<row|<cell|b>>>>>>,z)>|<cell|=>|<cell|0,>>|<row|<cell|E<rsub|x>(x,y,<with|math-level|1|<tabular|<tformat|<table|<row|<cell|0>>|<row|<cell|c>>>>>>)>|<cell|=>|<cell|0,>>>>
</eqnarray*>
and so on, as well as
<\eqnarray*>
<tformat|<table|<row|<cell|H<rsub|x>(<with|math-level|1|<tabular|<tformat|<table|<row|<cell|0>>|<row|<cell|a>>>>>>,y,z)>|<cell|=>|<cell|0.>>>>
</eqnarray*>
So
<\equation*>
E<rsub|x>=\<alpha\><rsub|x>exp(i\<beta\><rsub|x>x)sin<left|(><frac|n\<pi\>y|b><right|)>sin<left|(><frac|o\<pi\>z|c><right|)>exp(-i\<omega\>t)=\<alpha\><rsub|x>e<rsub|x>s<rsub|y>s<rsub|z>
</equation*>
and analogous terms for <with|mode|math|E<rsub|y>> and
<with|mode|math|E<rsub|z>> satisfy the first batch of boundary conditions.
Because of the Helmholtz equation, we find that
<with|mode|math|\<beta\><rsub|x>=m\<pi\>/a>; otherwise, not all vector
components would share the same eigenvalue, which would not solve the
equation.
<with|prog-language|axiom|prog-session|default|<\session>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)clear all
</input>
<\output>
\ \ \ All user variables and function definitions have been cleared.
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)library )dir "/home/andreas/axiom"
</input>
<\output>
\ \ \ TexFormat is already explicitly exposed in frame initial\
\ \ \ TexFormat will be automatically loaded when needed from\
\ \ \ \ \ \ /home/andreas/axiom/TEX.NRLIB/code
\ \ \ TexFormat1 is already explicitly exposed in frame initial\
\ \ \ TexFormat1 will be automatically loaded when needed from\
\ \ \ \ \ \ /home/andreas/axiom/TEX1.NRLIB/code
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
factors:=[f,g,h];
</input>
<\output>
<axiomtype|List OrderedVariableList [f,g,h] >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
coord := [x,y,z];
</input>
<\output>
<axiomtype|List OrderedVariableList [x,y,z] >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
curl(v)== vector [D(v.3,y)-D(v.2,z),D(v.1,z)-D(v.3,x),D(v.2,x)-D(v.1,y)];
</input>
<\output>
<axiomtype|Void >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
c:=1/sqrt(epsilon*mu);
</input>
<\output>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
sines(i) == sin(factors.i*coord.i);
</input>
<\output>
<axiomtype|Void >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
cosines(i) == cos(factors.i*coord.i);
</input>
<\output>
<axiomtype|Void >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
k:=sqrt(f^2+g^2+h^2);omega:=k*c;
</input>
<\output>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
zdep1:=exp(%i*h*z); zdep2:=exp(-%i*h*z);
</input>
<\output>
<axiomtype|Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
zf1:=1; zf2:=-1;
</input>
<\output>
<axiomtype|Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
zdep:=zf1*zdep1 + zf2*zdep2;
</input>
<\output>
<axiomtype|Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
C:=%i/(f^2+g^2);
</input>
<\output>
<axiomtype|Fraction Polynomial Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
efield := vector [
C*f*h*cosines(1)* \ sines(2)*(zf1*zdep1-zf2*zdep2),
C*g*h* \ sines(1)*cosines(2)*(zf1*zdep1-zf2*zdep2),
\ \ \ \ \ \ \ \ sines(1)* \ sines(2)*zdep];
</input>
<\output>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
hfield:=1/(-%i*omega*mu)*(-curl efield);
</input>
<\output>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
efield2:=1/(-%i*omega*epsilon)*(curl hfield);
</input>
<\output>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
efield2-efield
</input>
<\output>
<with|mode|math|math-display|true|<left|[>0,<space|0.5spc>0,<space|0.5spc>0<right|]><leqno>(71)>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
hfield
</input>
<\output>
<with|mode|math|math-display|true|<left|[><frac|<left|(><left|(>-i*g*h<rsup|2>-i*g<rsup|3>-i*f<rsup|2>g<right|)>cos
<left|(>g*y<right|)>e<rsup|<left|(>i*h*z<right|)>>+<left|(>i*g*h<rsup|2>+i*g<rsup|3>+i*f<rsup|2>g<right|)>cos
<left|(>g*y<right|)>e<rsup|<left|(>-i*h*z<right|)>><right|)>sin
<left|(>f*x<right|)><sqrt|\<epsilon\>\<mu\>>|<left|(>g<rsup|2>+f<rsup|2><right|)>\<mu\><sqrt|h<rsup|2>+g<rsup|2>+f<rsup|2>>>,<space|0.5spc><frac|<left|(><left|(>i*f*h<rsup|2>+i*f*g<rsup|2>+i*f<rsup|3><right|)>cos
<left|(>f*x<right|)>e<rsup|<left|(>i*h*z<right|)>>+<left|(>-i*f*h<rsup|2>-i*f*g<rsup|2>-i*f<rsup|3><right|)>cos
<left|(>f*x<right|)>e<rsup|<left|(>-i*h*z<right|)>><right|)>sin
<left|(>g*y<right|)><sqrt|\<epsilon\>\<mu\>>|<left|(>g<rsup|2>+f<rsup|2><right|)>\<mu\><sqrt|h<rsup|2>+g<rsup|2>+f<rsup|2>>>,<space|0.5spc>0<right|]><leqno>(72)>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
hfield2:=vector [
-%i*g/(f^2+g^2)*epsilon*omega*sines(1)*cosines(2)*(zf1*zdep1+zf2*zdep2),
\ %i*f/(f^2+g^2)*epsilon*omega*cosines(1)*sines(2)*(zf1*zdep1+zf2*zdep2),
0]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><frac|<left|(>-i\<epsilon\>g*cos
<left|(>g*y<right|)>e<rsup|<left|(>i*h*z<right|)>>+i\<epsilon\>g*cos
<left|(>g*y<right|)>e<rsup|<left|(>-i*h*z<right|)>><right|)>sin
<left|(>f*x<right|)><sqrt|h<rsup|2>+g<rsup|2>+f<rsup|2>>|<left|(>g<rsup|2>+f<rsup|2><right|)><sqrt|\<epsilon\>\<mu\>>>,<space|0.5spc><frac|<left|(>i\<epsilon\>f*cos
<left|(>f*x<right|)>e<rsup|<left|(>i*h*z<right|)>>-i\<epsilon\>f*cos
<left|(>f*x<right|)>e<rsup|<left|(>-i*h*z<right|)>><right|)>sin
<left|(>g*y<right|)><sqrt|h<rsup|2>+g<rsup|2>+f<rsup|2>>|<left|(>g<rsup|2>+f<rsup|2><right|)><sqrt|\<epsilon\>\<mu\>>>,<space|0.5spc>0<right|]><leqno>(73)>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
hfield-hfield2
</input>
<\output>
<with|mode|math|math-display|true|<left|[>0,<space|0.5spc>0,<space|0.5spc>0<right|]><leqno>(74)>
<axiomtype|Vector Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
bcs:=[
eval(efield.1, z=0),
eval(efield.2, z=0),
eval(efield.3, z=0),
eval(hfield.1, x=0),
eval(hfield.2, y=0),
eval(hfield.3, z=0)
]
</input>
<\output>
<with|mode|math|math-display|true|<left|[>0,<space|0.5spc><frac|2i*g*h*cos
<left|(>g*y<right|)>sin <left|(>f*x<right|)>|g<rsup|2>+f<rsup|2>>,<space|0.5spc><frac|2i*f*h*cos
<left|(>f*x<right|)>sin <left|(>g*y<right|)>|g<rsup|2>+f<rsup|2>>,<space|0.5spc>0,<space|0.5spc>0,<space|0.5spc>0<right|]><leqno>(76)>
<axiomtype|List Expression Complex Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
\;
</input>
</session>>
\;
</body>
<\initial>
<\collection>
<associate|page-type|letter>
</collection>
</initial>
<\references>
<\collection>
<associate|auto-1|<tuple|1|1>>
<associate|auto-2|<tuple|2|1>>
<associate|auto-3|<tuple|3|?>>
</collection>
</references>
<\auxiliary>
<\collection>
<\associate|toc>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|Cylindrical
TM Maxwell Cavity Mode> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-1><vspace|0.5fn>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|Rectangular
Cavity Mode> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-2><vspace|0.5fn>
</associate>
</collection>
</auxiliary>
\ No newline at end of file
examples/poisson/helmholtz.py deleted 100644 → 0
View file @ 90acbf2e
# 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 grudge.tools import Reflection, Rotation
class ResidualPrinter:
def __init__(self, compute_resid=None):
self.count = 0
self.compute_resid = compute_resid
def __call__(self, cur_sol):
import sys
if cur_sol is not None:
if self.count % 20 == 0:
sys.stdout.write("IT %8d %g \r" % (
self.count, la.norm(self.compute_resid(cur_sol))))
else:
sys.stdout.write("IT %8d \r" % self.count)
self.count += 1
sys.stdout.flush()
def main(write_output=True):
from grudge.data import GivenFunction, ConstantGivenFunction
from grudge.backends import guess_run_context
rcon = guess_run_context()
dim = 2
def boundary_tagger(fvi, el, fn, points):
from math import atan2, pi
normal = el.face_normals[fn]
if -90/180*pi < atan2(normal[1], normal[0]) < 90/180*pi:
return ["neumann"]
else:
return ["dirichlet"]
def dirichlet_boundary_tagger(fvi, el, fn, points):
return ["dirichlet"]
if dim == 2:
if rcon.is_head_rank:
from grudge.mesh.generator import make_disk_mesh
mesh = make_disk_mesh(r=0.5,
boundary_tagger=dirichlet_boundary_tagger,
max_area=1e-3)
elif dim == 3:
if rcon.is_head_rank:
from grudge.mesh.generator import make_ball_mesh
mesh = make_ball_mesh(max_volume=0.0001,
boundary_tagger=lambda fvi, el, fn, points:
["dirichlet"])
else:
raise RuntimeError("bad number of dimensions")
if rcon.is_head_rank:
print("%d elements" % len(mesh.elements))
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
discr = rcon.make_discretization(mesh_data, order=5,
debug=[])
def dirichlet_bc(x, el):
from math import sin
return sin(10*x[0])
def rhs_c(x, el):
if la.norm(x) < 0.1:
return 1000
else:
return 0
def my_diff_tensor():
result = numpy.eye(dim)
result[0,0] = 0.1
return result
try:
from grudge.models.poisson import (
PoissonOperator,
HelmholtzOperator)
from grudge.second_order import \
IPDGSecondDerivative, LDGSecondDerivative, \
StabilizedCentralSecondDerivative
k = 1
from grudge.mesh import BTAG_NONE, BTAG_ALL
op = HelmholtzOperator(k, discr.dimensions,
#diffusion_tensor=my_diff_tensor(),
#dirichlet_tag="dirichlet",
#neumann_tag="neumann",
dirichlet_tag=BTAG_ALL,
neumann_tag=BTAG_NONE,
#dirichlet_tag=BTAG_ALL,
#neumann_tag=BTAG_NONE,
#dirichlet_bc=GivenFunction(dirichlet_bc),
dirichlet_bc=ConstantGivenFunction(0),
neumann_bc=ConstantGivenFunction(-10),
scheme=StabilizedCentralSecondDerivative(),
#scheme=LDGSecondDerivative(),
#scheme=IPDGSecondDerivative(),
)
bound_op = op.bind(discr)
if False:
from grudge.iterative import parallel_cg
u = -parallel_cg(rcon, -bound_op,
bound_op.prepare_rhs(discr.interpolate_volume_function(rhs_c)),
debug=20, tol=5e-4,
dot=discr.nodewise_dot_product,
x=discr.volume_zeros())
else:
rhs = bound_op.prepare_rhs(discr.interpolate_volume_function(rhs_c))
def compute_resid(x):
return bound_op(x)-rhs
from scipy.sparse.linalg import minres, LinearOperator
u, info = minres(
LinearOperator(
(len(discr), len(discr)),
matvec=bound_op, dtype=bound_op.dtype),
rhs,
callback=ResidualPrinter(compute_resid),
tol=1e-5)
print()
if info != 0:
raise RuntimeError("gmres reported error %d" % info)
print("finished gmres")
print(la.norm(bound_op(u)-rhs)/la.norm(rhs))
if write_output:
from grudge.visualization import SiloVisualizer, VtkVisualizer
vis = VtkVisualizer(discr, rcon)
visf = vis.make_file("fld")
vis.add_data(visf, [ ("sol", discr.convert_volume(u, kind="numpy")), ])
visf.close()
finally:
discr.close()
if __name__ == "__main__":
main()
examples/poisson/poisson.py deleted 100644 → 0
View file @ 90acbf2e
# 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 grudge.tools import Reflection, Rotation
def main(write_output=True):
from grudge.data import GivenFunction, ConstantGivenFunction
from grudge.backends import guess_run_context
rcon = guess_run_context()
dim = 2
def boundary_tagger(fvi, el, fn, points):
from math import atan2, pi
normal = el.face_normals[fn]
if -90/180*pi < atan2(normal[1], normal[0]) < 90/180*pi:
return ["neumann"]
else:
return ["dirichlet"]
if dim == 2:
if rcon.is_head_rank:
from grudge.mesh.generator import make_disk_mesh
mesh = make_disk_mesh(r=0.5, boundary_tagger=boundary_tagger,
max_area=1e-2)
elif dim == 3:
if rcon.is_head_rank:
from grudge.mesh.generator import make_ball_mesh
mesh = make_ball_mesh(max_volume=0.0001,
boundary_tagger=lambda fvi, el, fn, points:
["dirichlet"])
else:
raise RuntimeError("bad number of dimensions")
if rcon.is_head_rank:
print("%d elements" % len(mesh.elements))
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
discr = rcon.make_discretization(mesh_data, order=5,
debug=[])
def dirichlet_bc(x, el):
from math import sin
return sin(10*x[0])
def rhs_c(x, el):
if la.norm(x) < 0.1:
return 1000
else:
return 0
def my_diff_tensor():
result = numpy.eye(dim)
result[0,0] = 0.1
return result
try:
from grudge.models.poisson import PoissonOperator
from grudge.second_order import \
IPDGSecondDerivative, LDGSecondDerivative, \
StabilizedCentralSecondDerivative
from grudge.mesh import BTAG_NONE, BTAG_ALL
op = PoissonOperator(discr.dimensions,
diffusion_tensor=my_diff_tensor(),
#dirichlet_tag="dirichlet",
#neumann_tag="neumann",
dirichlet_tag=BTAG_ALL,
neumann_tag=BTAG_NONE,
#dirichlet_tag=BTAG_ALL,
#neumann_tag=BTAG_NONE,
dirichlet_bc=GivenFunction(dirichlet_bc),
neumann_bc=ConstantGivenFunction(-10),
scheme=StabilizedCentralSecondDerivative(),
#scheme=LDGSecondDerivative(),
#scheme=IPDGSecondDerivative(),
)
bound_op = op.bind(discr)
from grudge.iterative import parallel_cg
u = -parallel_cg(rcon, -bound_op,
bound_op.prepare_rhs(discr.interpolate_volume_function(rhs_c)),
debug=20, tol=5e-4,
dot=discr.nodewise_dot_product,
x=discr.volume_zeros())
if write_output:
from grudge.visualization import SiloVisualizer, VtkVisualizer
vis = VtkVisualizer(discr, rcon)
visf = vis.make_file("fld")
vis.add_data(visf, [ ("sol", discr.convert_volume(u, kind="numpy")), ])
visf.close()
finally:
discr.close()
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
from pytools.test import mark_test
@mark_test.long
def test_poisson():
main(write_output=False)
examples/wave/var-propagation-speed.py
View file @ 49b9b7c7
"""Variable-coefficient wave propagation."""
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
from six.moves import range
__copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
__copyright__ = """
Copyright (C) 2015 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
......@@ -28,172 +24,151 @@ THE SOFTWARE.
"""
import logging
import numpy as np
from grudge.mesh import BTAG_ALL, BTAG_NONE
def main(write_output=True,
dir_tag=BTAG_NONE,
neu_tag=BTAG_NONE,
rad_tag=BTAG_ALL,
flux_type_arg="upwind"):
from math import sin, cos, pi, exp, sqrt # noqa
from grudge.backends import guess_run_context
rcon = guess_run_context()
dim = 2
if dim == 1:
if rcon.is_head_rank:
from grudge.mesh.generator import make_uniform_1d_mesh
mesh = make_uniform_1d_mesh(-10, 10, 500)
elif dim == 2:
from grudge.mesh.generator import make_rect_mesh
if rcon.is_head_rank:
mesh = make_rect_mesh(a=(-1, -1), b=(1, 1), max_area=0.003)
elif dim == 3:
if rcon.is_head_rank:
from grudge.mesh.generator import make_ball_mesh
mesh = make_ball_mesh(max_volume=0.0005)
else:
raise RuntimeError("bad number of dimensions")
if rcon.is_head_rank:
print("%d elements" % len(mesh.elements))
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
discr = rcon.make_discretization(mesh_data, order=4)
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper()
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "fld")
source_center = np.array([0.7, 0.4])
source_width = 1/16
source_omega = 3
import grudge.symbolic as sym
sym_x = sym.nodes(2)
sym_source_center_dist = sym_x - source_center
from grudge.models.wave import VariableVelocityStrongWaveOperator
op = VariableVelocityStrongWaveOperator(
c=sym.If(sym.Comparison(
np.dot(sym_x, sym_x), "<", 0.4**2),
1, 0.5),
dimensions=discr.dimensions,
source=
sym.CFunction("sin")(source_omega*sym.ScalarParameter("t"))
* sym.CFunction("exp")(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2),
dirichlet_tag=dir_tag,
neumann_tag=neu_tag,
radiation_tag=rad_tag,
flux_type=flux_type_arg
import pyopencl as cl
import pyopencl.tools as cl_tools
from pytools.obj_array import flat_obj_array
from grudge import op
from grudge.array_context import PyOpenCLArrayContext
from grudge.discretization import make_discretization_collection
from grudge.shortcuts import set_up_rk4
logger = logging.getLogger(__name__)
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)
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dim,
b=(0.5,)*dim,
nelements_per_axis=(20,)*dim)
dcoll = make_discretization_collection(actx, mesh, order=order)
def source_f(actx, dcoll, t=0):
source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
source_width = 0.05
source_omega = 3
nodes = actx.thaw(dcoll.nodes())
source_center_dist = flat_obj_array(
[nodes[i] - source_center[i] for i in range(dcoll.dim)]
)
return (
np.sin(source_omega*t)
* actx.np.exp(
-np.dot(source_center_dist, source_center_dist)
/ source_width**2
)
)
from grudge.tools import join_fields
fields = join_fields(discr.volume_zeros(),
[discr.volume_zeros() for i in range(discr.dimensions)])
# {{{ diagnostics setup
from pytools.log import LogManager, \
add_general_quantities, \
add_simulation_quantities, \
add_run_info
if write_output:
log_file_name = "wave.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 pytools.log import IntervalTimer
vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
logmgr.add_quantity(vis_timer)
stepper.add_instrumentation(logmgr)
from grudge.log import LpNorm
u_getter = lambda: fields[0]
logmgr.add_quantity(LpNorm(u_getter, discr, 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.stability import \
approximate_rk4_relative_imag_stability_region
max_dt = (
1/discr.compile(op.max_eigenvalue_expr())()
* discr.dt_non_geometric_factor()
* discr.dt_geometric_factor()
* approximate_rk4_relative_imag_stability_region(stepper))
if flux_type_arg == "central":
max_dt *= 0.25
from grudge.timestep import times_and_steps
step_it = times_and_steps(final_time=3, logmgr=logmgr,
max_dt_getter=lambda t: max_dt)
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", fields[0]),
("v", fields[1:]),
],
time=t,
step=step)
visf.close()
x = actx.thaw(dcoll.nodes())
ones = dcoll.zeros(actx) + 1
c = actx.np.where(np.dot(x, x) < 0.15, 0.1 * ones, 0.2 * ones)
fields = stepper(fields, t, dt, rhs)
from meshmode.mesh import BTAG_ALL, BTAG_NONE
assert discr.norm(fields) < 1
finally:
if write_output:
vis.close()
from grudge.models.wave import VariableCoefficientWeakWaveOperator
logmgr.close()
discr.close()
wave_op = VariableCoefficientWeakWaveOperator(
dcoll,
actx.freeze(c),
source_f=source_f,
dirichlet_tag=BTAG_NONE,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_ALL,
flux_type="upwind"
)
# }}}
fields = flat_obj_array(
dcoll.zeros(actx),
[dcoll.zeros(actx) for i in range(dcoll.dim)]
)
if __name__ == "__main__":
main(flux_type_arg="upwind")
wave_op.check_bc_coverage(mesh)
def rhs(t, w):
return wave_op.operator(t, w)
dt = actx.to_numpy(
2/3 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
dt_stepper = set_up_rk4("w", dt, fields, rhs)
final_t = 1
nsteps = int(final_t/dt) + 1
logger.info("dt=%g nsteps=%d", dt, nsteps)
from grudge.shortcuts import make_visualizer
vis = make_visualizer(dcoll)
step = 0
def norm(u):
return op.norm(dcoll, u, 2)
# entry points for py.test ----------------------------------------------------
def test_var_velocity_wave():
from pytools.test import mark_test
mark_long = mark_test.long
from time import time
t_last_step = time()
for flux_type in ["upwind", "central"]:
yield ("dirichlet var-v wave equation with %s flux" % flux_type,
mark_long(main),
False, BTAG_ALL, BTAG_NONE, TAG_NONE, flux_type)
yield ("neumann var-v wave equation", mark_long(main),
False, BTAG_NONE, BTAG_ALL, TAG_NONE)
yield ("radiation-bc var-v wave equation", mark_long(main),
False, BTAG_NONE, TAG_NONE, BTAG_ALL)
if visualize:
u = fields[0]
v = fields[1:]
vis.write_vtk_file(
f"fld-var-propagation-speed-{step:04d}.vtu",
[
("u", u),
("v", v),
("c", c),
]
)
# vim: foldmethod=marker
for event in dt_stepper.run(t_end=final_t):
if isinstance(event, dt_stepper.StateComputed):
assert event.component_id == "w"
step += 1
if step % 10 == 0:
logger.info("step: %d t: %.8e L2: %.8e",
step, time() - t_last_step,
actx.to_numpy(norm(event.state_component[0])))
if visualize:
vis.write_vtk_file(
f"fld-var-propagation-speed-{step:04d}.vtu",
[
("u", event.state_component[0]),
("v", event.state_component[1:]),
("c", c),
]
)
t_last_step = time()
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert norm(u=event.state_component[0]) < 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/wave/wave-min.py → examples/wave/wave-min-mpi.py
View file @ 49b9b7c7
"""Minimal example of a grudge driver."""
from __future__ import division, print_function
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__copyright__ = """
Copyright (C) 2015 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
......@@ -25,101 +24,175 @@ THE SOFTWARE.
"""
import logging
import numpy as np
from mpi4py import MPI
import pyopencl as cl
import pyopencl.tools as cl_tools
from pytools.obj_array import flat_obj_array
import grudge.op as op
from grudge import make_discretization_collection
from grudge.array_context import MPIPyOpenCLArrayContext
from grudge.shortcuts import set_up_rk4
from grudge import sym, bind, Discretization
def main(write_output=True, order=4):
logger = logging.getLogger(__name__)
class WaveTag:
pass
def main(dim=2, order=4, visualize=True):
comm = MPI.COMM_WORLD
num_parts = comm.size
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)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
actx = MPIPyOpenCLArrayContext(comm, queue, allocator=allocator)
from meshmode.distributed import get_partition_by_pymetis, membership_list_to_map
from meshmode.mesh.processing import partition_mesh
if comm.rank == 0:
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dim,
b=(0.5,)*dim,
nelements_per_axis=(16,)*dim)
logger.info("%d elements", mesh.nelements)
part_id_to_part = partition_mesh(mesh,
membership_list_to_map(
get_partition_by_pymetis(mesh, num_parts)))
parts = [part_id_to_part[i] for i in range(num_parts)]
local_mesh = comm.scatter(parts)
del mesh
else:
local_mesh = comm.scatter(None)
dcoll = make_discretization_collection(actx, local_mesh, order=order)
def source_f(actx, dcoll, t=0):
source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
source_width = 0.05
source_omega = 3
nodes = actx.thaw(dcoll.nodes())
source_center_dist = flat_obj_array(
[nodes[i] - source_center[i] for i in range(dcoll.dim)]
)
return (
np.sin(source_omega*t)
* actx.np.exp(
-np.dot(source_center_dist, source_center_dist)
/ source_width**2
)
)
discr = Discretization(cl_ctx, mesh, order=order)
from meshmode.mesh import BTAG_ALL, BTAG_NONE
source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
source_width = 0.05
source_omega = 3
from grudge.models.wave import WeakWaveOperator
sym_x = sym.nodes(mesh.dim)
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
wave_op = WeakWaveOperator(
dcoll,
0.1,
source_f=source_f,
dirichlet_tag=BTAG_NONE,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_ALL,
flux_type="upwind"
)
from grudge.models.wave import StrongWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
op = StrongWaveOperator(-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
fields = flat_obj_array(
dcoll.zeros(actx),
[dcoll.zeros(actx) for i in range(dcoll.dim)]
)
def rhs(t, w):
return bound_op(queue, t=t, w=w)
dt = actx.to_numpy(
2/3 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
if mesh.dim == 2:
dt = 0.04
elif mesh.dim == 3:
dt = 0.02
wave_op.check_bc_coverage(local_mesh)
def rhs(t, w):
return wave_op.operator(t, w)
dt_stepper = set_up_rk4("w", dt, fields, rhs)
final_t = 10
nsteps = int(final_t/dt)
print("dt=%g nsteps=%d" % (dt, nsteps))
nsteps = int(final_t/dt) + 1
if comm.rank == 0:
logger.info("dt=%g nsteps=%d", dt, nsteps)
from grudge.shortcuts import make_visualizer
vis = make_visualizer(discr, vis_order=order)
vis = make_visualizer(dcoll)
step = 0
norm = bind(discr, sym.norm(2, sym.var("u")))
def norm(u):
return op.norm(dcoll, u, 2)
from time import time
t_last_step = time()
if visualize:
u = fields[0]
v = fields[1:]
vis.write_parallel_vtk_file(
comm,
f"fld-wave-min-mpi-{{rank:03d}}-{step:04d}.vtu",
[
("u", u),
("v", v),
]
)
for event in dt_stepper.run(t_end=final_t):
if isinstance(event, dt_stepper.StateComputed):
assert event.component_id == "w"
step += 1
l2norm = norm(u=event.state_component[0])
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,
if comm.rank == 0:
logger.info("step: %d t: %.8e L2: %.8e",
step, time() - t_last_step, l2norm)
if visualize:
vis.write_parallel_vtk_file(
comm,
f"fld-wave-min-mpi-{{rank:03d}}-{step:04d}.vtu",
[
("u", event.state_component[0]),
("v", event.state_component[1:]),
])
]
)
t_last_step = time()
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert l2norm < 1
if __name__ == "__main__":
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(
dim=args.dim,
order=args.order,
visualize=args.visualize)
examples/wave/wave-op-mpi.py 0 → 100644
View file @ 49b9b7c7
__copyright__ = """
Copyright (C) 2020 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
from dataclasses import dataclass
import numpy as np
import pyopencl as cl
import pyopencl.tools as cl_tools
from arraycontext import dataclass_array_container, with_container_arithmetic
from meshmode.dof_array import DOFArray
from meshmode.mesh import BTAG_ALL
from pytools.obj_array import flat_obj_array, make_obj_array
import grudge.geometry as geo
import grudge.op as op
from grudge.discretization import make_discretization_collection
from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, as_dofdesc
from grudge.shortcuts import compiled_lsrk45_step, make_visualizer
from grudge.trace_pair import TracePair
logger = logging.getLogger(__name__)
from mpi4py import MPI
# {{{ wave equation bits
@with_container_arithmetic(bcasts_across_obj_array=True,
rel_comparison=True,
_cls_has_array_context_attr=True,
)
@dataclass_array_container
@dataclass(frozen=True)
class WaveState:
u: DOFArray
v: np.ndarray # [object array]
# NOTE: disable numpy doing any array math
__array_ufunc__ = None
def __post_init__(self):
assert isinstance(self.v, np.ndarray) and self.v.dtype.char == "O"
@property
def array_context(self):
return self.u.array_context
def wave_flux(actx, dcoll, c, w_tpair):
u = w_tpair.u
v = w_tpair.v
dd = w_tpair.dd
normal = geo.normal(actx, dcoll, dd)
flux_weak = WaveState(
u=v.avg @ normal,
v=u.avg * normal
)
# upwind
v_jump = v.diff @ normal
flux_weak += WaveState(
u=0.5 * u.diff,
v=0.5 * v_jump * normal,
)
return op.project(dcoll, dd, dd.with_domain_tag("all_faces"), c*flux_weak)
class _WaveStateTag:
pass
def wave_operator(actx, dcoll, c, w, quad_tag=None):
dd_base = as_dofdesc("vol")
dd_vol = as_dofdesc("vol", quad_tag)
dd_faces = as_dofdesc("all_faces", quad_tag)
dd_btag = as_dofdesc(BTAG_ALL).with_discr_tag(quad_tag)
def interp_to_surf_quad(utpair):
local_dd = utpair.dd
local_dd_quad = local_dd.with_discr_tag(quad_tag)
return TracePair(
local_dd_quad,
interior=op.project(dcoll, local_dd, local_dd_quad, utpair.int),
exterior=op.project(dcoll, local_dd, local_dd_quad, utpair.ext)
)
w_quad = op.project(dcoll, dd_base, dd_vol, w)
u = w_quad.u
v = w_quad.v
dir_w = op.project(dcoll, dd_base, dd_btag, w)
dir_u = dir_w.u
dir_v = dir_w.v
dir_bval = WaveState(u=dir_u, v=dir_v)
dir_bc = WaveState(u=-dir_u, v=dir_v)
return (
op.inverse_mass(
dcoll,
WaveState(
u=-c*op.weak_local_div(dcoll, dd_vol, v),
v=-c*op.weak_local_grad(dcoll, dd_vol, u)
)
+ op.face_mass(
dcoll,
dd_faces,
wave_flux(
actx,
dcoll, c=c,
w_tpair=op.bdry_trace_pair(dcoll,
dd_btag,
interior=dir_bval,
exterior=dir_bc)
) + sum(
wave_flux(actx, dcoll, c=c, w_tpair=interp_to_surf_quad(tpair))
for tpair in op.interior_trace_pairs(dcoll, w,
comm_tag=_WaveStateTag)
)
)
)
)
# }}}
def estimate_rk4_timestep(actx, dcoll, c):
from grudge.dt_utils import characteristic_lengthscales
local_dts = characteristic_lengthscales(actx, dcoll) / c
return op.nodal_min(dcoll, "vol", local_dts)
def bump(actx, dcoll, t=0):
source_center = np.array([0.2, 0.35, 0.1])[:dcoll.dim]
source_width = 0.05
source_omega = 3
nodes = actx.thaw(dcoll.nodes())
center_dist = flat_obj_array([
nodes[i] - source_center[i]
for i in range(dcoll.dim)
])
return (
np.cos(source_omega*t)
* actx.np.exp(
-np.dot(center_dist, center_dist)
/ source_width**2))
def main(ctx_factory, dim=2, order=3,
visualize=False, lazy=False, numpy=False, use_quad=False,
use_nonaffine_mesh=False, no_diagnostics=False):
comm = MPI.COMM_WORLD
num_parts = comm.size
from grudge.array_context import get_reasonable_array_context_class
actx_class = get_reasonable_array_context_class(lazy=lazy,
distributed=True, numpy=numpy)
if numpy:
actx = actx_class(comm)
else:
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
if lazy:
actx = actx_class(comm, queue, allocator=allocator, mpi_base_tag=15000)
else:
actx = actx_class(comm, queue, allocator=allocator)
from meshmode.distributed import get_partition_by_pymetis, membership_list_to_map
from meshmode.mesh.processing import partition_mesh
nel_1d = 16
if comm.rank == 0:
if use_nonaffine_mesh:
from meshmode.mesh.generation import generate_warped_rect_mesh
# FIXME: *generate_warped_rect_mesh* in meshmode warps a
# rectangle domain with hard-coded vertices at (-0.5,)*dim
# and (0.5,)*dim. Should extend the function interface to allow
# for specifying the end points directly.
mesh = generate_warped_rect_mesh(dim=dim, order=order,
nelements_side=nel_1d)
else:
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dim,
b=(0.5,)*dim,
nelements_per_axis=(nel_1d,)*dim)
logger.info("%d elements", mesh.nelements)
part_id_to_part = partition_mesh(mesh,
membership_list_to_map(
get_partition_by_pymetis(mesh, num_parts)))
parts = [part_id_to_part[i] for i in range(num_parts)]
local_mesh = comm.scatter(parts)
del mesh
else:
local_mesh = comm.scatter(None)
from meshmode.discretization.poly_element import (
QuadratureSimplexGroupFactory,
default_simplex_group_factory,
)
dcoll = make_discretization_collection(
actx, local_mesh,
discr_tag_to_group_factory={
DISCR_TAG_BASE: default_simplex_group_factory(base_dim=dim, order=order),
# High-order quadrature to integrate inner products of polynomials
# on warped geometry exactly (metric terms are also polynomial)
DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(3*order),
})
if use_quad:
quad_tag = DISCR_TAG_QUAD
else:
quad_tag = None
fields = WaveState(
u=bump(actx, dcoll),
v=make_obj_array([dcoll.zeros(actx) for i in range(dcoll.dim)])
)
c = 1
# FIXME: Sketchy, empirically determined fudge factor
# 5/4 to account for larger LSRK45 stability region
dt = actx.to_numpy(0.45 * estimate_rk4_timestep(actx, dcoll, c)) * 5/4
vis = make_visualizer(dcoll)
def rhs(t, w):
return wave_operator(actx, dcoll, c=c, w=w, quad_tag=quad_tag)
compiled_rhs = actx.compile(rhs)
if comm.rank == 0:
logger.info("dt = %g", dt)
import time
start = time.time()
t = 0
t_final = 3
istep = 0
while t < t_final:
start = time.time()
fields = compiled_lsrk45_step(actx, fields, t, dt, compiled_rhs)
if istep % 10 == 0:
stop = time.time()
if no_diagnostics:
if comm.rank == 0:
logger.info("step: %d t: %.8e wall: %.8es",
istep, t, stop - start)
else:
l2norm = actx.to_numpy(op.norm(dcoll, fields.u, 2))
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert l2norm < 1
linfnorm = actx.to_numpy(op.norm(dcoll, fields.u, np.inf))
nodalmax = actx.to_numpy(op.nodal_max(dcoll, "vol", fields.u))
nodalmin = actx.to_numpy(op.nodal_min(dcoll, "vol", fields.u))
if comm.rank == 0:
logger.info("step: %d t: %.8e L2: %.8e Linf: %.8e "
"sol max: %.8e sol min: %.8e wall: %.8e",
istep, t, l2norm, linfnorm, nodalmax, nodalmin,
stop - start)
if visualize:
vis.write_parallel_vtk_file(
comm,
f"fld-wave-eager-mpi-{{rank:03d}}-{istep:04d}.vtu",
[
("u", fields.u),
("v", fields.v),
]
)
start = stop
t += dt
istep += 1
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--dim", default=2, type=int)
parser.add_argument("--order", default=3, type=int)
parser.add_argument("--visualize", action="store_true")
parser.add_argument("--lazy", action="store_true",
help="switch to a lazy computation mode")
parser.add_argument("--numpy", action="store_true",
help="switch to numpy-based array context")
parser.add_argument("--quad", action="store_true")
parser.add_argument("--nonaffine", action="store_true")
parser.add_argument("--no-diagnostics", 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,
lazy=args.lazy,
numpy=args.numpy,
use_quad=args.quad,
use_nonaffine_mesh=args.nonaffine,
no_diagnostics=args.no_diagnostics)
# vim: foldmethod=marker
examples/wave/wave-op-var-velocity.py 0 → 100644
View file @ 49b9b7c7
__copyright__ = """
Copyright (C) 2020 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 numpy as np
import numpy.linalg as la # noqa
import pyopencl as cl
import pyopencl.tools as cl_tools
from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
from pytools.obj_array import flat_obj_array
import grudge.geometry as geo
import grudge.op as op
from grudge.array_context import PyOpenCLArrayContext
from grudge.discretization import make_discretization_collection
from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, as_dofdesc
from grudge.shortcuts import make_visualizer, rk4_step
logger = logging.getLogger(__name__)
# {{{ wave equation bits
def wave_flux(actx, dcoll, c, w_tpair):
dd = w_tpair.dd
dd_quad = dd.with_discr_tag(DISCR_TAG_QUAD)
u = w_tpair[0]
v = w_tpair[1:]
normal = geo.normal(actx, dcoll, dd)
flux_weak = flat_obj_array(
np.dot(v.avg, normal),
normal*u.avg,
)
# upwind
flux_weak += flat_obj_array(
0.5*(u.ext-u.int),
0.5*normal*np.dot(normal, v.ext-v.int),
)
# FIXME this flux is only correct for continuous c
dd_allfaces_quad = dd_quad.with_domain_tag("all_faces")
c_quad = op.project(dcoll, "vol", dd_quad, c)
flux_quad = op.project(dcoll, dd, dd_quad, flux_weak)
return op.project(dcoll, dd_quad, dd_allfaces_quad, c_quad*flux_quad)
def wave_operator(actx, dcoll, c, w):
u = w[0]
v = w[1:]
dir_u = op.project(dcoll, "vol", BTAG_ALL, u)
dir_v = op.project(dcoll, "vol", BTAG_ALL, v)
dir_bval = flat_obj_array(dir_u, dir_v)
dir_bc = flat_obj_array(-dir_u, dir_v)
dd_quad = as_dofdesc("vol", DISCR_TAG_QUAD)
c_quad = op.project(dcoll, "vol", dd_quad, c)
w_quad = op.project(dcoll, "vol", dd_quad, w)
u_quad = w_quad[0]
v_quad = w_quad[1:]
dd_allfaces_quad = as_dofdesc("all_faces", DISCR_TAG_QUAD)
return (
op.inverse_mass(
dcoll,
flat_obj_array(
-op.weak_local_div(dcoll, dd_quad, c_quad*v_quad),
-op.weak_local_grad(dcoll, dd_quad, c_quad*u_quad)
# pylint: disable=invalid-unary-operand-type
) + op.face_mass(
dcoll,
dd_allfaces_quad,
wave_flux(
actx,
dcoll, c=c,
w_tpair=op.bdry_trace_pair(dcoll,
BTAG_ALL,
interior=dir_bval,
exterior=dir_bc)
) + sum(
wave_flux(actx, dcoll, c=c, w_tpair=tpair)
for tpair in op.interior_trace_pairs(dcoll, w)
)
)
)
)
# }}}
def estimate_rk4_timestep(actx, dcoll, c):
from grudge.dt_utils import characteristic_lengthscales
local_dts = characteristic_lengthscales(actx, dcoll) / c
return op.nodal_min(dcoll, "vol", local_dts)
def bump(actx, dcoll, t=0, width=0.05, center=None):
if center is None:
center = np.array([0.2, 0.35, 0.1])
center = center[:dcoll.dim]
source_omega = 3
nodes = actx.thaw(dcoll.nodes())
center_dist = flat_obj_array([
nodes[i] - center[i]
for i in range(dcoll.dim)
])
return (
np.cos(source_omega*t)
* actx.np.exp(
-np.dot(center_dist, center_dist)
/ width**2))
def main(ctx_factory, dim=2, order=3, visualize=False):
cl_ctx = ctx_factory()
queue = cl.CommandQueue(cl_ctx)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
actx = PyOpenCLArrayContext(queue, allocator=allocator)
nel_1d = 16
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dim,
b=(0.5,)*dim,
nelements_per_axis=(nel_1d,)*dim)
logger.info("%d elements", mesh.nelements)
from meshmode.discretization.poly_element import (
QuadratureSimplexGroupFactory,
default_simplex_group_factory,
)
dcoll = make_discretization_collection(
actx, mesh,
discr_tag_to_group_factory={
DISCR_TAG_BASE: default_simplex_group_factory(base_dim=dim, order=order),
DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(3*order),
}
)
# bounded above by 1
c = 0.2 + 0.8*bump(actx, dcoll, center=np.zeros(3), width=0.5)
dt = actx.to_numpy(0.5 * estimate_rk4_timestep(actx, dcoll, c=1))
fields = flat_obj_array(
bump(actx, dcoll, ),
[dcoll.zeros(actx) for i in range(dcoll.dim)]
)
vis = make_visualizer(dcoll)
def rhs(t, w):
return wave_operator(actx, dcoll, c=c, w=w)
logger.info("dt = %g", dt)
t = 0
t_final = 3
istep = 0
while t < t_final:
fields = rk4_step(fields, t, dt, rhs)
l2norm = actx.to_numpy(op.norm(dcoll, fields[0], 2))
if istep % 10 == 0:
linfnorm = actx.to_numpy(op.norm(dcoll, fields[0], np.inf))
nodalmax = actx.to_numpy(op.nodal_max(dcoll, "vol", fields[0]))
nodalmin = actx.to_numpy(op.nodal_min(dcoll, "vol", fields[0]))
logger.info("step: %d t: %.8e L2: %.8e Linf: %.8e "
"sol max: %.8e sol min: %.8e",
istep, t, l2norm, linfnorm, nodalmax, nodalmin)
if visualize:
vis.write_vtk_file(
f"fld-wave-eager-var-velocity-{istep:04d}.vtu",
[
("c", c),
("u", fields[0]),
("v", fields[1:]),
]
)
t += dt
istep += 1
# NOTE: These are here to ensure the solution is bounded for the
# time interval specified
assert l2norm < 1
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--dim", default=2, type=int)
parser.add_argument("--order", default=3, 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)
# vim: foldmethod=marker
examples/wave/wave.py deleted 100644 → 0
View file @ 90acbf2e
# 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 as np
from grudge.mesh import BTAG_ALL, BTAG_NONE
from six.moves import range
def main(write_output=True,
dir_tag=BTAG_NONE, neu_tag=TAG_NONE, rad_tag=BTAG_ALL,
flux_type_arg="upwind", dtype=np.float64, debug=[]):
from math import sin, cos, pi, exp, sqrt # noqa
from grudge.backends import guess_run_context
rcon = guess_run_context()
dim = 2
if dim == 1:
if rcon.is_head_rank:
from grudge.mesh.generator import make_uniform_1d_mesh
mesh = make_uniform_1d_mesh(-10, 10, 500)
elif dim == 2:
from grudge.mesh.generator import make_rect_mesh
if rcon.is_head_rank:
mesh = make_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), max_area=0.008)
elif dim == 3:
if rcon.is_head_rank:
from grudge.mesh.generator import make_ball_mesh
mesh = make_ball_mesh(max_volume=0.0005)
else:
raise RuntimeError("bad number of dimensions")
if rcon.is_head_rank:
print("%d elements" % len(mesh.elements))
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper(dtype=dtype)
from grudge.models.wave import StrongWaveOperator
from grudge.mesh import BTAG_ALL, BTAG_NONE # noqa
source_center = np.array([0.1, 0.22])
source_width = 0.05
source_omega = 3
import grudge.symbolic as sym
sym_x = sym.nodes(2)
sym_source_center_dist = sym_x - source_center
op = StrongWaveOperator(-1, dim,
source_f=
sym.CFunction("sin")(source_omega*sym.ScalarParameter("t"))
* sym.CFunction("exp")(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2),
dirichlet_tag=dir_tag,
neumann_tag=neu_tag,
radiation_tag=rad_tag,
flux_type=flux_type_arg
)
discr = rcon.make_discretization(mesh_data, order=4, debug=debug,
default_scalar_type=dtype,
tune_for=op.sym_operator())
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "fld")
from grudge.tools import join_fields
fields = join_fields(discr.volume_zeros(dtype=dtype),
[discr.volume_zeros(dtype=dtype) for i in range(discr.dimensions)])
# {{{ diagnostics setup
from pytools.log import LogManager, \
add_general_quantities, \
add_simulation_quantities, \
add_run_info
if write_output:
log_file_name = "wave.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 pytools.log import IntervalTimer
vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
logmgr.add_quantity(vis_timer)
stepper.add_instrumentation(logmgr)
from grudge.log import LpNorm
u_getter = lambda: fields[0]
logmgr.add_quantity(LpNorm(u_getter, discr, 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=4, 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 % 10 == 0 and write_output:
visf = vis.make_file("fld-%04d" % step)
vis.add_data(visf,
[
("u", discr.convert_volume(fields[0], kind="numpy")),
("v", discr.convert_volume(fields[1:], kind="numpy")),
],
time=t,
step=step)
visf.close()
fields = stepper(fields, t, dt, rhs)
assert discr.norm(fields) < 1
assert fields[0].dtype == dtype
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
# }}}
if __name__ == "__main__":
main(True, BTAG_ALL, BTAG_NONE, TAG_NONE, "upwind", np.float64,
debug=["cuda_no_plan", "dump_optemplate_stages"])
# {{{ entry points for py.test
def test_wave():
from pytools.test import mark_test
mark_long = mark_test.long
yield ("dirichlet wave equation with SP data", mark_long(main),
False, BTAG_ALL, BTAG_NONE, TAG_NONE, "upwind", np.float64)
yield ("dirichlet wave equation with SP complex data", mark_long(main),
False, BTAG_ALL, BTAG_NONE, TAG_NONE, "upwind", np.complex64)
yield ("dirichlet wave equation with DP complex data", mark_long(main),
False, BTAG_ALL, BTAG_NONE, TAG_NONE, "upwind", np.complex128)
for flux_type in ["upwind", "central"]:
yield ("dirichlet wave equation with %s flux" % flux_type,
mark_long(main),
False, BTAG_ALL, BTAG_NONE, TAG_NONE, flux_type)
yield ("neumann wave equation", mark_long(main),
False, BTAG_NONE, BTAG_ALL, TAG_NONE)
yield ("radiation-bc wave equation", mark_long(main),
False, BTAG_NONE, TAG_NONE, BTAG_ALL)
# }}}
# ij
examples/wave/wiggly.py deleted 100644 → 0
View file @ 90acbf2e
"""Wiggly geometry wave propagation."""
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
from six.moves import range
__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.mesh import BTAG_ALL, BTAG_NONE # noqa
def main(write_output=True,
flux_type_arg="upwind", dtype=np.float64, debug=[]):
from math import sin, cos, pi, exp, sqrt # noqa
from grudge.backends import guess_run_context
rcon = guess_run_context()
if rcon.is_head_rank:
from grudge.mesh.reader.gmsh import generate_gmsh
mesh = generate_gmsh(GEOMETRY, 2,
allow_internal_boundaries=True,
force_dimension=2)
print("%d elements" % len(mesh.elements))
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
discr = rcon.make_discretization(mesh_data, order=4, debug=debug,
default_scalar_type=dtype)
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper(dtype=dtype)
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "fld")
source_center = 0
source_width = 0.05
source_omega = 3
import grudge.symbolic as sym
sym_x = sym.nodes(2)
sym_source_center_dist = sym_x - source_center
from grudge.models.wave import StrongWaveOperator
op = StrongWaveOperator(-1, discr.dimensions,
source_f=
sym.CFunction("sin")(source_omega*sym.ScalarParameter("t"))
* sym.CFunction("exp")(
-np.dot(sym_source_center_dist, sym_source_center_dist)
/ source_width**2),
dirichlet_tag="boundary",
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_NONE,
flux_type=flux_type_arg
)
from grudge.tools import join_fields
fields = join_fields(discr.volume_zeros(dtype=dtype),
[discr.volume_zeros(dtype=dtype) for i in range(discr.dimensions)])
# diagnostics setup -------------------------------------------------------
from pytools.log import LogManager, \
add_general_quantities, \
add_simulation_quantities, \
add_run_info
if write_output:
log_file_name = "wiggly.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)
logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
# timestep loop -----------------------------------------------------------
rhs = op.bind(discr)
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=4, 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 % 10 == 0 and write_output:
visf = vis.make_file("fld-%04d" % step)
vis.add_data(visf,
[
("u", fields[0]),
("v", fields[1:]),
],
time=t,
step=step)
visf.close()
fields = stepper(fields, t, dt, rhs)
assert discr.norm(fields) < 1
assert fields[0].dtype == dtype
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
GEOMETRY = """
w = 1;
dx = 0.2;
ch_width = 0.2;
rows = 4;
Point(0) = {0,0,0};
Point(1) = {w,0,0};
bottom_line = newl;
Line(bottom_line) = {0,1};
left_pts[] = { 0 };
right_pts[] = { 1 };
left_pts[] = { };
emb_lines[] = {};
For row In {1:rows}
If (row % 2 == 0)
// left
rp = newp; Point(rp) = {w,dx*row, 0};
right_pts[] += {rp};
mp = newp; Point(mp) = {ch_width,dx*row, 0};
emb_line = newl; Line(emb_line) = {mp,rp};
emb_lines[] += {emb_line};
EndIf
If (row % 2)
// right
lp = newp; Point(lp) = {0,dx*row, 0};
left_pts[] += {lp};
mp = newp; Point(mp) = { w-ch_width,dx*row, 0};
emb_line = newl; Line(emb_line) = {mp,lp};
emb_lines[] += {emb_line};
EndIf
EndFor
lep = newp; Point(lep) = {0,(rows+1)*dx,0};
rep = newp; Point(rep) = {w,(rows+1)*dx,0};
top_line = newl; Line(top_line) = {lep,rep};
left_pts[] += { lep };
right_pts[] += { rep };
lines[] = {bottom_line};
For i In {0:#right_pts[]-2}
l = newl; Line(l) = {right_pts[i], right_pts[i+1]};
lines[] += {l};
EndFor
lines[] += {-top_line};
For i In {#left_pts[]-1:0:-1}
l = newl; Line(l) = {left_pts[i], left_pts[i-1]};
lines[] += {l};
EndFor
Line Loop (1) = lines[];
Plane Surface (1) = {1};
Physical Surface(1) = {1};
For i In {0:#emb_lines[]-1}
Line { emb_lines[i] } In Surface { 1 };
EndFor
boundary_lines[] = {};
boundary_lines[] += lines[];
boundary_lines[] += emb_lines[];
Physical Line ("boundary") = boundary_lines[];
Mesh.CharacteristicLengthFactor = 0.4;
"""
if __name__ == "__main__":
main()
grudge/__init__.py
View file @ 49b9b7c7
from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__license__ = """
......@@ -22,6 +20,12 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import grudge.sym # noqa
from grudge.execution import bind # noqa
from grudge.discretization import Discretization # noqa
from grudge.discretization import (
DiscretizationCollection,
make_discretization_collection,
)
__all__ = [
"DiscretizationCollection", "make_discretization_collection"
]
grudge/array_context.py 0 → 100644
View file @ 49b9b7c7
"""
.. autoclass:: PyOpenCLArrayContext
.. autoclass:: PytatoPyOpenCLArrayContext
.. autoclass:: MPIBasedArrayContext
.. autoclass:: MPIPyOpenCLArrayContext
.. autoclass:: MPINumpyArrayContext
.. class:: MPIPytatoArrayContext
.. autofunction:: get_reasonable_array_context_class
"""
from __future__ import annotations
__copyright__ = "Copyright (C) 2020 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.
"""
# {{{ imports
import logging
from collections.abc import Callable, Mapping
from dataclasses import dataclass
from typing import TYPE_CHECKING, Any
from warnings import warn
from typing_extensions import Self
from meshmode.array_context import (
PyOpenCLArrayContext as _PyOpenCLArrayContextBase,
PytatoPyOpenCLArrayContext as _PytatoPyOpenCLArrayContextBase,
)
from pytools import to_identifier
from pytools.tag import Tag
logger = logging.getLogger(__name__)
try:
# FIXME: temporary workaround while FusionContractorArrayContext
# is not available in meshmode's main branch
from meshmode.array_context import FusionContractorArrayContext
try:
# Crude check if we have the correct loopy branch
# (https://github.com/kaushikcfd/loopy/tree/pytato-array-context-transforms)
from loopy.transform.loop_fusion import get_kennedy_unweighted_fusion_candidates # noqa
except ImportError:
warn("Your loopy and meshmode branches are mismatched. "
"Please make sure that you have the "
"https://github.com/kaushikcfd/loopy/tree/pytato-array-context-transforms "
"branch of loopy.", stacklevel=1)
_HAVE_FUSION_ACTX = False
else:
_HAVE_FUSION_ACTX = True
except ImportError:
_HAVE_FUSION_ACTX = False
from arraycontext import ArrayContext, NumpyArrayContext
from arraycontext.container import ArrayContainer
from arraycontext.impl.pytato.compile import LazilyPyOpenCLCompilingFunctionCaller
from arraycontext.pytest import (
_PytestNumpyArrayContextFactory,
_PytestPyOpenCLArrayContextFactoryWithClass,
_PytestPytatoPyOpenCLArrayContextFactory,
register_pytest_array_context_factory,
)
if TYPE_CHECKING:
import pytato as pt
from mpi4py import MPI
from pytato import DistributedGraphPartition
from pytato.partition import PartId
import pyopencl
import pyopencl.tools
class PyOpenCLArrayContext(_PyOpenCLArrayContextBase):
"""Inherits from :class:`meshmode.array_context.PyOpenCLArrayContext`. Extends it
to understand :mod:`grudge`-specific transform metadata. (Of which there isn't
any, for now.)
"""
def __init__(self, queue: pyopencl.CommandQueue,
allocator: pyopencl.tools.AllocatorBase | None = None,
wait_event_queue_length: int | None = None,
force_device_scalars: bool | None = None) -> None:
if allocator is None:
warn("No memory allocator specified, please pass one. "
"(Preferably a pyopencl.tools.MemoryPool in order "
"to reduce device allocations)", stacklevel=2)
super().__init__(queue, allocator,
wait_event_queue_length, force_device_scalars)
# }}}
# {{{ pytato
class PytatoPyOpenCLArrayContext(_PytatoPyOpenCLArrayContextBase):
"""Inherits from :class:`meshmode.array_context.PytatoPyOpenCLArrayContext`.
Extends it to understand :mod:`grudge`-specific transform metadata. (Of
which there isn't any, for now.)
"""
def __init__(self, queue, allocator=None,
*,
compile_trace_callback: Callable[[Any, str, Any], None] | None
= None) -> None:
"""
:arg compile_trace_callback: A function of three arguments
*(what, stage, ir)*, where *what* identifies the object
being compiled, *stage* is a string describing the compilation
pass, and *ir* is an object containing the intermediate
representation. This interface should be considered
unstable.
"""
if allocator is None:
warn("No memory allocator specified, please pass one. "
"(Preferably a pyopencl.tools.MemoryPool in order "
"to reduce device allocations)", stacklevel=2)
super().__init__(queue, allocator,
compile_trace_callback=compile_trace_callback)
# }}}
class MPIBasedArrayContext(ArrayContext):
mpi_communicator: MPI.Intracomm
# {{{ distributed + pytato
@dataclass(frozen=True)
class _DistributedPartProgramID:
f: Callable[..., Any]
part_id: Any
def __str__(self):
name = getattr(self.f, "__name__", "anonymous")
if not name.isidentifier():
name = to_identifier(name)
part = to_identifier(str(self.part_id))
if part:
return f"{name}_part{part}"
else:
return name
class _DistributedLazilyPyOpenCLCompilingFunctionCaller(
LazilyPyOpenCLCompilingFunctionCaller):
def _dag_to_compiled_func(self, dict_of_named_arrays,
input_id_to_name_in_program, output_id_to_name_in_program,
output_template):
import pytato as pt
from pytools import ProcessLogger
self.actx._compile_trace_callback(self.f, "pre_deduplicate_data_wrappers",
dict_of_named_arrays)
with ProcessLogger(logger, "deduplicate_data_wrappers[pre-partition]"):
dict_of_named_arrays = pt.transform.deduplicate_data_wrappers(
dict_of_named_arrays)
self.actx._compile_trace_callback(self.f, "post_deduplicate_data_wrappers",
dict_of_named_arrays)
self.actx._compile_trace_callback(self.f, "pre_materialize",
dict_of_named_arrays)
with ProcessLogger(logger, "materialize_with_mpms[pre-partition]"):
dict_of_named_arrays = pt.transform.materialize_with_mpms(
dict_of_named_arrays)
self.actx._compile_trace_callback(self.f, "post_materialize",
dict_of_named_arrays)
# FIXME: Remove the import failure handling once this is in upstream grudge
self.actx._compile_trace_callback(self.f, "pre_infer_axes_tags",
dict_of_named_arrays)
with ProcessLogger(logger,
"transform_dag.infer_axes_tags[pre-partition]"):
from meshmode.transform_metadata import DiscretizationEntityAxisTag
dict_of_named_arrays = pt.unify_axes_tags(
dict_of_named_arrays,
tag_t=DiscretizationEntityAxisTag,
)
self.actx._compile_trace_callback(self.f, "post_infer_axes_tags",
dict_of_named_arrays)
self.actx._compile_trace_callback(self.f, "pre_find_distributed_partition",
dict_of_named_arrays)
distributed_partition = pt.find_distributed_partition(
# pylint-ignore-reason:
# '_BasePytatoArrayContext' has no
# 'mpi_communicator' member
# pylint: disable=no-member
self.actx.mpi_communicator, dict_of_named_arrays)
if __debug__:
# pylint-ignore-reason:
# '_BasePytatoArrayContext' has no 'mpi_communicator' member
pt.verify_distributed_partition(
self.actx.mpi_communicator, # pylint: disable=no-member
distributed_partition)
self.actx._compile_trace_callback(self.f, "post_find_distributed_partition",
distributed_partition)
# {{{ turn symbolic tags into globally agreed-upon integers
self.actx._compile_trace_callback(self.f, "pre_number_distributed_tags",
distributed_partition)
from pytato import number_distributed_tags
prev_mpi_base_tag = self.actx.mpi_base_tag
# type-ignore-reason: 'PytatoPyOpenCLArrayContext' has no 'mpi_communicator'
# pylint: disable=no-member
distributed_partition, new_mpi_base_tag = number_distributed_tags(
self.actx.mpi_communicator,
distributed_partition,
base_tag=prev_mpi_base_tag)
assert prev_mpi_base_tag == self.actx.mpi_base_tag
# FIXME: Updating stuff inside the array context from here is *cough*
# not super pretty.
self.actx.mpi_base_tag = new_mpi_base_tag
self.actx._compile_trace_callback(self.f, "post_number_distributed_tags",
distributed_partition)
# }}}
part_id_to_prg = {}
name_in_program_to_tags = {}
name_in_program_to_axes = {}
from pytato import make_dict_of_named_arrays
for part in distributed_partition.parts.values():
d = make_dict_of_named_arrays(
{var_name: distributed_partition.name_to_output[var_name]
for var_name in part.output_names
})
(
part_id_to_prg[part.pid],
part_prg_name_to_tags,
part_prg_name_to_axes
) = self._dag_to_transformed_pytato_prg(
d, prg_id=_DistributedPartProgramID(self.f, part.pid))
assert not (set(name_in_program_to_tags.keys())
& set(part_prg_name_to_tags.keys()))
assert not (set(name_in_program_to_axes.keys())
& set(part_prg_name_to_axes.keys()))
name_in_program_to_tags.update(part_prg_name_to_tags)
name_in_program_to_axes.update(part_prg_name_to_axes)
from constantdict import constantdict
return _DistributedCompiledFunction(
actx=self.actx,
distributed_partition=distributed_partition,
part_id_to_prg=part_id_to_prg,
input_id_to_name_in_program=input_id_to_name_in_program,
output_id_to_name_in_program=output_id_to_name_in_program,
name_in_program_to_tags=constantdict(name_in_program_to_tags),
name_in_program_to_axes=constantdict(name_in_program_to_axes),
output_template=output_template)
@dataclass(frozen=True)
class _DistributedCompiledFunction:
"""
A callable which captures the :class:`pytato.target.BoundProgram` resulting
from calling :attr:`~LazilyPyOpenCLCompilingFunctionCaller.f` with a given
set of input types, and generating :mod:`loopy` IR from it.
.. attribute:: pytato_program
.. attribute:: input_id_to_name_in_program
A mapping from input id to the placeholder name in
:attr:`CompiledFunction.pytato_program`. Input id is represented as the
position of :attr:`~LazilyPyOpenCLCompilingFunctionCaller.f`'s argument
augmented with the leaf array's key if the argument is an array
container.
.. attribute:: output_id_to_name_in_program
A mapping from output id to the name of
:class:`pytato.array.NamedArray` in
:attr:`CompiledFunction.pytato_program`. Output id is represented by
the key of a leaf array in the array container
:attr:`CompiledFunction.output_template`.
.. attribute:: output_template
An instance of :class:`arraycontext.ArrayContainer` that is the return
type of the callable.
"""
actx: MPIBasedArrayContext
distributed_partition: DistributedGraphPartition
part_id_to_prg: Mapping[PartId, pt.target.BoundProgram]
input_id_to_name_in_program: Mapping[tuple[Any, ...], str]
output_id_to_name_in_program: Mapping[tuple[Any, ...], str]
name_in_program_to_tags: Mapping[str, frozenset[Tag]]
name_in_program_to_axes: Mapping[str, tuple[pt.Axis, ...]]
output_template: ArrayContainer
def __call__(self, arg_id_to_arg) -> ArrayContainer:
"""
:arg arg_id_to_arg: Mapping from input id to the passed argument. See
:attr:`CompiledFunction.input_id_to_name_in_program` for input id's
representation.
"""
from arraycontext.impl.pyopencl.taggable_cl_array import to_tagged_cl_array
from arraycontext.impl.pytato.compile import _args_to_device_buffers
from arraycontext.impl.pytato.utils import get_cl_axes_from_pt_axes
input_args_for_prg = _args_to_device_buffers(
self.actx, self.input_id_to_name_in_program, arg_id_to_arg)
from pytato import execute_distributed_partition
assert isinstance(self.actx, PytatoPyOpenCLArrayContext | PyOpenCLArrayContext)
out_dict = execute_distributed_partition(
self.distributed_partition, self.part_id_to_prg,
self.actx.queue, self.actx.mpi_communicator, # pylint: disable=no-member
allocator=self.actx.allocator, # pylint: disable=no-member
input_args=input_args_for_prg)
def to_output_template(keys, _):
ary_name_in_prg = self.output_id_to_name_in_program[keys]
return self.actx.thaw(to_tagged_cl_array(
out_dict[ary_name_in_prg],
axes=get_cl_axes_from_pt_axes(
self.name_in_program_to_axes[ary_name_in_prg]),
tags=self.name_in_program_to_tags[ary_name_in_prg]))
from arraycontext.container.traversal import rec_keyed_map_array_container
return rec_keyed_map_array_container(to_output_template,
self.output_template)
# }}}
# {{{ distributed + pyopencl
class MPIPyOpenCLArrayContext(PyOpenCLArrayContext, MPIBasedArrayContext):
"""An array context for using distributed computation with :mod:`pyopencl`
eager evaluation.
.. autofunction:: __init__
"""
def __init__(self,
mpi_communicator,
queue: pyopencl.CommandQueue,
*, allocator: pyopencl.tools.AllocatorBase | None = None,
wait_event_queue_length: int | None = None,
force_device_scalars: bool | None = None) -> None:
"""
See :class:`arraycontext.impl.pyopencl.PyOpenCLArrayContext` for most
arguments.
"""
super().__init__(queue, allocator=allocator,
wait_event_queue_length=wait_event_queue_length,
force_device_scalars=force_device_scalars)
self.mpi_communicator = mpi_communicator
def clone(self) -> Self:
return type(self)(self.mpi_communicator, self.queue,
allocator=self.allocator,
wait_event_queue_length=self._wait_event_queue_length,
force_device_scalars=self._force_device_scalars)
# }}}
# {{{ distributed + numpy
class MPINumpyArrayContext(NumpyArrayContext, MPIBasedArrayContext):
"""An array context for using distributed computation with :mod:`numpy`
eager evaluation.
.. autofunction:: __init__
"""
def __init__(self, mpi_communicator) -> None:
super().__init__()
self.mpi_communicator = mpi_communicator
def clone(self) -> Self:
return type(self)(self.mpi_communicator)
# }}}
# {{{ distributed + pytato array context subclasses
class MPIBasePytatoPyOpenCLArrayContext(
MPIBasedArrayContext, PytatoPyOpenCLArrayContext):
"""
.. autofunction:: __init__
"""
def __init__(
self, mpi_communicator, queue, *, mpi_base_tag, allocator=None,
compile_trace_callback: Callable[[Any, str, Any], None] | None = None,
) -> None:
"""
:arg compile_trace_callback: A function of three arguments
*(what, stage, ir)*, where *what* identifies the object
being compiled, *stage* is a string describing the compilation
pass, and *ir* is an object containing the intermediate
representation. This interface should be considered
unstable.
"""
if allocator is None:
warn("No memory allocator specified, please pass one. "
"(Preferably a pyopencl.tools.MemoryPool in order "
"to reduce device allocations)", stacklevel=2)
super().__init__(queue, allocator,
compile_trace_callback=compile_trace_callback)
self.mpi_communicator = mpi_communicator
self.mpi_base_tag = mpi_base_tag
# FIXME: implement distributed-aware freeze
def compile(self, f: Callable[..., Any]) -> Callable[..., Any]:
return _DistributedLazilyPyOpenCLCompilingFunctionCaller(self, f)
def clone(self) -> Self:
return type(self)(self.mpi_communicator, self.queue,
mpi_base_tag=self.mpi_base_tag,
allocator=self.allocator)
MPIPytatoArrayContext: type[MPIBasedArrayContext] = MPIBasePytatoPyOpenCLArrayContext
if _HAVE_FUSION_ACTX:
class MPIFusionContractorArrayContext(
MPIBasePytatoPyOpenCLArrayContext, FusionContractorArrayContext):
"""
.. autofunction:: __init__
"""
MPIPytatoArrayContext = MPIFusionContractorArrayContext
else:
MPIPytatoArrayContext = MPIBasePytatoPyOpenCLArrayContext
# }}}
# {{{ pytest actx factory
class PytestPyOpenCLArrayContextFactory(
_PytestPyOpenCLArrayContextFactoryWithClass):
actx_class = PyOpenCLArrayContext
def __call__(self):
from pyopencl.tools import ImmediateAllocator, MemoryPool
_ctx, queue = self.get_command_queue()
alloc = MemoryPool(ImmediateAllocator(queue))
return self.actx_class(queue, allocator=alloc)
class PytestPytatoPyOpenCLArrayContextFactory(
_PytestPytatoPyOpenCLArrayContextFactory):
actx_class = PytatoPyOpenCLArrayContext
def __call__(self):
_ctx, queue = self.get_command_queue()
from pyopencl.tools import ImmediateAllocator, MemoryPool
alloc = MemoryPool(ImmediateAllocator(queue))
return self.actx_class(queue, allocator=alloc)
class PytestNumpyArrayContextFactory(_PytestNumpyArrayContextFactory):
actx_class = NumpyArrayContext
def __call__(self):
return self.actx_class()
register_pytest_array_context_factory("grudge.pyopencl",
PytestPyOpenCLArrayContextFactory)
register_pytest_array_context_factory("grudge.pytato-pyopencl",
PytestPytatoPyOpenCLArrayContextFactory)
register_pytest_array_context_factory("grudge.numpy",
PytestNumpyArrayContextFactory)
# }}}
# {{{ actx selection
def _get_single_grid_pytato_actx_class(distributed: bool) -> type[ArrayContext]:
warn("No device-parallel actx available, execution will be slow.",
stacklevel=1)
# lazy, non-distributed
if not distributed:
return PytatoPyOpenCLArrayContext
else:
# distributed+lazy:
return MPIBasePytatoPyOpenCLArrayContext
def get_reasonable_array_context_class(
lazy: bool = True, distributed: bool = True,
fusion: bool | None = None, numpy: bool = False,
) -> type[ArrayContext]:
"""Returns a reasonable :class:`~arraycontext.ArrayContext` currently
supported given the constraints of *lazy*, *distributed*, and *numpy*."""
if fusion is None:
fusion = lazy
if numpy:
assert not (lazy or fusion)
if distributed:
actx_class: type[ArrayContext] = MPINumpyArrayContext
else:
actx_class = NumpyArrayContext
return actx_class
if lazy:
if fusion:
if not _HAVE_FUSION_ACTX:
warn("No device-parallel actx available, execution will be slow. "
"Please make sure you have the right branches for loopy "
"(https://github.com/kaushikcfd/loopy/tree/pytato-array-context-transforms) " # noqa
"and meshmode "
"(https://github.com/kaushikcfd/meshmode/tree/pytato-array-context-transforms).",
stacklevel=1)
# lazy+fusion, non-distributed
if _HAVE_FUSION_ACTX:
if distributed:
actx_class = MPIFusionContractorArrayContext
else:
actx_class = FusionContractorArrayContext
else:
actx_class = _get_single_grid_pytato_actx_class(distributed)
else:
actx_class = _get_single_grid_pytato_actx_class(distributed)
else:
if fusion:
raise ValueError("No eager actx's support op-fusion.")
if distributed:
actx_class = MPIPyOpenCLArrayContext
else:
actx_class = PyOpenCLArrayContext
logger.info("get_reasonable_array_context_class: %s lazy=%r distributed=%r "
"device-parallel=%r",
actx_class.__name__, lazy, distributed,
# eager is always device-parallel:
(_HAVE_FUSION_ACTX or not lazy))
return actx_class
# }}}
# vim: foldmethod=marker
grudge/discretization.py
View file @ 49b9b7c7
from __future__ import division, absolute_import
"""
.. autoclass:: DiscretizationTag
.. currentmodule:: grudge
.. autoclass:: DiscretizationCollection
.. autofunction:: make_discretization_collection
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
.. currentmodule:: grudge.discretization
"""
__copyright__ = """
Copyright (C) 2015-2017 Andreas Kloeckner, 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
......@@ -22,202 +33,832 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
from collections.abc import Mapping
from typing import TYPE_CHECKING, Any, Optional, cast
from warnings import warn
import numpy as np
from arraycontext import ArrayContext
from meshmode.discretization import Discretization, ElementGroupFactory
from meshmode.discretization.connection import (
FACE_RESTR_ALL,
FACE_RESTR_INTERIOR,
DirectDiscretizationConnection,
DiscretizationConnection,
make_face_restriction,
)
from meshmode.discretization.poly_element import (
InterpolatoryEdgeClusteredGroupFactory,
ModalGroupFactory,
)
from meshmode.dof_array import DOFArray
from meshmode.mesh import BTAG_PARTITION, Mesh, ModepyElementGroup
from pytools import memoize_method, single_valued
from pytools import memoize_method
from grudge import sym
from grudge.dof_desc import (
DD_VOLUME_ALL,
DISCR_TAG_BASE,
DISCR_TAG_MODAL,
VTAG_ALL,
BoundaryDomainTag,
ConvertibleToDOFDesc,
DiscretizationTag,
DOFDesc,
DomainTag,
VolumeDomainTag,
VolumeTag,
as_dofdesc,
)
class Discretization(object):
"""
.. attribute :: volume_discr
if TYPE_CHECKING:
import mpi4py.MPI
MeshOrDiscr = Mesh | Discretization
TagToElementGroupFactory = Mapping[DiscretizationTag, ElementGroupFactory]
# {{{ discr_tag_to_group_factory normalization
def _normalize_discr_tag_to_group_factory(
discr_tag_to_group_factory: TagToElementGroupFactory | None,
order: int | None
) -> TagToElementGroupFactory:
if discr_tag_to_group_factory is None:
if order is None:
raise TypeError(
"one of 'order' and 'discr_tag_to_group_factory' must be given"
)
discr_tag_to_group_factory = {
DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order=order)}
else:
discr_tag_to_group_factory = dict(discr_tag_to_group_factory)
if order is not None:
if DISCR_TAG_BASE in discr_tag_to_group_factory:
raise ValueError(
"if 'order' is given, 'discr_tag_to_group_factory' must "
"not have a key of DISCR_TAG_BASE"
)
discr_tag_to_group_factory[DISCR_TAG_BASE] = \
InterpolatoryEdgeClusteredGroupFactory(order)
assert discr_tag_to_group_factory is not None
# Modal discr should always come from the base discretization
if DISCR_TAG_MODAL not in discr_tag_to_group_factory and order is not None:
discr_tag_to_group_factory[DISCR_TAG_MODAL] = ModalGroupFactory(order)
return discr_tag_to_group_factory
# }}}
class DiscretizationCollection:
"""A collection of discretizations, defined on the same underlying
:class:`~meshmode.mesh.Mesh`, corresponding to various mesh entities
(volume, interior facets, boundaries) and associated element
groups.
.. automethod :: boundary_connection
.. automethod :: boundary_discr
.. note::
.. automethod :: interior_faces_connection
.. automethod :: interior_faces_discr
Do not call the constructor directly. Use
:func:`make_discretization_collection` instead.
.. automethod :: all_faces_connection
.. automethod :: all_faces_discr
.. autoattribute:: dim
.. autoattribute:: ambient_dim
.. autoattribute:: real_dtype
.. autoattribute:: complex_dtype
.. automethod:: discr_from_dd
.. automethod:: connection_from_dds
.. autoattribute :: cl_context
.. autoattribute :: dim
.. autoattribute :: ambient_dim
.. autoattribute :: mesh
.. automethod:: empty
.. automethod:: zeros
.. automethod:: nodes
.. automethod:: normal
.. rubric:: Internal functionality
.. automethod :: empty
.. automethod :: zeros
.. automethod:: _base_to_geoderiv_connection
"""
def __init__(self, cl_ctx, mesh, order, quad_min_degrees=None):
# {{{ constructor
def __init__(self, array_context: ArrayContext,
volume_discrs: Mesh | Mapping[VolumeTag, Discretization],
order: int | None = None,
discr_tag_to_group_factory: TagToElementGroupFactory | None = None,
mpi_communicator: Optional["mpi4py.MPI.Intracomm"] = None,
) -> None:
"""
:param quad_min_degrees: A mapping from quadrature tags to the degrees
to which the desired quadrature is supposed to be exact.
:arg discr_tag_to_group_factory: A mapping from discretization tags
(typically one of: :class:`grudge.dof_desc.DISCR_TAG_BASE`,
:class:`grudge.dof_desc.DISCR_TAG_MODAL`, or
:class:`grudge.dof_desc.DISCR_TAG_QUAD`) to a
:class:`~meshmode.discretization.ElementGroupFactory`
indicating with which type of discretization the operations are
to be carried out, or *None* to indicate that operations with this
discretization tag should be carried out with the standard volume
discretization.
"""
if quad_min_degrees is None:
quad_min_degrees = {}
self._setup_actx = array_context.clone()
# {{{ process mpi_communicator argument
if mpi_communicator is not None:
warn("Passing 'mpi_communicator' is deprecated. This will stop working "
"in 2023. Instead, pass an MPIBasedArrayContext.",
DeprecationWarning, stacklevel=2)
from grudge.array_context import MPIBasedArrayContext
if (isinstance(array_context, MPIBasedArrayContext)
and mpi_communicator is not array_context.mpi_communicator):
raise ValueError("mpi_communicator passed to "
"DiscretizationCollection and the MPI communicator "
"used to created the MPIBasedArrayContext must be "
"idetical, which they aren't.")
else:
from grudge.array_context import MPIBasedArrayContext
if isinstance(self._setup_actx, MPIBasedArrayContext):
mpi_communicator = self._setup_actx.mpi_communicator
self._mpi_communicator = mpi_communicator
# }}}
from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
self.volume_discr = Discretization(cl_ctx, mesh,
PolynomialWarpAndBlendGroupFactory(order=order))
if isinstance(volume_discrs, Mesh):
# {{{ deprecated backward compatibility yuck
self.order = order
self.quad_min_degrees = quad_min_degrees
warn("Calling the DiscretizationCollection constructor directly "
"is deprecated, call make_discretization_collection "
"instead. This will stop working in 2023.",
DeprecationWarning, stacklevel=2)
# {{{ management of discretization-scoped common subexpressions
mesh = volume_discrs
from pytools import UniqueNameGenerator
self._discr_scoped_name_gen = UniqueNameGenerator()
discr_tag_to_group_factory = _normalize_discr_tag_to_group_factory(
discr_tag_to_group_factory=discr_tag_to_group_factory,
order=order)
self._discr_tag_to_group_factory = discr_tag_to_group_factory
self._discr_scoped_subexpr_to_name = {}
self._discr_scoped_subexpr_name_to_value = {}
volume_discr = Discretization(
array_context, mesh,
self.group_factory_for_discretization_tag(DISCR_TAG_BASE))
volume_discrs = {VTAG_ALL: volume_discr}
# }}}
del mesh
# {{{ boundary
# }}}
else:
assert discr_tag_to_group_factory is not None
self._discr_tag_to_group_factory = discr_tag_to_group_factory
@memoize_method
def boundary_connection(self, boundary_tag, quadrature_tag=None):
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
self._volume_discrs = volume_discrs
self._dist_boundary_connections = {
vtag: self._set_up_distributed_communication(
vtag, mpi_communicator, array_context)
for vtag in self._volume_discrs.keys()}
if quadrature_tag is not sym.QTAG_NONE:
# }}}
@property
def mpi_communicator(self):
warn("Accessing DiscretizationCollection.mpi_communicator is deprecated. "
"This will stop working in 2023. "
"Instead, use an MPIBasedArrayContext, and obtain the communicator "
"from that.", DeprecationWarning, stacklevel=2)
return self._mpi_communicator
def get_management_rank_index(self):
return 0
def is_management_rank(self):
if self.mpi_communicator is None:
return True
else:
return self.mpi_communicator.Get_rank() \
== self.get_management_rank_index()
# {{{ distributed
def _set_up_distributed_communication(
self, vtag, mpi_communicator, array_context):
from_dd = DOFDesc(VolumeDomainTag(vtag), DISCR_TAG_BASE)
boundary_connections = {}
from meshmode.distributed import get_connected_parts
connected_parts = get_connected_parts(self._volume_discrs[vtag].mesh)
if connected_parts:
if mpi_communicator is None:
raise RuntimeError("must supply an MPI communicator when using a "
"distributed mesh")
grp_factory = \
self.group_factory_for_discretization_tag(DISCR_TAG_BASE)
local_boundary_connections = {}
for i_remote_part in connected_parts:
local_boundary_connections[i_remote_part] = self.connection_from_dds(
from_dd, from_dd.trace(BTAG_PARTITION(i_remote_part)))
from meshmode.distributed import MPIBoundaryCommSetupHelper
with MPIBoundaryCommSetupHelper(mpi_communicator, array_context,
local_boundary_connections, grp_factory) as bdry_setup_helper:
while True:
conns = bdry_setup_helper.complete_some()
if not conns:
break
for i_remote_part, conn in conns.items():
boundary_connections[i_remote_part] = conn
return boundary_connections
def distributed_boundary_swap_connection(self, dd):
"""Provides a mapping from the base volume discretization
to the exterior boundary restriction on a parallel boundary
partition described by *dd*. This connection is used to
communicate across element boundaries in different parallel
partitions during distributed runs.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
convertible to one. The domain tag must be a subclass
of :class:`grudge.dof_desc.BoundaryDomainTag` with an
associated :class:`meshmode.mesh.BTAG_PARTITION`
corresponding to a particular communication rank.
"""
if dd.discretization_tag is not DISCR_TAG_BASE:
# FIXME
raise NotImplementedError("quadrature")
raise NotImplementedError(
"Distributed communication with discretization tag "
f"{dd.discretization_tag} is not implemented."
)
from meshmode.discretization.connection import make_face_restriction
return make_face_restriction(
self.volume_discr,
PolynomialWarpAndBlendGroupFactory(order=self.order),
boundary_tag=boundary_tag)
assert isinstance(dd.domain_tag, BoundaryDomainTag)
assert isinstance(dd.domain_tag.tag, BTAG_PARTITION)
vtag = dd.domain_tag.volume_tag
def boundary_discr(self, boundary_tag, quadrature_tag=None):
return self.boundary_connection(boundary_tag, quadrature_tag).to_discr
return self._dist_boundary_connections[vtag][dd.domain_tag.tag.part_id]
# }}}
# {{{ interior faces
# {{{ discr_from_dd
@memoize_method
def interior_faces_connection(self, quadrature_tag=None):
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
def discr_from_dd(self, dd: "ConvertibleToDOFDesc") -> Discretization:
"""Provides a :class:`meshmode.discretization.Discretization`
object from *dd*.
"""
dd = as_dofdesc(dd)
if quadrature_tag is not sym.QTAG_NONE:
# FIXME
raise NotImplementedError("quadrature")
discr_tag = dd.discretization_tag
from meshmode.discretization.connection import (
make_face_restriction, FRESTR_INTERIOR_FACES)
return make_face_restriction(
self.volume_discr,
PolynomialWarpAndBlendGroupFactory(order=self.order),
FRESTR_INTERIOR_FACES,
if discr_tag is DISCR_TAG_MODAL:
return self._modal_discr(dd.domain_tag)
if dd.is_volume():
return self._volume_discr_from_dd(dd)
# FIXME: This will need to change as soon as we support
# pyramids or other elements with non-identical face
# types.
per_face_groups=False)
if discr_tag is not DISCR_TAG_BASE:
base_discr = self.discr_from_dd(dd.with_discr_tag(DISCR_TAG_BASE))
def interior_faces_discr(self, quadrature_tag=None):
return self.interior_faces_connection(quadrature_tag).to_discr
from meshmode.discretization import Discretization
return Discretization(
self._setup_actx,
base_discr.mesh,
self.group_factory_for_discretization_tag(discr_tag)
)
assert discr_tag is DISCR_TAG_BASE
if isinstance(dd.domain_tag, BoundaryDomainTag):
if dd.domain_tag.tag in [FACE_RESTR_ALL, FACE_RESTR_INTERIOR]:
return self._faces_connection(dd.domain_tag).to_discr
else:
return self._boundary_connection(dd.domain_tag).to_discr
else:
raise ValueError(f"DOF desc not understood: {dd}")
# }}}
# {{{ _base_to_geoderiv_connection
@memoize_method
def opposite_face_connection(self, quadrature_tag):
if quadrature_tag is not sym.QTAG_NONE:
# FIXME
raise NotImplementedError("quadrature")
def _has_affine_groups(self, domain_tag: DomainTag) -> bool:
from modepy.shapes import Simplex
discr = self.discr_from_dd(DOFDesc(domain_tag, DISCR_TAG_BASE))
return any(
megrp.is_affine
and issubclass(cast(ModepyElementGroup, megrp).shape_cls, Simplex)
for megrp in discr.mesh.groups)
from meshmode.discretization.connection import \
make_opposite_face_connection
@memoize_method
def _base_to_geoderiv_connection(self, dd: DOFDesc) -> DiscretizationConnection:
r"""The "geometry derivatives" discretization for a given *dd* is
typically identical to the one returned by :meth:`discr_from_dd`,
however for affinely-mapped simplicial elements, it will use a
:math:`P^0` discretization having a single DOF per element.
As a result, :class:`~meshmode.dof_array.DOFArray`\ s on this
are broadcast-compatible with the discretizations returned by
:meth:`discr_from_dd`.
This is an internal function, not intended for use outside
:mod:`grudge`.
"""
base_discr = self.discr_from_dd(dd)
if not self._has_affine_groups(dd.domain_tag):
# no benefit to having another discretization that takes
# advantage of affine-ness
from meshmode.discretization.connection import (
IdentityDiscretizationConnection,
)
return IdentityDiscretizationConnection(base_discr)
base_group_factory = self.group_factory_for_discretization_tag(
dd.discretization_tag)
def geo_group_factory(megrp):
from meshmode.discretization.poly_element import (
PolynomialEquidistantSimplexElementGroup,
)
from modepy.shapes import Simplex
if megrp.is_affine and issubclass(megrp._modepy_shape_cls, Simplex):
return PolynomialEquidistantSimplexElementGroup(
megrp, order=0)
else:
return base_group_factory(megrp)
return make_opposite_face_connection(
self.interior_faces_connection(quadrature_tag))
from meshmode.discretization import Discretization
geo_deriv_discr = Discretization(
self._setup_actx, base_discr.mesh,
geo_group_factory)
from meshmode.discretization.connection.same_mesh import (
make_same_mesh_connection,
)
return make_same_mesh_connection(
self._setup_actx,
to_discr=geo_deriv_discr,
from_discr=base_discr)
# }}}
# {{{ all-faces
# {{{ connection_from_dds
@memoize_method
def all_faces_volume_connection(self, quadrature_tag=None):
if quadrature_tag is not sym.QTAG_NONE:
# FIXME
raise NotImplementedError("quadrature")
def connection_from_dds(
self, from_dd: "ConvertibleToDOFDesc", to_dd: "ConvertibleToDOFDesc"
) -> DiscretizationConnection:
"""Provides a mapping (connection) from one discretization to
another, e.g. from the volume to the boundary, or from the
base to the an overintegrated quadrature discretization, or from
a nodal representation to a modal representation.
:arg from_dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
convertible to one.
:arg to_dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value
convertible to one.
"""
from_dd = as_dofdesc(from_dd)
to_dd = as_dofdesc(to_dd)
to_discr_tag = to_dd.discretization_tag
from_discr_tag = from_dd.discretization_tag
# {{{ mapping between modal and nodal representations
if (from_discr_tag is DISCR_TAG_MODAL
and to_discr_tag is not DISCR_TAG_MODAL):
return self._modal_to_nodal_connection(to_dd)
if (to_discr_tag is DISCR_TAG_MODAL
and from_discr_tag is not DISCR_TAG_MODAL):
return self._nodal_to_modal_connection(from_dd)
# }}}
assert (to_discr_tag is not DISCR_TAG_MODAL
and from_discr_tag is not DISCR_TAG_MODAL)
if (isinstance(from_dd.domain_tag, BoundaryDomainTag)
and from_discr_tag == to_discr_tag
and isinstance(to_dd.domain_tag, BoundaryDomainTag)
and to_dd.domain_tag.tag is FACE_RESTR_ALL):
faces_conn = self.connection_from_dds(
DOFDesc(
VolumeDomainTag(from_dd.domain_tag.volume_tag),
DISCR_TAG_BASE),
from_dd.with_discr_tag(DISCR_TAG_BASE))
from meshmode.discretization.connection import (
make_face_to_all_faces_embedding,
)
assert isinstance(faces_conn, DirectDiscretizationConnection)
return make_face_to_all_faces_embedding(
self._setup_actx,
faces_conn, self.discr_from_dd(to_dd),
self.discr_from_dd(from_dd))
# {{{ simplify domain + discr_tag change into chained
if (from_dd.domain_tag != to_dd.domain_tag
and from_discr_tag is DISCR_TAG_BASE
and to_discr_tag is not DISCR_TAG_BASE):
from meshmode.discretization.connection import (
ChainedDiscretizationConnection,
)
intermediate_dd = to_dd.with_discr_tag(DISCR_TAG_BASE)
return ChainedDiscretizationConnection(
[
# first change domain
self.connection_from_dds(
from_dd,
intermediate_dd),
# then go to quad grid
self.connection_from_dds(
intermediate_dd,
to_dd
)])
# }}}
# {{{ generic to-quad
# DISCR_TAG_MODAL is handled above
if (from_dd.domain_tag == to_dd.domain_tag
and from_discr_tag is DISCR_TAG_BASE
and to_discr_tag is not DISCR_TAG_BASE):
from meshmode.discretization.connection.same_mesh import (
make_same_mesh_connection,
)
return make_same_mesh_connection(
self._setup_actx,
self.discr_from_dd(to_dd),
self.discr_from_dd(from_dd))
# }}}
if from_discr_tag is not DISCR_TAG_BASE:
raise ValueError("cannot get a connection *from* a "
f"(non-interpolatory) quadrature grid: '{from_dd}'")
assert to_discr_tag is DISCR_TAG_BASE
if isinstance(from_dd.domain_tag, VolumeDomainTag):
if isinstance(to_dd.domain_tag, BoundaryDomainTag):
if to_dd.domain_tag.volume_tag != from_dd.domain_tag.tag:
raise ValueError("cannot get a connection from one volume "
f"('{from_dd.domain_tag.tag}') "
"to the boundary of another volume "
f"('{to_dd.domain_tag.volume_tag}') ")
if to_dd.domain_tag.tag in [FACE_RESTR_ALL, FACE_RESTR_INTERIOR]:
return self._faces_connection(to_dd.domain_tag)
elif isinstance(to_dd.domain_tag, BoundaryDomainTag):
assert from_discr_tag is DISCR_TAG_BASE
return self._boundary_connection(to_dd.domain_tag)
elif to_dd.is_volume():
if to_dd.domain_tag != from_dd.domain_tag:
raise ValueError("cannot get a connection between "
"volumes of different tags: requested "
f"'{from_dd.domain_tag}' -> '{to_dd.domain_tag}'")
from meshmode.discretization.connection import make_same_mesh_connection
return make_same_mesh_connection(
self._setup_actx,
self._volume_discr_from_dd(to_dd),
self._volume_discr_from_dd(from_dd))
else:
raise ValueError(
f"cannot get a connection from volume to: '{to_dd}'")
else:
raise ValueError(f"cannot get a connection from: '{from_dd}'")
# }}}
# {{{ group_factory_for_discretization_tag
def group_factory_for_discretization_tag(self, discretization_tag):
if discretization_tag is None:
discretization_tag = DISCR_TAG_BASE
return self._discr_tag_to_group_factory[discretization_tag]
# }}}
# {{{ (internal) discretization getters
@memoize_method
def _volume_discr_from_dd(self, dd: DOFDesc) -> Discretization:
assert isinstance(dd.domain_tag, VolumeDomainTag)
try:
base_volume_discr = self._volume_discrs[dd.domain_tag.tag]
except KeyError:
raise ValueError("a volume discretization with volume tag "
f"'{dd.domain_tag.tag}' is not known") from None
# Refuse to re-make the volume discretization
if dd.discretization_tag is DISCR_TAG_BASE:
return base_volume_discr
from meshmode.discretization import Discretization
return Discretization(
self._setup_actx, base_volume_discr.mesh,
self.group_factory_for_discretization_tag(dd.discretization_tag)
)
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
@memoize_method
def _modal_discr(self, domain_tag) -> Discretization:
from meshmode.discretization import Discretization
discr_base = self.discr_from_dd(DOFDesc(domain_tag, DISCR_TAG_BASE))
return Discretization(
self._setup_actx, discr_base.mesh,
self.group_factory_for_discretization_tag(DISCR_TAG_MODAL)
)
# }}}
# {{{ connection factories: modal<->nodal
@memoize_method
def _modal_to_nodal_connection(self, to_dd: DOFDesc) -> DiscretizationConnection:
"""
:arg to_dd: a :class:`grudge.dof_desc.DOFDesc`
describing the dofs corresponding to the
*to_discr*
"""
from meshmode.discretization.connection import (
make_face_restriction, FRESTR_ALL_FACES)
ModalToNodalDiscretizationConnection,
)
return ModalToNodalDiscretizationConnection(
from_discr=self._modal_discr(to_dd.domain_tag),
to_discr=self.discr_from_dd(to_dd)
)
@memoize_method
def _nodal_to_modal_connection(
self, from_dd: DOFDesc) -> DiscretizationConnection:
"""
:arg from_dd: a :class:`grudge.dof_desc.DOFDesc`
describing the dofs corresponding to the
*from_discr*
"""
from meshmode.discretization.connection import (
NodalToModalDiscretizationConnection,
)
return NodalToModalDiscretizationConnection(
from_discr=self.discr_from_dd(from_dd),
to_discr=self._modal_discr(from_dd.domain_tag)
)
# }}}
# {{{ connection factories: boundary
@memoize_method
def _boundary_connection(
self, domain_tag: BoundaryDomainTag) -> DiscretizationConnection:
return make_face_restriction(
self.volume_discr,
PolynomialWarpAndBlendGroupFactory(order=self.order),
FRESTR_ALL_FACES,
self._setup_actx,
self._volume_discr_from_dd(
DOFDesc(VolumeDomainTag(domain_tag.volume_tag), DISCR_TAG_BASE)),
self.group_factory_for_discretization_tag(DISCR_TAG_BASE),
boundary_tag=domain_tag.tag)
# FIXME: This will need to change as soon as we support
# pyramids or other elements with non-identical face
# types.
per_face_groups=False)
# }}}
def all_faces_discr(self, quadrature_tag=None):
return self.all_faces_volume_connection(quadrature_tag).to_discr
# {{{ connection factories: faces
@memoize_method
def all_faces_connection(self, boundary_tag, quadrature_tag=None):
"""Return a
:class:`meshmode.discretization.connection.DiscretizationConnection`
that goes from either
:meth:`interior_faces_discr` (if *boundary_tag* is None)
or
:meth:`boundary_discr` (if *boundary_tag* is not None)
to a discretization containing all the faces of the volume
discretization.
def _faces_connection(
self, domain_tag: BoundaryDomainTag) -> DiscretizationConnection:
assert domain_tag.tag in [FACE_RESTR_INTERIOR, FACE_RESTR_ALL]
return make_face_restriction(
self._setup_actx,
self._volume_discr_from_dd(
DOFDesc(
VolumeDomainTag(domain_tag.volume_tag), DISCR_TAG_BASE)),
self.group_factory_for_discretization_tag(DISCR_TAG_BASE),
domain_tag.tag,
# FIXME: This will need to change as soon as we support
# pyramids or other elements with non-identical face
# types.
per_face_groups=False
)
@memoize_method
def opposite_face_connection(
self, domain_tag: BoundaryDomainTag) -> DiscretizationConnection:
"""Provides a mapping from the base volume discretization
to the exterior boundary restriction on a neighboring element.
This does not take into account parallel partitions.
"""
from meshmode.discretization.connection import \
make_face_to_all_faces_embedding
from meshmode.discretization.connection import make_opposite_face_connection
if boundary_tag is None:
faces_conn = self.interior_faces_connection(quadrature_tag)
else:
faces_conn = self.boundary_connection(boundary_tag, quadrature_tag)
assert domain_tag.tag is FACE_RESTR_INTERIOR
return make_face_to_all_faces_embedding(
faces_conn, self.all_faces_discr(quadrature_tag))
return make_opposite_face_connection(
self._setup_actx,
self._faces_connection(domain_tag))
# }}}
# {{{ properties
@property
def cl_context(self):
return self.volume_discr.cl_context
def dim(self) -> int:
"""Return the topological dimension."""
return single_valued(discr.dim for discr in self._volume_discrs.values())
@property
def dim(self):
return self.volume_discr.dim
def ambient_dim(self) -> int:
"""Return the dimension of the ambient space."""
return single_valued(
discr.ambient_dim for discr in self._volume_discrs.values())
@property
def ambient_dim(self):
return self.volume_discr.ambient_dim
def real_dtype(self) -> "np.dtype[Any]":
"""Return the data type used for real-valued arithmetic."""
return single_valued(
discr.real_dtype for discr in self._volume_discrs.values())
@property
def mesh(self):
return self.volume_discr.mesh
def complex_dtype(self) -> "np.dtype[Any]":
"""Return the data type used for complex-valued arithmetic."""
return single_valued(
discr.complex_dtype for discr in self._volume_discrs.values())
def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None):
return self.volume_discr.empty(queue, dtype, extra_dims=extra_dims,
allocator=allocator)
# }}}
def zeros(self, queue, dtype=None, extra_dims=None, allocator=None):
return self.volume_discr.zeros(queue, dtype, extra_dims=extra_dims,
allocator=allocator)
# {{{ array creators
def empty(self, array_context: ArrayContext, dtype=None,
*, dd: DOFDesc | None = None) -> DOFArray:
"""Return an empty :class:`~meshmode.dof_array.DOFArray` defined at
the volume nodes: :class:`grudge.dof_desc.DD_VOLUME_ALL`.
:arg array_context: an :class:`~arraycontext.context.ArrayContext`.
:arg dtype: type special value 'c' will result in a
vector of dtype :attr:`complex_dtype`. If
*None* (the default), a real vector will be returned.
"""
if dd is None:
dd = DD_VOLUME_ALL
return self.discr_from_dd(dd).empty(array_context, dtype)
def zeros(self, array_context: ArrayContext, dtype=None,
*, dd: DOFDesc | None = None) -> DOFArray:
"""Return a zero-initialized :class:`~meshmode.dof_array.DOFArray`
defined at the volume nodes, :class:`grudge.dof_desc.DD_VOLUME_ALL`.
:arg array_context: an :class:`~arraycontext.context.ArrayContext`.
:arg dtype: type special value 'c' will result in a
vector of dtype :attr:`complex_dtype`. If
*None* (the default), a real vector will be returned.
"""
if dd is None:
dd = DD_VOLUME_ALL
return self.discr_from_dd(dd).zeros(array_context, dtype)
def is_volume_where(self, where):
from grudge import sym
return (
where is None
or where == sym.VTAG_ALL)
return where is None or as_dofdesc(where).is_volume()
# }}}
# {{{ discretization-specific geometric fields
def nodes(self, dd=None):
r"""Return the nodes of a discretization specified by *dd*.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:returns: an object array of frozen :class:`~meshmode.dof_array.DOFArray`\ s
"""
if dd is None:
dd = DD_VOLUME_ALL
return self.discr_from_dd(dd).nodes()
def normal(self, dd):
r"""Get the unit normal to the specified surface discretization, *dd*.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
:returns: an object array of frozen :class:`~meshmode.dof_array.DOFArray`\ s.
"""
warn("DiscretizationCollection.normal is deprecated and "
"will stop working in 2024. "
"Use grudge.geometry.normal instead. "
"DiscretizationCollection.normal may provide non-P^0 normals "
"even when this would otherwise be possible, see "
"https://github.com/inducer/grudge/issues/314 for details.",
DeprecationWarning, stacklevel=2)
from grudge.geometry import normal
return self._setup_actx.freeze(normal(self._setup_actx, self, dd))
# }}}
# {{{ make_discretization_collection
def make_discretization_collection(
array_context: ArrayContext,
volumes: Mesh | Mapping[VolumeTag, Mesh],
order: int | None = None,
discr_tag_to_group_factory: TagToElementGroupFactory | None = None,
) -> DiscretizationCollection:
"""
:arg discr_tag_to_group_factory: A mapping from discretization tags
(typically one of: :class:`~grudge.dof_desc.DISCR_TAG_BASE`,
:class:`~grudge.dof_desc.DISCR_TAG_MODAL`, or
:class:`~grudge.dof_desc.DISCR_TAG_QUAD`) to a
:class:`~meshmode.discretization.ElementGroupFactory`
indicating with which type of discretization the operations are
to be carried out, or *None* to indicate that operations with this
discretization tag should be carried out with the standard volume
discretization.
.. note::
If passing a :class:`~meshmode.discretization.Discretization` in
*volumes*, it must be nodal and unisolvent, consistent with
:class:`~grudge.dof_desc.DISCR_TAG_BASE`.
.. note::
To use the resulting :class:`DiscretizationCollection` in a
distributed-memory manner, the *array_context* passed in
must be one of the distributed-memory array contexts
from :mod:`grudge.array_context`. Unlike the (now-deprecated,
for direct use) constructor of :class:`DiscretizationCollection`,
this function no longer accepts a separate MPI communicator.
.. note::
If the resulting :class:`DiscretizationCollection` is distributed
across multiple ranks, then this is an MPI-collective operation,
i.e. all ranks in the communicator must enter this function at the same
time.
"""
if not isinstance(volumes, Mesh):
volumes_dict = volumes
else:
volumes_dict = {VTAG_ALL: volumes}
from pytools import is_single_valued
assert len(volumes_dict) > 0
if not is_single_valued(mesh_or_discr.ambient_dim
for mesh_or_discr in volumes_dict.values()):
raise ValueError("all parts of a discretization collection must share "
"an ambient dimension")
discr_tag_to_group_factory = _normalize_discr_tag_to_group_factory(
discr_tag_to_group_factory=discr_tag_to_group_factory,
order=order)
del order
volume_discrs = {
vtag: Discretization(
array_context,
mesh,
discr_tag_to_group_factory[DISCR_TAG_BASE])
for vtag, mesh in volumes_dict.items()}
return DiscretizationCollection(
array_context=array_context,
volume_discrs=volume_discrs,
discr_tag_to_group_factory=discr_tag_to_group_factory)
# }}}
# vim: foldmethod=marker
grudge/dof_desc.py 0 → 100644
View file @ 49b9b7c7
"""
Volume tags
-----------
.. autoclass:: VolumeTag
.. autoclass:: VTAG_ALL
:mod:`grudge`-specific boundary tags
------------------------------------
Domain tags
-----------
A domain tag identifies a geometric part (or whole) of the domain described
by a :class:`grudge.DiscretizationCollection`. This can be a volume or a boundary.
.. autoclass:: DTAG_SCALAR
.. autoclass:: DTAG_VOLUME_ALL
.. autoclass:: VolumeDomainTag
.. autoclass:: BoundaryDomainTag
Discretization tags
-------------------
A discretization tag serves as a symbolic identifier of the manner in which
meaning is assigned to degrees of freedom.
.. autoclass:: DISCR_TAG_BASE
.. autoclass:: DISCR_TAG_QUAD
.. autoclass:: DISCR_TAG_MODAL
DOF Descriptor
--------------
.. autoclass:: DOFDesc
.. autofunction:: as_dofdesc
Shortcuts
---------
.. data:: DD_SCALAR
.. data:: DD_VOLUME_ALL
.. data:: DD_VOLUME_ALL_MODAL
Internal things that are visible due to type annotations
--------------------------------------------------------
.. class:: _DiscretizationTag
.. class:: ConvertibleToDOFDesc
Anything that is convertible to a :class:`DOFDesc` via :func:`as_dofdesc`.
"""
__copyright__ = """
Copyright (C) 2008 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.
"""
from collections.abc import Hashable
from dataclasses import dataclass, replace
from typing import Any
from warnings import warn
from meshmode.discretization.connection import FACE_RESTR_ALL, FACE_RESTR_INTERIOR
from meshmode.mesh import (
BTAG_ALL,
BTAG_NONE,
BTAG_PARTITION,
BTAG_REALLY_ALL,
BoundaryTag,
)
from pytools import to_identifier
# {{{ volume tags
class VTAG_ALL: # noqa: N801
pass
VolumeTag = Hashable
# }}}
# {{{ domain tag
@dataclass(frozen=True, eq=True)
class ScalarDomainTag:
"""A domain tag denoting scalar values."""
DTAG_SCALAR = ScalarDomainTag()
@dataclass(frozen=True, eq=True, init=True)
class VolumeDomainTag:
"""A domain tag referring to a volume identified by the volume tag :attr:`tag`.
.. attribute:: tag
.. automethod:: __init__
"""
tag: VolumeTag
DTAG_VOLUME_ALL = VolumeDomainTag(VTAG_ALL)
@dataclass(frozen=True, eq=True, init=True)
class BoundaryDomainTag:
"""A domain tag referring to a boundary identified by the
boundary tag :attr:`tag`.
.. attribute:: tag
.. attribute:: volume_tag
.. automethod:: __init__
"""
tag: BoundaryTag
volume_tag: VolumeTag = VTAG_ALL
DomainTag = ScalarDomainTag | VolumeDomainTag | BoundaryDomainTag
# }}}
# {{{ discretization tag
class _DiscretizationTag:
pass
DiscretizationTag = type[_DiscretizationTag]
class DISCR_TAG_BASE(_DiscretizationTag): # noqa: N801
"""A discretization tag indicating the use of a
nodal and unisolvent discretization. This tag is used
to distinguish the base discretization from quadrature
(e.g. overintegration) or modal (:class:`DISCR_TAG_MODAL`)
discretizations.
"""
class DISCR_TAG_QUAD(_DiscretizationTag): # noqa: N801
"""A discretization tag indicating the use of a quadrature discretization
grid, which typically affords higher quadrature accuracy (e.g. for
nonlinear terms) at the expense of unisolvency. This tag is used to
distinguish the quadrature discretization (e.g. overintegration) from modal
(:class:`DISCR_TAG_MODAL`) or base (:class:`DISCR_TAG_BASE`)
discretizations.
For working with multiple quadrature grids, it is
recommended to create appropriate subclasses of
:class:`DISCR_TAG_QUAD` and define appropriate
:class:`DOFDesc` objects corresponding to each
subclass. For example:
.. code-block:: python
class CustomQuadTag(DISCR_TAG_QUAD):
"A custom quadrature discretization tag."
dd = DOFDesc(DTAG_VOLUME_ALL, CustomQuadTag)
"""
class DISCR_TAG_MODAL(_DiscretizationTag): # noqa: N801
"""A discretization tag indicating the use of unisolvent modal degrees of
freedom. This tag is used to distinguish the modal discretization from the
base (nodal) discretization (e.g. :class:`DISCR_TAG_BASE`) or
discretizations on quadrature grids (:class:`DISCR_TAG_QUAD`).
"""
# }}}
# {{{ DOF descriptor
@dataclass(frozen=True, eq=True)
class DOFDesc:
"""Describes the meaning of degrees of freedom.
.. attribute:: domain_tag
.. attribute:: discretization_tag
.. automethod:: __init__
.. automethod:: is_scalar
.. automethod:: is_discretized
.. automethod:: is_volume
.. automethod:: is_boundary_or_partition_interface
.. automethod:: is_trace
.. automethod:: uses_quadrature
.. automethod:: with_domain_tag
.. automethod:: with_discr_tag
.. automethod:: trace
.. automethod:: untrace
.. automethod:: with_boundary_tag
.. automethod:: __eq__
.. automethod:: __ne__
.. automethod:: __hash__
.. automethod:: as_identifier
"""
domain_tag: DomainTag
discretization_tag: DiscretizationTag
def __init__(self,
domain_tag: Any,
discretization_tag: DiscretizationTag | None = None) -> None:
if (
not isinstance(domain_tag,
ScalarDomainTag | BoundaryDomainTag | VolumeDomainTag)
or discretization_tag is None
or (
not isinstance(discretization_tag, type)
or not issubclass(discretization_tag, _DiscretizationTag))):
warn("Sloppy construction of DOFDesc is deprecated. "
"This will stop working in 2023. "
"Call as_dofdesc instead, with the same arguments. ",
DeprecationWarning, stacklevel=2)
domain_tag, discretization_tag = _normalize_domain_and_discr_tag(
domain_tag, discretization_tag)
object.__setattr__(self, "domain_tag", domain_tag)
object.__setattr__(self, "discretization_tag", discretization_tag)
def is_scalar(self) -> bool:
return isinstance(self.domain_tag, ScalarDomainTag)
def is_discretized(self) -> bool:
return not self.is_scalar()
def is_volume(self) -> bool:
return isinstance(self.domain_tag, VolumeDomainTag)
def is_boundary_or_partition_interface(self) -> bool:
return (isinstance(self.domain_tag, BoundaryDomainTag)
and self.domain_tag.tag not in [
FACE_RESTR_ALL,
FACE_RESTR_INTERIOR])
def is_trace(self) -> bool:
return isinstance(self.domain_tag, BoundaryDomainTag)
def uses_quadrature(self) -> bool:
# FIXME: String tags are deprecated
if isinstance(self.discretization_tag, str):
# All strings are interpreted as quadrature-related tags
return True
elif isinstance(self.discretization_tag, type):
if issubclass(self.discretization_tag, DISCR_TAG_QUAD):
return True
elif issubclass(self.discretization_tag,
DISCR_TAG_BASE | DISCR_TAG_MODAL):
return False
raise ValueError(
f"Invalid discretization tag: {self.discretization_tag}")
def with_domain_tag(self, dtag) -> "DOFDesc":
return replace(self, domain_tag=dtag)
def with_discr_tag(self, discr_tag) -> "DOFDesc":
return replace(self, discretization_tag=discr_tag)
def trace(self, btag: BoundaryTag) -> "DOFDesc":
"""Return a :class:`DOFDesc` for the restriction of the volume
descriptor *self* to the boundary named by *btag*.
An error is raised if this method is called on a non-volume instance of
:class:`DOFDesc`.
"""
if not isinstance(self.domain_tag, VolumeDomainTag):
raise ValueError(f"must originate on volume, got '{self.domain_tag}'")
return replace(self,
domain_tag=BoundaryDomainTag(btag, volume_tag=self.domain_tag.tag))
def untrace(self) -> "DOFDesc":
"""Return a :class:`DOFDesc` for the volume associated with the boundary
descriptor *self*.
An error is raised if this method is called on a non-boundary instance of
:class:`DOFDesc`.
"""
if not isinstance(self.domain_tag, BoundaryDomainTag):
raise ValueError(f"must originate on boundary, got '{self.domain_tag}'")
return replace(self,
domain_tag=VolumeDomainTag(self.domain_tag.volume_tag))
def with_boundary_tag(self, btag: BoundaryTag) -> "DOFDesc":
"""Return a :class:`DOFDesc` representing a boundary named by *btag*
on the same volume as *self*.
An error is raised if this method is called on a non-boundary instance of
:class:`DOFDesc`.
"""
if not isinstance(self.domain_tag, BoundaryDomainTag):
raise ValueError(f"must originate on boundary, got '{self.domain_tag}'")
return replace(self,
domain_tag=replace(self.domain_tag, tag=btag))
def as_identifier(self) -> str:
"""Returns a descriptive string for this :class:`DOFDesc` that is usable
in Python identifiers.
"""
if self.domain_tag is DTAG_SCALAR:
dom_id = "sc"
elif self.domain_tag is DTAG_VOLUME_ALL:
dom_id = "vol"
elif self.domain_tag is FACE_RESTR_ALL:
dom_id = "f_all"
elif self.domain_tag is FACE_RESTR_INTERIOR:
dom_id = "f_int"
elif isinstance(self.domain_tag, VolumeDomainTag):
vtag = self.domain_tag.tag
if isinstance(vtag, type):
vtag = vtag.__name__.replace("VTAG_", "").lower()
elif isinstance(vtag, str):
vtag = to_identifier(vtag)
else:
vtag = to_identifier(str(vtag))
dom_id = f"v_{vtag}"
elif isinstance(self.domain_tag, BoundaryDomainTag):
btag = self.domain_tag.tag
if isinstance(btag, type):
btag = btag.__name__.replace("BTAG_", "").lower()
elif isinstance(btag, str):
btag = to_identifier(btag)
else:
btag = to_identifier(str(btag))
dom_id = f"b_{btag}"
else:
raise ValueError(f"unexpected domain tag: '{self.domain_tag}'")
if isinstance(self.discretization_tag, str):
discr_id = to_identifier(self.discretization_tag)
elif issubclass(self.discretization_tag, DISCR_TAG_QUAD):
discr_id = "_quad"
elif self.discretization_tag is DISCR_TAG_BASE:
discr_id = ""
elif self.discretization_tag is DISCR_TAG_MODAL:
discr_id = "_modal"
else:
raise ValueError(
f"Unexpected discretization tag: {self.discretization_tag}"
)
return f"{dom_id}{discr_id}"
DD_SCALAR = DOFDesc(DTAG_SCALAR, DISCR_TAG_BASE)
DD_VOLUME_ALL = DOFDesc(DTAG_VOLUME_ALL, DISCR_TAG_BASE)
DD_VOLUME_ALL_MODAL = DOFDesc(DTAG_VOLUME_ALL, DISCR_TAG_MODAL)
def _normalize_domain_and_discr_tag(
domain: Any,
discretization_tag: DiscretizationTag | None = None,
*, _contextual_volume_tag: VolumeTag | None = None
) -> tuple[DomainTag, DiscretizationTag]:
contextual_volume_tag = _contextual_volume_tag
del _contextual_volume_tag
if contextual_volume_tag is None:
contextual_volume_tag = VTAG_ALL
if domain == "scalar":
domain = DTAG_SCALAR
elif isinstance(domain, ScalarDomainTag | BoundaryDomainTag | VolumeDomainTag):
pass
elif domain in [VTAG_ALL, "vol"]:
domain = DTAG_VOLUME_ALL
elif domain in [FACE_RESTR_ALL, "all_faces"]:
domain = BoundaryDomainTag(FACE_RESTR_ALL, contextual_volume_tag)
elif domain in [FACE_RESTR_INTERIOR, "int_faces"]:
domain = BoundaryDomainTag(FACE_RESTR_INTERIOR, contextual_volume_tag)
elif isinstance(domain, BTAG_PARTITION):
domain = BoundaryDomainTag(domain, contextual_volume_tag)
elif domain in [BTAG_ALL, BTAG_REALLY_ALL, BTAG_NONE]:
domain = BoundaryDomainTag(domain, contextual_volume_tag)
else:
raise ValueError(f"domain tag not understood: {domain}")
if domain is DTAG_SCALAR and discretization_tag is not None:
raise ValueError("cannot have nontrivial discretization tag on scalar")
if discretization_tag is None:
discretization_tag = DISCR_TAG_BASE
return domain, discretization_tag
ConvertibleToDOFDesc = Any
def as_dofdesc(
domain: "ConvertibleToDOFDesc",
discretization_tag: DiscretizationTag | None = None,
*, _contextual_volume_tag: VolumeTag | None = None) -> DOFDesc:
"""
:arg domain_tag: One of the following:
:class:`DTAG_SCALAR` (or the string ``"scalar"``),
:class:`DTAG_VOLUME_ALL` (or the string ``"vol"``)
for the default volume discretization,
:data:`~meshmode.discretization.connection.FACE_RESTR_ALL`
(or the string ``"all_faces"``), or
:data:`~meshmode.discretization.connection.FACE_RESTR_INTERIOR`
(or the string ``"int_faces"``), or one of
:class:`~meshmode.mesh.BTAG_ALL`,
:class:`~meshmode.mesh.BTAG_NONE`,
:class:`~meshmode.mesh.BTAG_REALLY_ALL`,
:class:`~meshmode.mesh.BTAG_PARTITION`,
or *None* to indicate that the geometry is not yet known.
:arg discretization_tag:
*None* or :class:`DISCR_TAG_BASE` to indicate the use of the basic
discretization grid, :class:`DISCR_TAG_MODAL` to indicate a
modal discretization, or :class:`DISCR_TAG_QUAD` to indicate
the use of a quadrature grid.
"""
if isinstance(domain, DOFDesc):
return domain
domain, discretization_tag = _normalize_domain_and_discr_tag(
domain, discretization_tag,
_contextual_volume_tag=_contextual_volume_tag)
return DOFDesc(domain, discretization_tag)
# }}}
# {{{ deprecations
_deprecated_name_to_new_name = {
"DTAG_VOLUME": "VolumeDomainTag",
"DTAG_BOUNDARY": "BoundaryDomainTag",
"DD_VOLUME": "DD_VOLUME_ALL",
"DD_VOLUME_MODAL": "DD_VOLUME_ALL_MODAL"
}
def __getattr__(name):
if name in _deprecated_name_to_new_name:
warn(f"'{name}' is deprecated and will be dropped "
f"in version 2023.x. Use '{_deprecated_name_to_new_name[name]}' "
"instead.",
DeprecationWarning, stacklevel=2)
return globals()[_deprecated_name_to_new_name[name]]
raise AttributeError(f"module {__name__} has no attribute {name}")
# }}}
# vim: foldmethod=marker
grudge/dt_utils.py 0 → 100644
View file @ 49b9b7c7
"""Helper functions for estimating stable time steps for RKDG methods.
Characteristic lengthscales
---------------------------
.. autofunction:: characteristic_lengthscales
Non-geometric quantities
------------------------
.. autofunction:: dt_non_geometric_factors
Mesh size utilities
-------------------
.. autofunction:: dt_geometric_factors
.. autofunction:: h_max_from_volume
.. autofunction:: h_min_from_volume
"""
__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.
"""
from collections.abc import Sequence
from typing import cast
import numpy as np
from arraycontext import ArrayContext, Scalar, tag_axes
from arraycontext.metadata import NameHint
from meshmode.discretization import NodalElementGroupBase
from meshmode.dof_array import DOFArray
from meshmode.transform_metadata import (
DiscretizationDOFAxisTag,
DiscretizationElementAxisTag,
DiscretizationFaceAxisTag,
FirstAxisIsElementsTag,
)
from pytools import memoize_in, memoize_on_first_arg
import grudge.op as op
from grudge.discretization import DiscretizationCollection
from grudge.dof_desc import (
DD_VOLUME_ALL,
FACE_RESTR_ALL,
BoundaryDomainTag,
DOFDesc,
ScalarDomainTag,
as_dofdesc,
)
def characteristic_lengthscales(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None) -> DOFArray:
r"""Computes the characteristic length scale :math:`h_{\text{loc}}` at
each node. The characteristic length scale is mainly useful for estimating
the stable time step size. E.g. for a hyperbolic system, an estimate of the
stable time step can be estimated as :math:`h_{\text{loc}} / c`, where
:math:`c` is the characteristic wave speed. The estimate is obtained using
the following formula:
.. math::
h_{\text{loc}} = \operatorname{min}\left(\Delta r_i\right) r_D
where :math:`\operatorname{min}\left(\Delta r_i\right)` is the minimum
node distance on the reference cell (see :func:`dt_non_geometric_factors`),
and :math:`r_D` is the inradius of the cell (see :func:`dt_geometric_factors`).
:returns: a :class:`~meshmode.dof_array.DOFArray` containing a characteristic
lengthscale for each element, at each nodal location.
.. note::
While a prediction of stability is only meaningful in relation to a given
time integrator with a known stability region, the lengthscale provided here
is not intended to be specific to any one time integrator, though the
stability region of standard four-stage, fourth-order Runge-Kutta
methods has been used as a guide. Any concrete time integrator will
likely require scaling of the values returned by this routine.
"""
@memoize_in(dcoll, (characteristic_lengthscales, dd,
"compute_characteristic_lengthscales"))
def _compute_characteristic_lengthscales():
return actx.freeze(
actx.tag(NameHint("char_lscales"),
DOFArray(
actx,
data=tuple(
# Scale each group array of geometric factors by the
# corresponding group non-geometric factor
cng * geo_facts
for cng, geo_facts in zip(
dt_non_geometric_factors(dcoll, dd),
actx.thaw(dt_geometric_factors(dcoll, dd)),
strict=True)))))
return actx.thaw(_compute_characteristic_lengthscales())
@memoize_on_first_arg
def dt_non_geometric_factors(
dcoll: DiscretizationCollection, dd: DOFDesc | None = None
) -> Sequence[float | np.floating]:
r"""Computes the non-geometric scale factors following [Hesthaven_2008]_,
section 6.4, for each element group in the *dd* discretization:
.. math::
c_{ng} = \operatorname{min}\left( \Delta r_i \right),
where :math:`\Delta r_i` denotes the distance between two distinct
nodal points on the reference element.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization if not provided.
:returns: a :class:`list` of :class:`float` values denoting the minimum
node distance on the reference element for each group.
"""
if dd is None:
dd = DD_VOLUME_ALL
discr = dcoll.discr_from_dd(dd)
min_delta_rs: list[np.floating | float] = []
for grp in discr.groups:
assert isinstance(grp, NodalElementGroupBase)
nodes = np.asarray(list(zip(*grp.unit_nodes, strict=True)))
nnodes = grp.nunit_dofs
# NOTE: order 0 elements have 1 node located at the centroid of
# the reference element and is equidistant from each vertex
if grp.order == 0:
assert nnodes == 1
min_delta_rs.append(
np.linalg.norm(
nodes[0] - grp.mesh_el_group.vertex_unit_coordinates()[0]
)
)
else:
min_delta_rs.append(
min(
float(np.linalg.norm(nodes[i] - nodes[j]))
for i in range(nnodes) for j in range(nnodes) if i != j
)
)
return min_delta_rs
# {{{ Mesh size functions
@memoize_on_first_arg
def h_max_from_volume(
dcoll: DiscretizationCollection,
dim: int | None = None,
dd: DOFDesc | None = None) -> Scalar:
"""Returns a (maximum) characteristic length based on the volume of the
elements. This length may not be representative if the elements have very
high aspect ratios.
:arg dim: an integer denoting topological dimension. If *None*, the
spatial dimension specified by
:attr:`grudge.DiscretizationCollection.dim` is used.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization if not provided.
:returns: a scalar denoting the maximum characteristic length.
"""
from grudge.reductions import elementwise_sum, nodal_max
if dd is None:
dd = DD_VOLUME_ALL
dd = as_dofdesc(dd)
if dim is None:
dim = dcoll.dim
ones = dcoll.discr_from_dd(dd).zeros(dcoll._setup_actx) + 1.0
return nodal_max(
dcoll,
dd,
elementwise_sum(dcoll, op.mass(dcoll, dd, ones))
) ** (1.0 / dim)
@memoize_on_first_arg
def h_min_from_volume(
dcoll: DiscretizationCollection,
dim: int | None = None,
dd: DOFDesc | None = None) -> Scalar:
"""Returns a (minimum) characteristic length based on the volume of the
elements. This length may not be representative if the elements have very
high aspect ratios.
:arg dim: an integer denoting topological dimension. If *None*, the
spatial dimension specified by
:attr:`grudge.DiscretizationCollection.dim` is used.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization if not provided.
:returns: a scalar denoting the minimum characteristic length.
"""
from grudge.reductions import elementwise_sum, nodal_min
if dd is None:
dd = DD_VOLUME_ALL
dd = as_dofdesc(dd)
if dim is None:
dim = dcoll.dim
ones = dcoll.discr_from_dd(dd).zeros(dcoll._setup_actx) + 1.0
return nodal_min(
dcoll,
dd,
elementwise_sum(dcoll, op.mass(dcoll, dd, ones))
) ** (1.0 / dim)
def dt_geometric_factors(
dcoll: DiscretizationCollection, dd: DOFDesc | None = None) -> DOFArray:
r"""Computes a geometric scaling factor for each cell following
[Hesthaven_2008]_, section 6.4, defined as the inradius (radius of an
inscribed circle/sphere).
Specifically, the inradius for each element is computed using the following
formula from [Shewchuk_2002]_, Table 1, for simplicial cells
(triangles/tetrahedra):
.. math::
r_D = \frac{d V}{\sum_{i=1}^{N_{faces}} F_i},
where :math:`d` is the topological dimension, :math:`V` is the cell volume,
and :math:`F_i` are the areas of each face of the cell.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization if not provided.
:returns: a frozen :class:`~meshmode.dof_array.DOFArray` containing the
geometric factors for each cell and at each nodal location.
"""
from meshmode.discretization.poly_element import SimplexElementGroupBase
if dd is None:
dd = DD_VOLUME_ALL
actx = dcoll._setup_actx
volm_discr = dcoll.discr_from_dd(dd)
if isinstance(dd.domain_tag, ScalarDomainTag):
raise TypeError("not sensible for scalar domains")
if any(not isinstance(grp, SimplexElementGroupBase)
for grp in volm_discr.groups):
raise NotImplementedError(
"Geometric factors are only implemented for simplex element groups"
)
if volm_discr.dim != volm_discr.ambient_dim:
from warnings import warn
warn("The geometric factor for the characteristic length scale in "
"time step estimation is not necessarily valid for non-volume-"
"filling discretizations. Continuing anyway.", stacklevel=3)
cell_vols: DOFArray = abs(
cast(DOFArray, op.elementwise_integral(
dcoll, dd, volm_discr.zeros(actx) + 1.0
))
)
if dcoll.dim == 1:
# Inscribed "circle" radius is half the cell size
return actx.freeze(cell_vols/2)
dd_face = dd.with_domain_tag(
BoundaryDomainTag(FACE_RESTR_ALL, dd.domain_tag.tag))
face_discr = dcoll.discr_from_dd(dd_face)
# Compute areas of each face
face_areas: DOFArray = abs(
cast(DOFArray, op.elementwise_integral(
dcoll, dd_face, face_discr.zeros(actx) + 1.0
))
)
if actx.supports_nonscalar_broadcasting:
# Compute total surface area of an element by summing over the
# individual face areas
surface_areas = DOFArray(
actx,
data=tuple(
actx.einsum(
"fej->e",
tag_axes(actx, {
0: DiscretizationFaceAxisTag(),
1: DiscretizationElementAxisTag(),
2: DiscretizationDOFAxisTag()
},
face_ae_i.reshape(
vgrp.mesh_el_group.nfaces,
vgrp.nelements,
face_ae_i.shape[-1])),
tagged=(FirstAxisIsElementsTag(),))
for vgrp, face_ae_i in zip(volm_discr.groups, face_areas, strict=True)))
else:
surface_areas = DOFArray(
actx,
data=tuple(
# NOTE: Whenever the array context can't perform nonscalar
# broadcasting, elementwise reductions
# (like `elementwise_integral`) repeat the *same* scalar value of
# the reduction at each degree of freedom. To get a single
# value for the face area (per face),
# we simply average over the nodes, which gives the desired result.
actx.einsum(
"fej->e",
face_ae_i.reshape(
vgrp.mesh_el_group.nfaces,
vgrp.nelements,
face_ae_i.shape[-1]
) / afgrp.nunit_dofs,
tagged=(FirstAxisIsElementsTag(),))
for vgrp, afgrp, face_ae_i in zip(volm_discr.groups,
face_discr.groups,
face_areas, strict=True)
)
)
return actx.freeze(
actx.tag(NameHint(f"dt_geometric_{dd.as_identifier()}"),
DOFArray(actx,
data=tuple(
actx.einsum(
"e,ei->ei",
1/sae_i,
actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i),
tagged=(FirstAxisIsElementsTag(),)) * dcoll.dim
for cv_i, sae_i in zip(cell_vols, surface_areas,
strict=True)))))
# }}}
# vim: foldmethod=marker
grudge/execution.py deleted 100644 → 0
View file @ 90acbf2e
from __future__ import division, absolute_import
__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 six
import numpy as np
import loopy as lp
import pyopencl as cl
import pyopencl.array # noqa
from pytools import memoize_in
import grudge.symbolic.mappers as mappers
from grudge import sym
import logging
logger = logging.getLogger(__name__)
# {{{ exec mapper
class ExecutionMapper(mappers.Evaluator,
mappers.BoundOpMapperMixin,
mappers.LocalOpReducerMixin):
def __init__(self, queue, context, bound_op):
super(ExecutionMapper, self).__init__(context)
self.discr = bound_op.discr
self.bound_op = bound_op
self.queue = queue
def get_discr(self, dd):
qtag = dd.quadrature_tag
if qtag is None:
# FIXME: Remove once proper quadrature support arrives
qtag = sym.QTAG_NONE
if dd.is_volume():
if qtag is not sym.QTAG_NONE:
# FIXME
raise NotImplementedError("quadrature")
return self.discr.volume_discr
elif dd.domain_tag is sym.FRESTR_ALL_FACES:
return self.discr.all_faces_discr(qtag)
elif dd.domain_tag is sym.FRESTR_INTERIOR_FACES:
return self.discr.interior_faces_discr(qtag)
elif dd.is_boundary():
return self.discr.boundary_discr(dd.domain_tag, qtag)
else:
raise ValueError("DOF desc tag not understood: " + str(dd))
# {{{ expression mappings -------------------------------------------------
def map_ones(self, expr):
discr = self.get_discr(expr.dd)
result = discr.empty(self.queue, allocator=self.bound_op.allocator)
result.fill(1)
return result
def map_node_coordinate_component(self, expr):
discr = self.get_discr(expr.dd)
return discr.nodes()[expr.axis].with_queue(self.queue)
def map_boundarize(self, op, field_expr):
return self.discr.boundarize_volume_field(
self.rec(field_expr), tag=op.tag,
kind=self.discr.compute_kind)
def map_grudge_variable(self, expr):
return self.context[expr.name]
def map_call(self, expr):
from pymbolic.primitives import Variable
assert isinstance(expr.function, Variable)
# 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:
func = getattr(np, expr.function.name)
return func(*pars)
def map_nodal_sum(self, op, field_expr):
return cl.array.sum(self.rec(field_expr))
def map_nodal_max(self, op, field_expr):
return cl.array.max(self.rec(field_expr))
def map_nodal_min(self, op, field_expr):
return cl.array.min(self.rec(field_expr))
def map_if(self, expr):
bool_crit = self.rec(expr.condition)
then = self.rec(expr.then)
else_ = self.rec(expr.else_)
result = cl.array.empty_like(then, queue=self.queue,
allocator=self.bound_op.allocator)
cl.array.if_positive(bool_crit, then, else_, out=result,
queue=self.queue)
return result
def map_ref_diff_base(self, op, field_expr):
raise NotImplementedError(
"differentiation should be happening in batched form")
def map_elementwise_linear(self, op, field_expr):
field = self.rec(field_expr)
from grudge.tools import is_zero
if is_zero(field):
return 0
@memoize_in(self.bound_op, "elwise_linear_knl")
def knl():
knl = lp.make_kernel(
"""{[k,i,j]:
0<=k<nelements and
0<=i,j<ndiscr_nodes}""",
"result[k,i] = sum(j, mat[i, j] * vec[k, j])",
default_offset=lp.auto, name="diff")
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
discr = self.get_discr(op.dd_in)
# FIXME: This shouldn't really assume that it's dealing with a volume
# input. What about quadrature? What about boundaries?
result = discr.empty(
queue=self.queue,
dtype=field.dtype, allocator=self.bound_op.allocator)
for grp in discr.groups:
cache_key = "elwise_linear", grp, op, field.dtype
try:
matrix = self.bound_op.operator_data_cache[cache_key]
except KeyError:
matrix = (
cl.array.to_device(
self.queue,
np.asarray(op.matrix(grp), dtype=field.dtype))
.with_queue(None))
self.bound_op.operator_data_cache[cache_key] = matrix
knl()(self.queue, mat=matrix, result=grp.view(result),
vec=grp.view(field))
return result
def map_elementwise_max(self, op, field_expr):
from grudge._internal import perform_elwise_max
field = self.rec(field_expr)
out = self.discr.volume_zeros(dtype=field.dtype)
for eg in self.discr.element_groups:
perform_elwise_max(eg.ranges, field, out)
return out
def map_interpolation(self, op, field_expr):
if op.dd_in.quadrature_tag not in [None, sym.QTAG_NONE]:
raise ValueError("cannot interpolate *from* a quadrature grid")
dd_in = op.dd_in
dd_out = op.dd_out
qtag = dd_out.quadrature_tag
if qtag is None:
# FIXME: Remove once proper quadrature support arrives
qtag = sym.QTAG_NONE
if dd_in.is_volume():
if dd_out.domain_tag is sym.FRESTR_ALL_FACES:
conn = self.discr.all_faces_connection(qtag)
elif dd_out.domain_tag is sym.FRESTR_INTERIOR_FACES:
conn = self.discr.interior_faces_connection(qtag)
elif dd_out.is_boundary():
conn = self.discr.boundary_connection(dd_out.domain_tag, qtag)
else:
raise ValueError("cannot interpolate from volume to: " + str(dd_out))
elif dd_in.domain_tag is sym.FRESTR_INTERIOR_FACES:
if dd_out.domain_tag is sym.FRESTR_ALL_FACES:
conn = self.discr.all_faces_connection(None, qtag)
else:
raise ValueError(
"cannot interpolate from interior faces to: "
+ str(dd_out))
elif dd_in.is_boundary():
if dd_out.domain_tag is sym.FRESTR_ALL_FACES:
conn = self.discr.all_faces_connection(dd_in.domain_tag, qtag)
else:
raise ValueError(
"cannot interpolate from interior faces to: "
+ str(dd_out))
else:
raise ValueError("cannot interpolate from: " + str(dd_in))
return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
def map_opposite_interior_face_swap(self, op, field_expr):
dd = op.dd_in
qtag = dd.quadrature_tag
if qtag is None:
# FIXME: Remove once proper quadrature support arrives
qtag = sym.QTAG_NONE
return self.discr.opposite_face_connection(qtag)(
self.queue, self.rec(field_expr)).with_queue(self.queue)
# {{{ face mass operator
def map_ref_face_mass_operator(self, op, field_expr):
field = self.rec(field_expr)
from grudge.tools import is_zero
if is_zero(field):
return 0
@memoize_in(self.bound_op, "face_mass_knl")
def knl():
knl = lp.make_kernel(
"""{[k,i,f,j]:
0<=k<nelements and
0<=f<nfaces and
0<=i<nvol_nodes and
0<=j<nface_nodes}""",
"result[k,i] = sum(f, sum(j, mat[i, f, j] * vec[f, k, j]))",
default_offset=lp.auto, name="face_mass")
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
qtag = op.dd_in.quadrature_tag
if qtag is None:
# FIXME: Remove once proper quadrature support arrives
qtag = sym.QTAG_NONE
all_faces_conn = self.discr.all_faces_volume_connection(qtag)
all_faces_discr = all_faces_conn.to_discr
vol_discr = all_faces_conn.from_discr
result = vol_discr.empty(
queue=self.queue,
dtype=field.dtype, allocator=self.bound_op.allocator)
assert len(all_faces_discr.groups) == len(vol_discr.groups)
for cgrp, afgrp, volgrp in zip(all_faces_conn.groups,
all_faces_discr.groups, vol_discr.groups):
cache_key = "face_mass", afgrp, op, field.dtype
nfaces = volgrp.mesh_el_group.nfaces
try:
matrix = self.bound_op.operator_data_cache[cache_key]
except KeyError:
matrix = op.matrix(afgrp, volgrp, field.dtype)
matrix = (
cl.array.to_device(self.queue, matrix)
.with_queue(None))
self.bound_op.operator_data_cache[cache_key] = matrix
input_view = afgrp.view(field).reshape(
nfaces, volgrp.nelements, afgrp.nunit_nodes)
knl()(self.queue, mat=matrix, result=volgrp.view(result),
vec=input_view)
return result
# }}}
# }}}
# {{{ code execution functions
def exec_assign(self, insn):
return [(name, self.rec(expr))
for name, expr in zip(insn.names, insn.exprs)], []
def exec_assign_to_discr_scoped(self, insn):
assignments = []
for name, expr in zip(insn.names, insn.exprs):
value = self.rec(expr)
self.discr._discr_scoped_subexpr_name_to_value[name] = value
assignments.append((name, value))
return assignments, []
def exec_assign_from_discr_scoped(self, insn):
return [(insn.name,
self.discr._discr_scoped_subexpr_name_to_value[insn.name])], []
def exec_diff_batch_assign(self, insn):
field = self.rec(insn.field)
repr_op = insn.operators[0]
# FIXME: There's no real reason why differentiation is special,
# execution-wise.
# This should be unified with map_elementwise_linear, which should
# be extended to support batching.
discr = self.get_discr(repr_op.dd_in)
# FIXME: Enable
# assert repr_op.dd_in == repr_op.dd_out
assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
@memoize_in(self.discr, "reference_derivative_knl")
def knl():
knl = lp.make_kernel(
"""{[imatrix,k,i,j]:
0<=imatrix<nmatrices and
0<=k<nelements and
0<=i,j<nunit_nodes}""",
"""
result[imatrix, k, i] = sum(
j, diff_mat[imatrix, i, j] * vec[k, j])
""",
default_offset=lp.auto, name="diff")
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
noperators = len(insn.operators)
result = discr.empty(
queue=self.queue,
dtype=field.dtype, extra_dims=(noperators,),
allocator=self.bound_op.allocator)
for grp in discr.groups:
if grp.nelements == 0:
continue
matrices = repr_op.matrices(grp)
# FIXME: Should transfer matrices to device and cache them
matrices_ary = np.empty((
noperators, grp.nunit_nodes, grp.nunit_nodes))
for i, op in enumerate(insn.operators):
matrices_ary[i] = matrices[op.rst_axis]
knl()(self.queue,
diff_mat=matrices_ary,
result=grp.view(result), vec=grp.view(field))
return [(name, result[i]) for i, name in enumerate(insn.names)], []
# }}}
# }}}
# {{{ bound operator
class BoundOperator(object):
def __init__(self, discr, discr_code, eval_code, debug_flags, allocator=None):
self.discr = discr
self.discr_code = discr_code
self.eval_code = eval_code
self.operator_data_cache = {}
self.debug_flags = debug_flags
self.allocator = allocator
def __str__(self):
sep = 75 * "=" + "\n"
return (
sep
+ "DISCRETIZATION-SCOPE CODE\n"
+ sep
+ str(self.discr_code) + "\n"
+ sep
+ "PER-EVALUATION CODE\n"
+ sep
+ str(self.eval_code))
def __call__(self, queue, **context):
import pyopencl.array as cl_array
def replace_queue(a):
if isinstance(a, cl_array.Array):
return a.with_queue(queue)
else:
return a
from pytools.obj_array import with_object_array_or_scalar
# {{{ discr-scope evaluation
if any(result_var.name not in self.discr._discr_scoped_subexpr_name_to_value
for result_var in self.discr_code.result):
# need to do discr-scope evaluation
discr_eval_context = {}
self.discr_code.execute(ExecutionMapper(queue, discr_eval_context, self))
# }}}
new_context = {}
for name, var in six.iteritems(context):
new_context[name] = with_object_array_or_scalar(replace_queue, var)
return self.eval_code.execute(ExecutionMapper(queue, new_context, self))
# }}}
# {{{ process_sym_operator function
def process_sym_operator(sym_operator, post_bind_mapper=None,
dumper=lambda name, sym_operator: None, mesh=None):
import grudge.symbolic.mappers as mappers
dumper("before-bind", sym_operator)
sym_operator = mappers.OperatorBinder()(sym_operator)
mappers.ErrorChecker(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-cfold", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
# Ordering restriction:
#
# - Must run constant fold before first type inference pass, because zeros,
# while allowed, violate typing constraints (because they can't be assigned
# a unique type), and need to be killed before the type inferrer sees them.
# FIXME: Reenable type inference
# from grudge.symbolic.mappers.type_inference import TypeInferrer
# dumper("before-specializer", sym_operator)
# sym_operator = mappers.OperatorSpecializer(
# TypeInferrer()(sym_operator)
# )(sym_operator)
# Ordering restriction:
#
# - Must run OperatorSpecializer before performing the GlobalToReferenceMapper,
# 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)
# Ordering restriction:
#
# - Must specialize quadrature operators before performing inverse mass
# contraction, because there are no inverse-mass-contracted variants of the
# quadrature operators.
dumper("before-imass", sym_operator)
sym_operator = mappers.InverseMassContractor()(sym_operator)
dumper("before-cfold-2", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
# FIXME: Reenable derivative joiner
# dumper("before-derivative-join", sym_operator)
# sym_operator = mappers.DerivativeJoiner()(sym_operator)
dumper("process-finished", sym_operator)
return sym_operator
# }}}
def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
debug_flags=set(), allocator=None):
# from grudge.symbolic.mappers import QuadratureUpsamplerRemover
# sym_operator = QuadratureUpsamplerRemover(self.quad_min_degrees)(
# sym_operator)
stage = [0]
def dump_optemplate(name, sym_operator):
if "dump_optemplate_stages" in debug_flags:
from grudge.tools import open_unique_debug_file
from grudge.optemplate import pretty
open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt").write(
pretty(sym_operator))
stage[0] += 1
sym_operator = process_sym_operator(sym_operator,
post_bind_mapper=post_bind_mapper,
dumper=dump_optemplate,
mesh=discr.mesh)
from grudge.symbolic.compiler import OperatorCompiler
discr_code, eval_code = OperatorCompiler(discr)(sym_operator)
bound_op = BoundOperator(discr, discr_code, eval_code,
debug_flags=debug_flags, allocator=allocator)
if "dump_op_code" in debug_flags:
from grudge.tools import open_unique_debug_file
open_unique_debug_file("op-code", ".txt").write(
str(bound_op))
if "dump_dataflow_graph" in debug_flags:
bound_op.code.dump_dataflow_graph()
return bound_op
# vim: foldmethod=marker
grudge/sym.py → grudge/geometry/__init__.py
View file @ 49b9b7c7
from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__copyright__ = """
Copyright (C) 2021 University of Illinois Board of Trustees
"""
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
......@@ -22,7 +22,38 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
from grudge.symbolic.primitives import * # noqa
from grudge.symbolic.operators import * # noqa
from grudge.symbolic.tools import pretty # noqa
from grudge.geometry.metrics import (
area_element,
first_fundamental_form,
forward_metric_derivative_mat,
forward_metric_nth_derivative,
inverse_first_fundamental_form,
inverse_metric_derivative_mat,
inverse_surface_metric_derivative,
inverse_surface_metric_derivative_mat,
mv_normal,
normal,
pseudoscalar,
second_fundamental_form,
shape_operator,
summed_curvature,
)
__all__ = (
"area_element",
"first_fundamental_form",
"forward_metric_derivative_mat",
"forward_metric_nth_derivative",
"inverse_first_fundamental_form",
"inverse_metric_derivative_mat",
"inverse_surface_metric_derivative",
"inverse_surface_metric_derivative_mat",
"mv_normal",
"normal",
"pseudoscalar",
"second_fundamental_form",
"shape_operator",
"summed_curvature",
)
grudge/geometry/metrics.py 0 → 100644
View file @ 49b9b7c7
"""
.. currentmodule:: grudge.geometry
Coordinate transformations
--------------------------
.. autofunction:: forward_metric_nth_derivative
.. autofunction:: forward_metric_derivative_mat
.. autofunction:: inverse_metric_derivative_mat
.. autofunction:: first_fundamental_form
.. autofunction:: inverse_first_fundamental_form
Geometry terms
--------------
.. autofunction:: inverse_surface_metric_derivative
.. autofunction:: inverse_surface_metric_derivative_mat
.. autofunction:: pseudoscalar
.. autofunction:: area_element
Normal vectors
--------------
.. autofunction:: mv_normal
.. autofunction:: normal
Curvature tensors
-----------------
.. autofunction:: second_fundamental_form
.. autofunction:: shape_operator
.. autofunction:: summed_curvature
"""
__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 numpy as np
from arraycontext import ArrayContext, register_multivector_as_array_container, tag_axes
from arraycontext.metadata import NameHint
from meshmode.discretization.connection import DirectDiscretizationConnection
from meshmode.dof_array import DOFArray
from meshmode.transform_metadata import (
DiscretizationAmbientDimAxisTag,
DiscretizationTopologicalDimAxisTag,
)
from pymbolic.geometric_algebra import MultiVector
from pytools import memoize_in
from pytools.obj_array import make_obj_array
import grudge.dof_desc as dof_desc
from grudge import DiscretizationCollection
from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE, DOFDesc
register_multivector_as_array_container()
def _geometry_to_quad_if_requested(
dcoll, inner_dd, dd, vec, _use_geoderiv_connection):
def to_quad(vec):
if not dd.uses_quadrature():
return vec
return dcoll.connection_from_dds(inner_dd, dd)(vec)
# FIXME: At least for eager evaluation, this is somewhat inefficient, as
# all element groups vectors are upsampled to the quadrature grid, but then
# only the data for the non-affinely-mapped ones is used.
all_quad_vec = to_quad(vec)
if not _use_geoderiv_connection:
return all_quad_vec
return DOFArray(
vec.array_context,
tuple(
geoderiv_vec_i if megrp.is_affine else all_quad_vec_i
for megrp, geoderiv_vec_i, all_quad_vec_i in zip(
dcoll.discr_from_dd(inner_dd).mesh.groups,
dcoll._base_to_geoderiv_connection(inner_dd)(vec),
all_quad_vec, strict=True)))
# {{{ Metric computations
def forward_metric_nth_derivative(
actx: ArrayContext, dcoll: DiscretizationCollection,
xyz_axis: int, ref_axes: int | tuple[tuple[int, int], ...],
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False) -> DOFArray:
r"""Pointwise metric derivatives representing repeated derivatives of the
physical coordinate enumerated by *xyz_axis*: :math:`x_{\mathrm{xyz\_axis}}`
with respect to the coordinates on the reference element :math:`\xi_i`:
.. math::
D^\alpha x_{\mathrm{xyz\_axis}} =
\frac{\partial^{|\alpha|} x_{\mathrm{xyz\_axis}} }{
\partial \xi_1^{\alpha_1}\cdots \partial \xi_m^{\alpha_m}}
where :math:`\alpha` is a multi-index described by *ref_axes*.
:arg xyz_axis: an integer denoting which physical coordinate to
differentiate.
:arg ref_axes: a :class:`tuple` of tuples indicating indices of
coordinate axes of the reference element to the number of derivatives
which will be taken. For example, the value ``((0, 2), (1, 1))``
indicates taking the second derivative with respect to the first
axis and the first derivative with respect to the second
axis. Each axis must occur only once and the tuple must be sorted
by the axis index.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: If *True*, process returned
:class:`~meshmode.dof_array.DOFArray`\ s through
:meth:`~grudge.DiscretizationCollection._base_to_geoderiv_connection`.
This should be set based on whether the code using the result of this
function is able to make use of these arrays. (This is an internal
argument and not intended for use outside :mod:`grudge`.)
:returns: a :class:`~meshmode.dof_array.DOFArray` containing the pointwise
metric derivative at each nodal coordinate.
"""
if dd is None:
dd = DD_VOLUME_ALL
inner_dd = dd.with_discr_tag(DISCR_TAG_BASE)
if isinstance(ref_axes, int):
ref_axes = ((ref_axes, 1),)
if not isinstance(ref_axes, tuple):
raise ValueError("ref_axes must be a tuple")
if tuple(sorted(ref_axes)) != ref_axes:
raise ValueError("ref_axes must be sorted")
if len(set(ref_axes)) != len(ref_axes):
raise ValueError("ref_axes must not contain an axis more than once")
from pytools import flatten
flat_ref_axes = flatten([rst_axis] * n for rst_axis, n in ref_axes)
from meshmode.discretization import num_reference_derivative
vec = num_reference_derivative(
dcoll.discr_from_dd(inner_dd),
flat_ref_axes,
actx.thaw(dcoll.discr_from_dd(inner_dd).nodes())[xyz_axis]
)
return _geometry_to_quad_if_requested(
dcoll, inner_dd, dd, vec, _use_geoderiv_connection)
def forward_metric_derivative_vector(
actx: ArrayContext, dcoll: DiscretizationCollection,
rst_axis: int | tuple[tuple[int, int], ...],
dd: DOFDesc | None = None, *, _use_geoderiv_connection=False
) -> np.ndarray:
r"""Computes an object array containing the forward metric derivatives
of each physical coordinate.
:arg rst_axis: a :class:`tuple` of tuples indicating indices of
coordinate axes of the reference element to the number of derivatives
which will be taken.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
containing the pointwise metric derivatives at each nodal coordinate.
"""
return make_obj_array([
forward_metric_nth_derivative(
actx, dcoll, i, rst_axis, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection)
for i in range(dcoll.ambient_dim)
]
)
def forward_metric_derivative_mv(
actx: ArrayContext, dcoll: DiscretizationCollection,
rst_axis: int | tuple[tuple[int, int], ...],
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False) -> MultiVector:
r"""Computes a :class:`pymbolic.geometric_algebra.MultiVector` containing
the forward metric derivatives of each physical coordinate.
:arg rst_axis: a :class:`tuple` of tuples indicating indices of
coordinate axes of the reference element to the number of derivatives
which will be taken.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a :class:`pymbolic.geometric_algebra.MultiVector` containing
the forward metric derivatives in each physical coordinate.
"""
return MultiVector(
forward_metric_derivative_vector(
actx, dcoll, rst_axis, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection)
)
def forward_metric_derivative_mat(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False) -> np.ndarray:
r"""Computes the forward metric derivative matrix, also commonly
called the Jacobian matrix, with entries defined as the
forward metric derivatives:
.. math::
J = \left\lbrack
\frac{\partial x_i}{\partial \xi_j}
\right\rbrack_{(0, 0) \leq (i, j) \leq (n, m)}
where :math:`x_1, \dots, x_n` denote the physical coordinates and
:math:`\xi_1, \dots, \xi_m` denote coordinates on the reference element.
Note that, in the case of immersed manifolds, `J` is not necessarily
a square matrix.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a matrix containing the evaluated forward metric derivatives
of each physical coordinate, with respect to each reference coordinate.
"""
ambient_dim = dcoll.ambient_dim
if dd is None:
dd = DD_VOLUME_ALL
dim = dcoll.discr_from_dd(dd).dim
result = np.zeros((ambient_dim, dim), dtype=object)
for j in range(dim):
result[:, j] = forward_metric_derivative_vector(
actx, dcoll, j, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection)
return result
def first_fundamental_form(actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None, *, _use_geoderiv_connection=False
) -> np.ndarray:
r"""Computes the first fundamental form using the Jacobian matrix:
.. math::
\begin{bmatrix}
E & F \\ F & G
\end{bmatrix} :=
\begin{bmatrix}
(\partial_u x)^2 & \partial_u x \partial_v x \\
\partial_u x \partial_v x & (\partial_v x)^2
\end{bmatrix} =
J^T \cdot J
where :math:`u, v` are coordinates on the parameterized surface and
:math:`x(u, v)` defines a parameterized region. Here, :math:`J` is the
corresponding Jacobian matrix.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a matrix containing coefficients of the first fundamental
form.
"""
if dd is None:
dd = DD_VOLUME_ALL
mder = forward_metric_derivative_mat(
actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection)
return mder.T.dot(mder)
def inverse_metric_derivative_mat(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False) -> np.ndarray:
r"""Computes the inverse metric derivative matrix, which is
the inverse of the Jacobian (forward metric derivative) matrix.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a matrix containing the evaluated inverse metric derivatives.
"""
ambient_dim = dcoll.ambient_dim
if dd is None:
dd = DD_VOLUME_ALL
dim = dcoll.discr_from_dd(dd).dim
result = np.zeros((ambient_dim, dim), dtype=object)
for i in range(dim):
for j in range(ambient_dim):
result[i, j] = inverse_metric_derivative(
actx, dcoll, i, j, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection
)
return result
def inverse_first_fundamental_form(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False) -> np.ndarray:
r"""Computes the inverse of the first fundamental form:
.. math::
\begin{bmatrix}
E & F \\ F & G
\end{bmatrix}^{-1} =
\frac{1}{E G - F^2}
\begin{bmatrix}
G & -F \\ -F & E
\end{bmatrix}
where :math:`E, F, G` are coefficients of the first fundamental form.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a matrix containing coefficients of the inverse of the
first fundamental form.
"""
if dd is None:
dd = DD_VOLUME_ALL
dim = dcoll.discr_from_dd(dd).dim
if dcoll.ambient_dim == dim:
inv_mder = inverse_metric_derivative_mat(
actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection)
inv_form1 = inv_mder.dot(inv_mder.T)
else:
form1 = first_fundamental_form(
actx, dcoll, dd=dd, _use_geoderiv_connection=_use_geoderiv_connection)
if dim == 1:
inv_form1 = 1.0 / form1
elif dim == 2:
(E, F), (_, G) = form1 # noqa: N806
inv_form1 = 1.0 / (E * G - F * F) * np.stack(
[make_obj_array([G, -F]),
make_obj_array([-F, E])]
)
else:
raise ValueError(f"{dim}D surfaces not supported" % dim)
return inv_form1
def inverse_metric_derivative(
actx: ArrayContext, dcoll: DiscretizationCollection,
rst_axis: int, xyz_axis: int, dd: DOFDesc,
*, _use_geoderiv_connection=False
) -> DOFArray:
r"""Computes the inverse metric derivative of the physical
coordinate enumerated by *xyz_axis* with respect to the
reference axis *rst_axis*.
:arg rst_axis: an integer denoting the reference coordinate axis.
:arg xyz_axis: an integer denoting the physical coordinate axis.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a :class:`~meshmode.dof_array.DOFArray` containing the
inverse metric derivative at each nodal coordinate.
"""
dim = dcoll.dim
if dim != dcoll.ambient_dim:
raise ValueError(
"Not clear what inverse_metric_derivative means if "
"the derivative matrix is not square!"
)
par_vecs = [
forward_metric_derivative_mv(
actx, dcoll, rst, dd,
_use_geoderiv_connection=_use_geoderiv_connection)
for rst in range(dim)]
# Yay Cramer's rule!
from functools import partial, reduce
from operator import xor as outerprod_op
outerprod = partial(reduce, outerprod_op)
def outprod_with_unit(i, at):
unit_vec = np.zeros(dim)
unit_vec[i] = 1
vecs = par_vecs[:]
vecs[at] = MultiVector(unit_vec)
return outerprod(vecs)
volume_pseudoscalar_inv = outerprod(
forward_metric_derivative_mv(
actx, dcoll, rst_axis, dd,
_use_geoderiv_connection=_use_geoderiv_connection)
for rst_axis in range(dim)
).inv()
result = (outprod_with_unit(xyz_axis, rst_axis)
* volume_pseudoscalar_inv).as_scalar()
return result
def inverse_surface_metric_derivative(
actx: ArrayContext, dcoll: DiscretizationCollection,
rst_axis, xyz_axis, dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False):
r"""Computes the inverse surface metric derivative of the physical
coordinate enumerated by *xyz_axis* with respect to the
reference axis *rst_axis*. These geometric terms are used in the
transformation of physical gradients.
This function does not cache its results.
:arg rst_axis: an integer denoting the reference coordinate axis.
:arg xyz_axis: an integer denoting the physical coordinate axis.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a :class:`~meshmode.dof_array.DOFArray` containing the
inverse metric derivative at each nodal coordinate.
"""
dim = dcoll.dim
ambient_dim = dcoll.ambient_dim
if dd is None:
dd = DD_VOLUME_ALL
dd = dof_desc.as_dofdesc(dd)
if ambient_dim == dim:
result = inverse_metric_derivative(
actx, dcoll, rst_axis, xyz_axis, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection
)
else:
inv_form1 = inverse_first_fundamental_form(actx, dcoll, dd=dd)
result = sum(
inv_form1[rst_axis, d]*forward_metric_nth_derivative(
actx, dcoll, xyz_axis, d, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection
) for d in range(dim))
return result
def inverse_surface_metric_derivative_mat(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None,
*, times_area_element=False, _use_geoderiv_connection=False):
r"""Computes the matrix of inverse surface metric derivatives, indexed by
``(xyz_axis, rst_axis)``. It returns all values of
:func:`inverse_surface_metric_derivative_mat` in cached matrix form.
This function caches its results.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg times_area_element: If *True*, each entry of the matrix is premultiplied
with the value of :func:`area_element`, reflecting the typical use
of the matrix in integrals evaluating weak derivatives.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a :class:`~meshmode.dof_array.DOFArray` containing the
inverse metric derivatives in per-group arrays of shape
``(xyz_dimension, rst_dimension, nelements, ndof)``.
"""
if dd is None:
dd = DD_VOLUME_ALL
dd = dof_desc.as_dofdesc(dd)
@memoize_in(dcoll, (inverse_surface_metric_derivative_mat, dd,
times_area_element, _use_geoderiv_connection))
def _inv_surf_metric_deriv():
if times_area_element:
multiplier = area_element(actx, dcoll, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection)
else:
multiplier = 1
mat = actx.np.stack([
actx.np.stack([
multiplier
* inverse_surface_metric_derivative(
actx, dcoll,
rst_axis, xyz_axis, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection)
for rst_axis in range(dcoll.dim)])
for xyz_axis in range(dcoll.ambient_dim)])
return actx.freeze(
actx.tag(NameHint(f"inv_metric_deriv_{dd.as_identifier()}"),
tag_axes(actx, {
0: DiscretizationAmbientDimAxisTag(),
1: DiscretizationTopologicalDimAxisTag()},
mat)))
return actx.thaw(_inv_surf_metric_deriv())
def _signed_face_ones(
actx: ArrayContext, dcoll: DiscretizationCollection, dd: DOFDesc
) -> DOFArray:
assert dd.is_trace()
# NOTE: ignore quadrature_tags on dd, since we only care about
# the face_id here
dd_base = dd.with_discr_tag(DISCR_TAG_BASE)
all_faces_conn = dcoll.connection_from_dds(
dd_base.untrace(), dd_base
)
assert isinstance(all_faces_conn, DirectDiscretizationConnection)
signed_ones = dcoll.discr_from_dd(dd.with_discr_tag(DISCR_TAG_BASE)).zeros(
actx, dtype=dcoll.real_dtype
) + 1
signed_face_ones_numpy = actx.to_numpy(signed_ones)
for igrp, grp in enumerate(all_faces_conn.groups):
for batch in grp.batches:
assert batch.to_element_face is not None
i = actx.to_numpy(actx.thaw(batch.to_element_indices))
grp_field = signed_face_ones_numpy[igrp].reshape(-1)
grp_field[i] = \
(2.0 * (batch.to_element_face % 2) - 1.0) * grp_field[i]
return actx.from_numpy(signed_face_ones_numpy)
def parametrization_derivative(
actx: ArrayContext, dcoll: DiscretizationCollection, dd: DOFDesc,
*, _use_geoderiv_connection=False) -> MultiVector:
r"""Computes the product of forward metric derivatives spanning the
tangent space with topological dimension *dim*.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a :class:`pymbolic.geometric_algebra.MultiVector` containing
the product of metric derivatives.
"""
if dd is None:
dd = DD_VOLUME_ALL
dim = dcoll.discr_from_dd(dd).dim
if dim == 0:
from pymbolic.geometric_algebra import get_euclidean_space
return MultiVector(
_signed_face_ones(actx, dcoll, dd),
space=get_euclidean_space(dcoll.ambient_dim)
)
from pytools import product
return product(
forward_metric_derivative_mv(
actx, dcoll, rst_axis, dd,
_use_geoderiv_connection=_use_geoderiv_connection)
for rst_axis in range(dim)
)
def pseudoscalar(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None, *, _use_geoderiv_connection=False
) -> MultiVector:
r"""Computes the field of pseudoscalars for the domain/discretization
identified by *dd*.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: A :class:`~pymbolic.geometric_algebra.MultiVector` of
:class:`~meshmode.dof_array.DOFArray`\ s.
"""
if dd is None:
dd = DD_VOLUME_ALL
return parametrization_derivative(
actx, dcoll, dd,
_use_geoderiv_connection=_use_geoderiv_connection).project_max_grade()
def area_element(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False
) -> DOFArray:
r"""Computes the scale factor used to transform integrals from reference
to global space.
This function caches its results.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
Defaults to the base volume discretization.
:arg _use_geoderiv_connection: For internal use. See
:func:`forward_metric_nth_derivative` for an explanation.
:returns: a :class:`~meshmode.dof_array.DOFArray` containing the transformed
volumes for each element.
"""
if dd is None:
dd = DD_VOLUME_ALL
@memoize_in(dcoll, (area_element, dd, _use_geoderiv_connection))
def _area_elements():
result = actx.np.sqrt(
pseudoscalar(
actx, dcoll, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection).norm_squared())
return actx.freeze(
actx.tag(NameHint(f"area_el_{dd.as_identifier()}"), result))
return actx.thaw(_area_elements())
# }}}
# {{{ surface normal vectors
def rel_mv_normal(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None,
*, _use_geoderiv_connection=False) -> MultiVector:
r"""Computes surface normals at each nodal location as a
:class:`~pymbolic.geometric_algebra.MultiVector` relative to the
pseudoscalar of the discretization described by *dd*.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
"""
dd = dof_desc.as_dofdesc(dd)
# NOTE: Don't be tempted to add a sign here. As it is, it produces
# exterior normals for positively oriented curves.
pder = pseudoscalar(
actx, dcoll, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection) \
/ area_element(
actx, dcoll, dd=dd,
_use_geoderiv_connection=_use_geoderiv_connection)
# Dorst Section 3.7.2
return pder << pder.I.inv()
def mv_normal(
actx: ArrayContext, dcoll: DiscretizationCollection, dd: DOFDesc,
*, _use_geoderiv_connection=False
) -> MultiVector:
r"""Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`.
This supports both volume discretizations
(where ambient == topological dimension) and surface discretizations
(where ambient == topological dimension + 1). In the latter case, extra
processing ensures that the returned normal is in the local tangent space
of the element at the point where the normal is being evaluated.
This function caches its results.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
:arg _use_geoderiv_connection: If *True*, process returned
:class:`~meshmode.dof_array.DOFArray`\ s through
:meth:`~grudge.DiscretizationCollection._base_to_geoderiv_connection`.
This should be set based on whether the code using the result of this
function is able to make use of these arrays. The default value agrees
with ``actx.supports_nonscalar_broadcasting``. (This is an internal
argument and not intended for use outside :mod:`grudge`.)
:returns: a :class:`~pymbolic.geometric_algebra.MultiVector`
containing the unit normals.
"""
dd = dof_desc.as_dofdesc(dd)
use_geoderiv_connection = _use_geoderiv_connection
del _use_geoderiv_connection
if use_geoderiv_connection is None:
use_geoderiv_connection = actx.supports_nonscalar_broadcasting
@memoize_in(dcoll, (mv_normal, dd, use_geoderiv_connection))
def _normal():
dim = dcoll.discr_from_dd(dd).dim
ambient_dim = dcoll.ambient_dim
if dim == ambient_dim:
raise ValueError("may only request normals on domains whose topological "
f"dimension ({dim}) differs from "
f"their ambient dimension ({ambient_dim})")
if dim == ambient_dim - 1:
result = rel_mv_normal(
actx, dcoll, dd=dd,
_use_geoderiv_connection=use_geoderiv_connection)
else:
# NOTE: In the case of (d - 2)-dimensional curves, we don't really have
# enough information on the face to decide what an "exterior face normal"
# is (e.g the "normal" to a 1D curve in 3D space is actually a
# "normal plane")
#
# The trick done here is that we take the surface normal, move it to the
# face and then take a cross product with the face tangent to get the
# correct exterior face normal vector.
assert dim == ambient_dim - 2
from grudge.op import project
volm_normal = MultiVector(
project(dcoll, dd.untrace(), dd,
rel_mv_normal(
actx, dcoll,
dd=dd.untrace(),
_use_geoderiv_connection=use_geoderiv_connection
).as_vector(dtype=object))
)
pder = pseudoscalar(
actx, dcoll, dd=dd,
_use_geoderiv_connection=use_geoderiv_connection)
mv = -(volm_normal ^ pder) << volm_normal.I.inv()
result = mv / actx.np.sqrt(mv.norm_squared())
return MultiVector({
bits: actx.freeze(
actx.tag(NameHint(f"normal_{bits}_{dd.as_identifier()}"), ary))
for bits, ary in result.data.items()
}, result.space)
return actx.thaw(_normal())
def normal(actx: ArrayContext, dcoll: DiscretizationCollection, dd: DOFDesc,
*, _use_geoderiv_connection=None):
"""Get the unit normal to the specified surface discretization, *dd*.
This supports both volume discretizations
(where ambient == topological dimension) and surface discretizations
(where ambient == topological dimension + 1). In the latter case, extra
processing ensures that the returned normal is in the local tangent space
of the element at the point where the normal is being evaluated.
This function may be treated as if it caches its results.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
:arg _use_geoderiv_connection: See :func:`mv_normal` for a full description.
(This is an internal argument and not intended for use outside
:mod:`grudge`.)
:returns: an object array of :class:`~meshmode.dof_array.DOFArray`
containing the unit normals at each nodal location.
"""
return mv_normal(
actx, dcoll, dd,
_use_geoderiv_connection=_use_geoderiv_connection
).as_vector(dtype=object)
# }}}
# {{{ Curvature computations
def second_fundamental_form(
actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None) -> np.ndarray:
r"""Computes the second fundamental form:
.. math::
S(x) = \begin{bmatrix}
\partial_{uu} x\cdot n & \partial_{uv} x\cdot n \\
\partial_{uv} x\cdot n & \partial_{vv} x\cdot n
\end{bmatrix}
where :math:`n` is the surface normal, :math:`x(u, v)` defines a parameterized
surface, and :math:`u,v` are coordinates on the parameterized surface.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
:returns: a rank-2 object array describing second fundamental form.
"""
if dd is None:
dd = DD_VOLUME_ALL
dim = dcoll.discr_from_dd(dd).dim
normal = rel_mv_normal(actx, dcoll, dd=dd).as_vector(dtype=object)
if dim == 1:
second_ref_axes: list[tuple[tuple[int, int], ...]] = [((0, 2),)]
elif dim == 2:
second_ref_axes = [((0, 2),), ((0, 1), (1, 1)), ((1, 2),)]
else:
raise ValueError(f"{dim}D surfaces not supported")
from pytools import flatten
form2 = np.empty((dim, dim), dtype=object)
for ref_axes in second_ref_axes:
i, j = flatten([rst_axis] * n for rst_axis, n in ref_axes)
ruv = make_obj_array(
[forward_metric_nth_derivative(actx, dcoll, xyz_axis, ref_axes, dd=dd)
for xyz_axis in range(dcoll.ambient_dim)]
)
form2[i, j] = form2[j, i] = normal.dot(ruv)
return form2
def shape_operator(actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None) -> np.ndarray:
r"""Computes the shape operator (also called the curvature tensor) containing
second order derivatives:
.. math::
C(x) = \begin{bmatrix}
\partial_{uu} x & \partial_{uv} x \\
\partial_{uv} x & \partial_{vv} x
\end{bmatrix}
where :math:`x(u, v)` defines a parameterized surface, and :math:`u,v` are
coordinates on the parameterized surface.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
:returns: a rank-2 object array describing the shape operator.
"""
inv_form1 = inverse_first_fundamental_form(actx, dcoll, dd=dd)
form2 = second_fundamental_form(actx, dcoll, dd=dd)
return -form2.dot(inv_form1)
def summed_curvature(actx: ArrayContext, dcoll: DiscretizationCollection,
dd: DOFDesc | None = None) -> DOFArray | float:
r"""Computes the sum of the principal curvatures:
.. math::
\kappa = \operatorname{Trace}(C(x))
where :math:`x(u, v)` defines a parameterized surface, :math:`u,v` are
coordinates on the parameterized surface, and :math:`C(x)` is the shape
operator.
:arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
:returns: a :class:`~meshmode.dof_array.DOFArray` containing the summed
curvature at each nodal coordinate.
"""
if dd is None:
dd = DD_VOLUME_ALL
dim = dcoll.ambient_dim - 1
if dcoll.ambient_dim == 1:
return 0.0
if dcoll.ambient_dim == dim:
return 0.0
return np.trace(shape_operator(actx, dcoll, dd=dd))
# }}}
# vim: foldmethod=marker
grudge/symbolic/__init__.py → grudge/interpolation.py
View file @ 49b9b7c7
"""Building blocks and mappers for operator expression trees."""
"""
.. currentmodule:: grudge.op
from __future__ import division, absolute_import
Interpolation
-------------
__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
.. autofunction:: interp
"""
__copyright__ = """
Copyright (C) 2021 University of Illinois Board of Trustees
"""
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
......@@ -23,3 +30,19 @@ 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 grudge.discretization import DiscretizationCollection
# FIXME: Should revamp interp and make clear distinctions
# between projection and interpolations.
# Related issue: https://github.com/inducer/grudge/issues/38
def interp(dcoll: DiscretizationCollection, src, tgt, vec):
from warnings import warn
warn("'interp' currently calls to 'project'",
UserWarning, stacklevel=2)
from grudge.projection import project
return project(dcoll, src, tgt, vec)
grudge/models/__init__.py
View file @ 49b9b7c7
"""Base classes for operators."""
from __future__ import division
from __future__ import absolute_import
__copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
__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
......@@ -25,8 +25,10 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
from abc import ABC, abstractmethod
class Operator(object):
class Operator(ABC): # noqa: B024
"""A base class for Discontinuous Galerkin operators.
You may derive your own operators from this class, but, at present
......@@ -34,30 +36,39 @@ class Operator(object):
documentation, to group related classes together in an inheritance
tree.
"""
pass
class TimeDependentOperator(Operator):
"""A base class for time-dependent Discontinuous Galerkin operators.
class HyperbolicOperator(Operator):
"""A base class for hyperbolic Discontinuous Galerkin operators."""
You may derive your own operators from this class, but, at present
this class provides no functionality. Its function is merely as
documentation, to group related classes together in an inheritance
tree.
"""
pass
@abstractmethod
def max_characteristic_velocity(self, actx, **kwargs):
r"""Return a maximum characteristic wavespeed for the operator.
:arg actx: a :class:`arraycontext.ArrayContext`.
:arg \**kwargs: Optional keyword arguments for determining the
max characteristic velocity of the operator.
class HyperbolicOperator(Operator):
"""A base class for hyperbolic Discontinuous Galerkin operators."""
Return a :class:`~meshmode.dof_array.DOFArray` or scalar
representing the (local or global) maximal characteristic velocity of
the operator.
"""
def estimate_rk4_timestep(self, actx, dcoll, **kwargs):
r"""Estimate the largest stable timestep for an RK4 method.
def estimate_rk4_timestep(self, discr, t=None, fields=None):
u"""Estimate the largest stable timestep for an RK4 method.
:arg actx: a :class:`arraycontext.ArrayContext`.
:arg dcoll: a :class:`grudge.discretization.DiscretizationCollection`.
:arg \**kwargs: Optional keyword arguments for determining the
max characteristic velocity of the operator. These are passed
to :meth:`max_characteristic_velocity`.
"""
import grudge.op as op
from grudge.dt_utils import characteristic_lengthscales
wavespeeds = self.max_characteristic_velocity(actx, **kwargs)
local_timesteps = (
characteristic_lengthscales(actx, dcoll) / wavespeeds
)
from grudge.dt_finding import (
dt_non_geometric_factor,
dt_geometric_factor)
return 1 / self.max_eigenvalue(t, fields, discr) \
* (dt_non_geometric_factor(discr)
* dt_geometric_factor(discr))
return op.nodal_min(dcoll, "vol", local_timesteps)