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
unported-examples/maxwell/maxwell-pml.py deleted 100644 → 0
View file @ 9af966c7
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 logpyle 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 logpyle 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)
unported-examples/maxwell/notes.tm deleted 100644 → 0
View file @ 9af966c7
<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
unported-examples/poisson/helmholtz.py deleted 100644 → 0
View file @ 9af966c7
__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.
"""
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()
unported-examples/poisson/poisson.py deleted 100644 → 0
View file @ 9af966c7
__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.
"""
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)
unported-examples/wave/wiggly.py deleted 100644 → 0
View file @ 9af966c7
"""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.FunctionSymbol("sin")(source_omega*sym.ScalarParameter("t"))
* sym.FunctionSymbol("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 logpyle 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()