# 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
import logging
logger = logging.getLogger(__name__)
def main(write_output=True, allow_features=None, flux_type_arg=1,
bdry_flux_type_arg=None, extra_discr_args={}):
from grudge.mesh.generator import make_cylinder_mesh, make_box_mesh
from grudge.tools import EOCRecorder, to_obj_array
from math import sqrt, pi # noqa
from analytic_solutions import ( # noqa
RealPartAdapter,
SplitComplexAdapter,
CylindricalFieldAdapter,
CylindricalCavityMode,
RectangularWaveguideMode,
RectangularCavityMode)
from grudge.models.em import MaxwellOperator
logging.basicConfig(level=logging.DEBUG)
from grudge.backends import guess_run_context
rcon = guess_run_context(allow_features)
epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
mu0 = 4*pi*1e-7 # N/A**2.
epsilon = 1*epsilon0
mu = 1*mu0
eoc_rec = EOCRecorder()
cylindrical = False
periodic = False
if cylindrical:
R = 1
d = 2
mode = CylindricalCavityMode(m=1, n=1, p=1,
radius=R, height=d,
epsilon=epsilon, mu=mu)
# r_sol = CylindricalFieldAdapter(RealPartAdapter(mode))
# c_sol = SplitComplexAdapter(CylindricalFieldAdapter(mode))
if rcon.is_head_rank:
mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
else:
if periodic:
mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1))
periodicity = (False, False, True)
else:
periodicity = None
mode = RectangularCavityMode(epsilon, mu, (1, 2, 2))
if rcon.is_head_rank:
mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
if rcon.is_head_rank:
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [4, 5, 6]:
#for order in [1,2,3,4,5,6]:
extra_discr_args.setdefault("debug", []).extend([
"cuda_no_plan", "cuda_dump_kernels"])
op = MaxwellOperator(epsilon, mu,
flux_type=flux_type_arg,
bdry_flux_type=bdry_flux_type_arg)
discr = rcon.make_discretization(mesh_data, order=order,
tune_for=op.sym_operator(),
**extra_discr_args)
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "em-%d" % order)
mode.set_time(0)
def get_true_field():
return discr.convert_volume(
to_obj_array(mode(discr)
.real.astype(discr.default_scalar_type).copy()),
kind=discr.compute_kind)
fields = get_true_field()
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper(dtype=discr.default_scalar_type, rcon=rcon)
#from grudge.timestep.dumka3 import Dumka3TimeStepper
#stepper = Dumka3TimeStepper(3, dtype=discr.default_scalar_type, rcon=rcon)
# {{{ diagnostics setup
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
if write_output:
log_file_name = "maxwell-%d.dat" % order
else:
log_file_name = None
logmgr = LogManager(log_file_name, "w", rcon.communicator)
add_run_info(logmgr)
add_general_quantities(logmgr)
add_simulation_quantities(logmgr)
discr.add_instrumentation(logmgr)
stepper.add_instrumentation(logmgr)
from pytools.log import IntervalTimer
vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
logmgr.add_quantity(vis_timer)
from grudge.log import EMFieldGetter, add_em_quantities
field_getter = EMFieldGetter(discr, op, lambda: fields)
add_em_quantities(logmgr, op, field_getter)
logmgr.add_watches(
["step.max", "t_sim.max",
("W_field", "W_el+W_mag"),
"t_step.max"]
)
# }}}
# {{{ timestep loop
rhs = op.bind(discr)
final_time = 0.5e-9
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=final_time, logmgr=logmgr,
max_dt_getter=lambda t: op.estimate_timestep(discr,
stepper=stepper, t=t, fields=fields))
for step, t, dt in step_it:
if step % 50 == 0 and write_output:
sub_timer = vis_timer.start_sub_timer()
e, h = op.split_eh(fields)
visf = vis.make_file("em-%d-%04d" % (order, step))
vis.add_data(
visf,
[
("e",
discr.convert_volume(e, kind="numpy")),
("h",
discr.convert_volume(h, kind="numpy")),
],
time=t, step=step)
visf.close()
sub_timer.stop().submit()
fields = stepper(fields, t, dt, rhs)
mode.set_time(final_time)
eoc_rec.add_data_point(order,
discr.norm(fields-get_true_field()))
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
if rcon.is_head_rank:
print()
print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
# }}}
assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
# {{{ entry points for py.test
from pytools.test import mark_test
@mark_test.long
def test_maxwell_cavities():
main(write_output=False)
@mark_test.long
def test_maxwell_cavities_lf():
main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1)
@mark_test.mpi
@mark_test.long
def test_maxwell_cavities_mpi():
from pytools.mpi import run_with_mpi_ranks
run_with_mpi_ranks(__file__, 2, main,
write_output=False, allow_features=["mpi"])
def test_cuda():
try:
from pycuda.tools import mark_cuda_test
except ImportError:
pass
yield "SP CUDA Maxwell", mark_cuda_test(
do_test_maxwell_cavities_cuda), np.float32
yield "DP CUDA Maxwell", mark_cuda_test(
do_test_maxwell_cavities_cuda), np.float64
def do_test_maxwell_cavities_cuda(dtype):
import pytest # noqa
main(write_output=False, allow_features=["cuda"],
extra_discr_args=dict(
debug=["cuda_no_plan"],
default_scalar_type=dtype,
))
# }}}
# entry point -----------------------------------------------------------------
if __name__ == "__main__":
from pytools.mpi import in_mpi_relaunch
if in_mpi_relaunch():
test_maxwell_cavities_mpi()
else:
do_test_maxwell_cavities_cuda(np.float32)