diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py deleted file mode 100644 index ab983681f41ba9b043aa7f18e4365d70cd0a93a2..0000000000000000000000000000000000000000 --- a/examples/wave/wave-min.py +++ /dev/null @@ -1,172 +0,0 @@ -"""Minimal example of a grudge driver.""" - -__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 -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in -all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -THE SOFTWARE. -""" - - -import numpy as np -import pyopencl as cl -import pyopencl.tools as cl_tools - -from arraycontext import PyOpenCLArrayContext, thaw - -from grudge.shortcuts import set_up_rk4 -from grudge import DiscretizationCollection - -from pytools.obj_array import flat_obj_array - -import grudge.op as op - -import logging -logger = logging.getLogger(__name__) - - -def main(ctx_factory, dim=2, order=4, visualize=False): - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext( - queue, - allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)) - ) - - 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) - - dcoll = DiscretizationCollection(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 = thaw(dcoll.nodes(), actx) - 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.models.wave import WeakWaveOperator - from meshmode.mesh import BTAG_ALL, BTAG_NONE - - 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" - ) - - fields = flat_obj_array( - dcoll.zeros(actx), - [dcoll.zeros(actx) for i in range(dcoll.dim)] - ) - - dt_scaling_const = 2/3 - dt = dt_scaling_const * wave_op.estimate_rk4_timestep(dcoll, fields=fields) - - wave_op.check_bc_coverage(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) + 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) - - from time import time - t_last_step = time() - - if visualize: - u = fields[0] - v = fields[1:] - vis.write_vtk_file( - f"fld-wave-min-{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 - - if step % 10 == 0: - logger.info(f"step: {step} t: {time()-t_last_step} " - f"L2: {norm(u=event.state_component[0])}") - if visualize: - vis.write_vtk_file( - f"fld-wave-min-{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 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) diff --git a/examples/wave/wave-op.py b/examples/wave/wave-op.py deleted file mode 100644 index 56a901716efd60f30f813aacc29539f2b40bf2d2..0000000000000000000000000000000000000000 --- a/examples/wave/wave-op.py +++ /dev/null @@ -1,222 +0,0 @@ -"""Minimal example of a grudge driver.""" - -__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 numpy as np -import numpy.linalg as la # noqa -import pyopencl as cl -import pyopencl.tools as cl_tools - -from arraycontext import PyOpenCLArrayContext, thaw - -from pytools.obj_array import flat_obj_array - -from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa - -from grudge.discretization import DiscretizationCollection -from grudge.shortcuts import make_visualizer - -import grudge.op as op - -import logging -logger = logging.getLogger(__name__) - - -# {{{ wave equation bits - -def wave_flux(dcoll, c, w_tpair): - u = w_tpair[0] - v = w_tpair[1:] - - normal = thaw(dcoll.normal(w_tpair.dd), u.int.array_context) - - 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), - ) - - return op.project(dcoll, w_tpair.dd, "all_faces", c*flux_weak) - - -def wave_operator(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) - - return ( - op.inverse_mass( - dcoll, - flat_obj_array( - -c*op.weak_local_div(dcoll, v), - -c*op.weak_local_grad(dcoll, u) - ) - + op.face_mass( - dcoll, - sum( - wave_flux(dcoll, c=c, w_tpair=tpair) - for tpair in op.interior_trace_pairs(dcoll, w) - ) - + wave_flux( - dcoll, c=c, - w_tpair=op.bdry_trace_pair(dcoll, - BTAG_ALL, - interior=dir_bval, - exterior=dir_bc) - ) - ) - ) - ) - -# }}} - - -def rk4_step(y, t, h, f): - k1 = f(t, y) - k2 = f(t+h/2, y + h/2*k1) - k3 = f(t+h/2, y + h/2*k2) - k4 = f(t+h, y + h*k3) - return y + h/6*(k1 + 2*k2 + 2*k3 + k4) - - -def estimate_rk4_timestep(dcoll, c): - from grudge.dt_utils import (dt_non_geometric_factor, - dt_geometric_factors) - - dt_factor = (dt_non_geometric_factor(dcoll) - * op.nodal_min(dcoll, "vol", dt_geometric_factors(dcoll))) - - return dt_factor * (1 / c) - - -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 = thaw(dcoll.nodes(), actx) - 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): - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext( - queue, - allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)) - ) - - 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) - - dcoll = DiscretizationCollection(actx, mesh, order=order) - - fields = flat_obj_array( - bump(actx, dcoll), - [dcoll.zeros(actx) for i in range(dcoll.dim)] - ) - - c = 1 - dt_scaling_const = 0.45 - dt = dt_scaling_const * estimate_rk4_timestep(dcoll, c) - - vis = make_visualizer(dcoll) - - def rhs(t, w): - return wave_operator(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) - - if istep % 10 == 0: - logger.info(f"step: {istep} t: {t} " - f"L2: {op.norm(dcoll, fields[0], 2)} " - f"Linf: {op.norm(dcoll, fields[0], np.inf)} " - f"sol max: {op.nodal_max(dcoll, 'vol', fields[0])} " - f"sol min: {op.nodal_min(dcoll, 'vol', fields[0])}") - if visualize: - vis.write_vtk_file( - f"fld-wave-eager-{istep:04d}.vtu", - [ - ("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 op.norm(dcoll, fields[0], 2) < 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