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
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 247 additions and 3555 deletions
examples/gas_dynamics/euler/vortex-sources.py deleted 100644 → 0
View file @ d78fe727
# Copyright (C) 2008 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
# modifies isentropic vortex solution so that rho->Arho, P->A^gamma rho^gamma
# this will be analytic solution if appropriate source terms are added
# to the RHS. coded for vel_x=1, vel_y=0
class Vortex:
def __init__(self, beta, gamma, center, velocity, densityA):
self.beta = beta
self.gamma = gamma
self.center = numpy.asarray(center)
self.velocity = numpy.asarray(velocity)
self.densityA = densityA
def __call__(self, t, x_vec):
vortex_loc = self.center + t*self.velocity
# coordinates relative to vortex center
x_rel = x_vec[0] - vortex_loc[0]
y_rel = x_vec[1] - vortex_loc[1]
# Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
# also JSH/TW Nodal DG Methods, p. 209
from math import pi
r = numpy.sqrt(x_rel**2+y_rel**2)
expterm = self.beta*numpy.exp(1-r**2)
u = self.velocity[0] - expterm*y_rel/(2*pi)
v = self.velocity[1] + expterm*x_rel/(2*pi)
rho = self.densityA*(1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
p = rho**self.gamma
e = p/(self.gamma-1) + rho/2*(u**2+v**2)
from grudge.tools import join_fields
return join_fields(rho, e, rho*u, rho*v)
def volume_interpolant(self, t, discr):
return discr.convert_volume(
self(t, discr.nodes.T
.astype(discr.default_scalar_type)),
kind=discr.compute_kind)
def boundary_interpolant(self, t, discr, tag):
return discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T
.astype(discr.default_scalar_type)),
tag=tag, kind=discr.compute_kind)
class SourceTerms:
def __init__(self, beta, gamma, center, velocity, densityA):
self.beta = beta
self.gamma = gamma
self.center = numpy.asarray(center)
self.velocity = numpy.asarray(velocity)
self.densityA = densityA
def __call__(self,t,x_vec,q):
vortex_loc = self.center + t*self.velocity
# coordinates relative to vortex center
x_rel = x_vec[0] - vortex_loc[0]
y_rel = x_vec[1] - vortex_loc[1]
# sources written in terms of A=1.0 solution
# (standard isentropic vortex)
from math import pi
r = numpy.sqrt(x_rel**2+y_rel**2)
expterm = self.beta*numpy.exp(1-r**2)
u = self.velocity[0] - expterm*y_rel/(2*pi)
v = self.velocity[1] + expterm*x_rel/(2*pi)
rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
p = rho**self.gamma
#computed necessary derivatives
expterm_t = 2*expterm*x_rel
expterm_x = -2*expterm*x_rel
expterm_y = -2*expterm*y_rel
u_x = -expterm*y_rel/(2*pi)*(-2*x_rel)
v_y = expterm*x_rel/(2*pi)*(-2*y_rel)
#derivatives for rho (A=1)
facG=self.gamma-1
rho_t = (1/facG)*(1-(facG)/(16*self.gamma*pi**2)*expterm**2)**(1/facG-1)* \
(-facG/(16*self.gamma*pi**2)*2*expterm*expterm_t)
rho_x = (1/facG)*(1-(facG)/(16*self.gamma*pi**2)*expterm**2)**(1/facG-1)* \
(-facG/(16*self.gamma*pi**2)*2*expterm*expterm_x)
rho_y = (1/facG)*(1-(facG)/(16*self.gamma*pi**2)*expterm**2)**(1/facG-1)* \
(-facG/(16*self.gamma*pi**2)*2*expterm*expterm_y)
#derivatives for rho (A=1) to the power of gamma
rho_gamma_t = self.gamma*rho**(self.gamma-1)*rho_t
rho_gamma_x = self.gamma*rho**(self.gamma-1)*rho_x
rho_gamma_y = self.gamma*rho**(self.gamma-1)*rho_y
factorA=self.densityA**self.gamma-self.densityA
#construct source terms
source_rho = x_vec[0]-x_vec[0]
source_e = (factorA/(self.gamma-1))*(rho_gamma_t + self.gamma*(u_x*rho**self.gamma+u*rho_gamma_x)+ \
self.gamma*(v_y*rho**self.gamma+v*rho_gamma_y))
source_rhou = factorA*rho_gamma_x
source_rhov = factorA*rho_gamma_y
from grudge.tools import join_fields
return join_fields(source_rho, source_e, source_rhou, source_rhov, x_vec[0]-x_vec[0])
def volume_interpolant(self,t,q,discr):
return discr.convert_volume(
self(t,discr.nodes.T,q),
kind=discr.compute_kind)
def main(write_output=True):
from grudge.backends import guess_run_context
rcon = guess_run_context(
#["cuda"]
)
gamma = 1.4
# at A=1 we have case of isentropic vortex, source terms
# arise for other values
densityA = 2.0
from grudge.tools import EOCRecorder, to_obj_array
eoc_rec = EOCRecorder()
if rcon.is_head_rank:
from grudge.mesh import \
make_rect_mesh, \
make_centered_regular_rect_mesh
refine = 1
mesh = make_centered_regular_rect_mesh((0,-5), (10,5), n=(9,9),
post_refine_factor=refine)
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [4,5]:
discr = rcon.make_discretization(mesh_data, order=order,
debug=[#"cuda_no_plan",
#"print_op_code"
],
default_scalar_type=numpy.float64)
from grudge.visualization import SiloVisualizer, VtkVisualizer
#vis = VtkVisualizer(discr, rcon, "vortex-%d" % order)
vis = SiloVisualizer(discr, rcon)
vortex = Vortex(beta=5, gamma=gamma,
center=[5,0],
velocity=[1,0], densityA=densityA)
fields = vortex.volume_interpolant(0, discr)
sources=SourceTerms(beta=5, gamma=gamma,
center=[5,0],
velocity=[1,0], densityA=densityA)
from grudge.models.gas_dynamics import (
GasDynamicsOperator, GammaLawEOS)
from grudge.mesh import BTAG_ALL
op = GasDynamicsOperator(dimensions=2,
mu=0.0, prandtl=0.72, spec_gas_const=287.1,
equation_of_state=GammaLawEOS(vortex.gamma),
bc_inflow=vortex, bc_outflow=vortex, bc_noslip=vortex,
inflow_tag=BTAG_ALL, source=sources)
euler_ex = op.bind(discr)
max_eigval = [0]
def rhs(t, q):
ode_rhs, speed = euler_ex(t, q)
max_eigval[0] = speed
return ode_rhs
rhs(0, fields)
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
# limiter setup -------------------------------------------------------
from grudge.models.gas_dynamics import SlopeLimiter1NEuler
limiter = SlopeLimiter1NEuler(discr, gamma, 2, op)
# time stepper --------------------------------------------------------
from grudge.timestep import SSPRK3TimeStepper, RK4TimeStepper
#stepper = SSPRK3TimeStepper(limiter=limiter)
#stepper = SSPRK3TimeStepper()
stepper = RK4TimeStepper()
# diagnostics setup ---------------------------------------------------
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
if write_output:
log_file_name = "euler-%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)
logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
# timestep loop -------------------------------------------------------
t = 0
#fields = limiter(fields)
try:
from grudge.timestep import times_and_steps
step_it = times_and_steps(
final_time=.1,
#max_steps=500,
logmgr=logmgr,
max_dt_getter=lambda t: 0.4*op.estimate_timestep(discr,
stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
for step, t, dt in step_it:
if step % 1 == 0 and write_output:
#if False:
visf = vis.make_file("vortex-%d-%04d" % (order, step))
true_fields = vortex.volume_interpolant(t, discr)
#rhs_fields = rhs(t, fields)
from pyvisfile.silo import DB_VARTYPE_VECTOR
vis.add_data(visf,
[
("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
("e", discr.convert_volume(op.e(fields), kind="numpy")),
("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
("u", discr.convert_volume(op.u(fields), kind="numpy")),
#("true_rho", discr.convert_volume(op.rho(true_fields), kind="numpy")),
#("true_e", discr.convert_volume(op.e(true_fields), kind="numpy")),
#("true_rho_u", discr.convert_volume(op.rho_u(true_fields), kind="numpy")),
#("true_u", discr.convert_volume(op.u(true_fields), kind="numpy")),
#("rhs_rho", discr.convert_volume(op.rho(rhs_fields), kind="numpy")),
#("rhs_e", discr.convert_volume(op.e(rhs_fields), kind="numpy")),
#("rhs_rho_u", discr.convert_volume(op.rho_u(rhs_fields), kind="numpy")),
],
expressions=[
#("diff_rho", "rho-true_rho"),
#("diff_e", "e-true_e"),
#("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
("p", "0.4*(e- 0.5*(rho_u*u))"),
],
time=t, step=step
)
visf.close()
fields = stepper(fields, t, dt, rhs)
true_fields = vortex.volume_interpolant(t, discr)
l2_error = discr.norm(fields-true_fields)
l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields))
l2_error_e = discr.norm(op.e(fields)-op.e(true_fields))
l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields))
l2_error_u = discr.norm(op.u(fields)-op.u(true_fields))
eoc_rec.add_data_point(order, l2_error_rho)
print()
print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
logmgr.set_constant("l2_error", l2_error)
logmgr.set_constant("l2_error_rho", l2_error_rho)
logmgr.set_constant("l2_error_e", l2_error_e)
logmgr.set_constant("l2_error_rhou", l2_error_rhou)
logmgr.set_constant("l2_error_u", l2_error_u)
logmgr.set_constant("refinement", refine)
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
# after order loop
#assert eoc_rec.estimate_order_of_convergence()[0,1] > 6
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
from pytools.test import mark_test
@mark_test.long
def test_euler_vortex():
main(write_output=False)
examples/gas_dynamics/euler/vortex.py deleted 100644 → 0
View file @ d78fe727
# Copyright (C) 2008 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
def main(write_output=True):
from pytools import add_python_path_relative_to_script
add_python_path_relative_to_script("..")
from grudge.backends import guess_run_context
rcon = guess_run_context()
from grudge.tools import EOCRecorder
eoc_rec = EOCRecorder()
if rcon.is_head_rank:
from grudge.mesh.generator import \
make_rect_mesh, \
make_centered_regular_rect_mesh
refine = 4
mesh = make_centered_regular_rect_mesh((0,-5), (10,5), n=(9,9),
post_refine_factor=refine)
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [3, 4, 5]:
from gas_dynamics_initials import Vortex
flow = Vortex()
from grudge.models.gas_dynamics import (
GasDynamicsOperator, PolytropeEOS, GammaLawEOS)
from grudge.mesh import BTAG_ALL
# works equally well for GammaLawEOS
op = GasDynamicsOperator(dimensions=2, mu=flow.mu,
prandtl=flow.prandtl, spec_gas_const=flow.spec_gas_const,
equation_of_state=PolytropeEOS(flow.gamma),
bc_inflow=flow, bc_outflow=flow, bc_noslip=flow,
inflow_tag=BTAG_ALL, source=None)
discr = rcon.make_discretization(mesh_data, order=order,
default_scalar_type=numpy.float64,
quad_min_degrees={
"gasdyn_vol": 3*order,
"gasdyn_face": 3*order,
},
tune_for=op.sym_operator(),
debug=["cuda_no_plan"])
from grudge.visualization import SiloVisualizer, VtkVisualizer
vis = VtkVisualizer(discr, rcon, "vortex-%d" % order)
#vis = SiloVisualizer(discr, rcon)
fields = flow.volume_interpolant(0, discr)
euler_ex = op.bind(discr)
max_eigval = [0]
def rhs(t, q):
ode_rhs, speed = euler_ex(t, q)
max_eigval[0] = speed
return ode_rhs
rhs(0, fields)
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
# limiter ------------------------------------------------------------
from grudge.models.gas_dynamics import SlopeLimiter1NEuler
limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op)
from grudge.timestep.runge_kutta import SSP3TimeStepper
#stepper = SSP3TimeStepper(limiter=limiter)
stepper = SSP3TimeStepper(
vector_primitive_factory=discr.get_vector_primitive_factory())
#from grudge.timestep import RK4TimeStepper
#stepper = RK4TimeStepper()
# diagnostics setup ---------------------------------------------------
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
if write_output:
log_file_name = "euler-%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)
logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
# timestep loop -------------------------------------------------------
try:
final_time = flow.final_time
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, max_eigenvalue=max_eigval[0]))
print("run until t=%g" % final_time)
for step, t, dt in step_it:
if step % 10 == 0 and write_output:
#if False:
visf = vis.make_file("vortex-%d-%04d" % (order, step))
#true_fields = vortex.volume_interpolant(t, discr)
from pyvisfile.silo import DB_VARTYPE_VECTOR
vis.add_data(visf,
[
("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
("e", discr.convert_volume(op.e(fields), kind="numpy")),
("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
("u", discr.convert_volume(op.u(fields), kind="numpy")),
#("true_rho", discr.convert_volume(op.rho(true_fields), kind="numpy")),
#("true_e", discr.convert_volume(op.e(true_fields), kind="numpy")),
#("true_rho_u", discr.convert_volume(op.rho_u(true_fields), kind="numpy")),
#("true_u", discr.convert_volume(op.u(true_fields), kind="numpy")),
#("rhs_rho", discr.convert_volume(op.rho(rhs_fields), kind="numpy")),
#("rhs_e", discr.convert_volume(op.e(rhs_fields), kind="numpy")),
#("rhs_rho_u", discr.convert_volume(op.rho_u(rhs_fields), kind="numpy")),
],
#expressions=[
#("diff_rho", "rho-true_rho"),
#("diff_e", "e-true_e"),
#("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
#("p", "0.4*(e- 0.5*(rho_u*u))"),
#],
time=t, step=step
)
visf.close()
fields = stepper(fields, t, dt, rhs)
#fields = limiter(fields)
assert not numpy.isnan(numpy.sum(fields[0]))
true_fields = flow.volume_interpolant(final_time, discr)
l2_error = discr.norm(fields-true_fields)
l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields))
l2_error_e = discr.norm(op.e(fields)-op.e(true_fields))
l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields))
l2_error_u = discr.norm(op.u(fields)-op.u(true_fields))
eoc_rec.add_data_point(order, l2_error)
print()
print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
logmgr.set_constant("l2_error", l2_error)
logmgr.set_constant("l2_error_rho", l2_error_rho)
logmgr.set_constant("l2_error_e", l2_error_e)
logmgr.set_constant("l2_error_rhou", l2_error_rhou)
logmgr.set_constant("l2_error_u", l2_error_u)
logmgr.set_constant("refinement", refine)
finally:
if write_output:
vis.close()
logmgr.close()
discr.close()
# after order loop
assert eoc_rec.estimate_order_of_convergence()[0,1] > 6
if __name__ == "__main__":
main()
# entry points for py.test ----------------------------------------------------
from pytools.test import mark_test
@mark_test.long
def test_euler_vortex():
main(write_output=False)
examples/gas_dynamics/gas_dynamics_initials.py deleted 100644 → 0
View file @ d78fe727
# Copyright (C) 2008 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
import numpy
import numpy.linalg as la
from six.moves import range
class UniformMachFlow:
def __init__(self, mach=0.1, p=1, rho=1, reynolds=100,
gamma=1.4, prandtl=0.72, char_length=1, spec_gas_const=287.1,
angle_of_attack=None, direction=None, gaussian_pulse_at=None,
pulse_magnitude=0.1):
"""
:param direction: is a vector indicating the direction of the
flow. Only one of angle_of_attack and direction may be
specified. Only the direction, not the magnitude, of
direction is taken into account.
:param angle_of_attack: if not None, specifies the angle of
the flow along the Y axis, where the flow is
directed along the X axis.
"""
if angle_of_attack is not None and direction is not None:
raise ValueError("Only one of angle_of_attack and "
"direction may be specified.")
if angle_of_attack is None and direction is None:
angle_of_attack = 0
if direction is not None:
self.direction = direction/la.norm(direction)
else:
self.direction = None
self.mach = mach
self.p = p
self.rho = rho
self.gamma = gamma
self.prandtl = prandtl
self.reynolds = reynolds
self.length = char_length
self.spec_gas_const = spec_gas_const
self.angle_of_attack = angle_of_attack
self.gaussian_pulse_at = gaussian_pulse_at
self.pulse_magnitude = pulse_magnitude
self.c = (self.gamma * p / rho)**0.5
u = self.velocity = mach * self.c
self.e = p / (self.gamma - 1) + rho / 2 * u**2
if numpy.isinf(self.reynolds):
self.mu = 0
else:
self.mu = u * self.length * rho / self.reynolds
def direction_vector(self, dimensions):
# this must be done here because dimensions is not known above
if self.direction is None:
assert self.angle_of_attack is not None
direction = numpy.zeros(dimensions, dtype=numpy.float64)
direction[0] = numpy.cos(
self.angle_of_attack / 180. * numpy.pi)
direction[1] = numpy.sin(
self.angle_of_attack / 180. * numpy.pi)
return direction
else:
return self.direction
def __call__(self, t, x_vec):
ones = numpy.ones_like(x_vec[0])
rho_field = ones*self.rho
if self.gaussian_pulse_at is not None:
rel_to_pulse = [x_vec[i] - self.gaussian_pulse_at[i]
for i in range(len(x_vec))]
rho_field += self.pulse_magnitude * self.rho * numpy.exp(
- sum(rtp_i**2 for rtp_i in rel_to_pulse)/2)
direction = self.direction_vector(x_vec.shape[0])
from grudge.tools import make_obj_array
u_field = make_obj_array([ones*self.velocity*dir_i
for dir_i in direction])
from grudge.tools import join_fields
return join_fields(rho_field, self.e*ones, self.rho*u_field)
def volume_interpolant(self, t, discr):
return discr.convert_volume(
self(t, discr.nodes.T),
kind=discr.compute_kind,
dtype=discr.default_scalar_type)
def boundary_interpolant(self, t, discr, tag):
return discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T),
tag=tag, kind=discr.compute_kind,
dtype=discr.default_scalar_type)
class Vortex:
def __init__(self):
self.beta = 5
self.gamma = 1.4
self.center = numpy.array([5, 0])
self.velocity = numpy.array([1, 0])
self.mu = 0
self.prandtl = 0.72
self.spec_gas_const = 287.1
def __call__(self, t, x_vec):
vortex_loc = self.center + t*self.velocity
# coordinates relative to vortex center
x_rel = x_vec[0] - vortex_loc[0]
y_rel = x_vec[1] - vortex_loc[1]
# Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
# also JSH/TW Nodal DG Methods, p. 209
from math import pi
r = numpy.sqrt(x_rel**2+y_rel**2)
expterm = self.beta*numpy.exp(1-r**2)
u = self.velocity[0] - expterm*y_rel/(2*pi)
v = self.velocity[1] + expterm*x_rel/(2*pi)
rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
p = rho**self.gamma
e = p/(self.gamma-1) + rho/2*(u**2+v**2)
from grudge.tools import join_fields
return join_fields(rho, e, rho*u, rho*v)
def volume_interpolant(self, t, discr):
return discr.convert_volume(
self(t, discr.nodes.T
.astype(discr.default_scalar_type)),
kind=discr.compute_kind)
def boundary_interpolant(self, t, discr, tag):
return discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T
.astype(discr.default_scalar_type)),
tag=tag, kind=discr.compute_kind)
class Vortex:
def __init__(self):
self.beta = 5
self.gamma = 1.4
self.center = numpy.array([5, 0])
self.velocity = numpy.array([1, 0])
self.final_time = 0.5
self.mu = 0
self.prandtl = 0.72
self.spec_gas_const = 287.1
def __call__(self, t, x_vec):
vortex_loc = self.center + t*self.velocity
# coordinates relative to vortex center
x_rel = x_vec[0] - vortex_loc[0]
y_rel = x_vec[1] - vortex_loc[1]
# Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
# also JSH/TW Nodal DG Methods, p. 209
from math import pi
r = numpy.sqrt(x_rel**2+y_rel**2)
expterm = self.beta*numpy.exp(1-r**2)
u = self.velocity[0] - expterm*y_rel/(2*pi)
v = self.velocity[1] + expterm*x_rel/(2*pi)
rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
p = rho**self.gamma
e = p/(self.gamma-1) + rho/2*(u**2+v**2)
from grudge.tools import join_fields
return join_fields(rho, e, rho*u, rho*v)
def volume_interpolant(self, t, discr):
return discr.convert_volume(
self(t, discr.nodes.T
.astype(discr.default_scalar_type)),
kind=discr.compute_kind)
def boundary_interpolant(self, t, discr, tag):
return discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T
.astype(discr.default_scalar_type)),
tag=tag, kind=discr.compute_kind)
examples/gas_dynamics/lbm-simple.py deleted 100644 → 0
View file @ d78fe727
# Copyright (C) 2011 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 numpy.linalg as la
from six.moves import range
def main(write_output=True, dtype=np.float32):
from grudge.backends import guess_run_context
rcon = guess_run_context()
from grudge.mesh.generator import make_rect_mesh
if rcon.is_head_rank:
h_fac = 1
mesh = make_rect_mesh(a=(0,0),b=(1,1), max_area=h_fac**2*1e-4,
periodicity=(True,True),
subdivisions=(int(70/h_fac), int(70/h_fac)))
from grudge.models.gas_dynamics.lbm import \
D2Q9LBMMethod, LatticeBoltzmannOperator
op = LatticeBoltzmannOperator(
D2Q9LBMMethod(), lbm_delta_t=0.001, nu=1e-4)
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=3,
default_scalar_type=dtype,
debug=["cuda_no_plan"])
from grudge.timestep.runge_kutta import LSRK4TimeStepper
stepper = LSRK4TimeStepper(dtype=dtype,
#vector_primitive_factory=discr.get_vector_primitive_factory()
)
from grudge.visualization import VtkVisualizer
if write_output:
vis = VtkVisualizer(discr, rcon, "fld")
from grudge.data import CompiledExpressionData
def ic_expr(t, x, fields):
from grudge.symbolic import CFunction
from pymbolic.primitives import IfPositive
from pytools.obj_array import make_obj_array
tanh = CFunction("tanh")
sin = CFunction("sin")
rho = 1
u0 = 0.05
w = 0.05
delta = 0.05
from grudge.symbolic.primitives import make_common_subexpression as cse
u = cse(make_obj_array([
IfPositive(x[1]-1/2,
u0*tanh(4*(3/4-x[1])/w),
u0*tanh(4*(x[1]-1/4)/w)),
u0*delta*sin(2*np.pi*(x[0]+1/4))]),
"u")
return make_obj_array([
op.method.f_equilibrium(rho, alpha, u)
for alpha in range(len(op.method))
])
# timestep loop -----------------------------------------------------------
stream_rhs = op.bind_rhs(discr)
collision_update = op.bind(discr, op.collision_update)
get_rho = op.bind(discr, op.rho)
get_rho_u = op.bind(discr, op.rho_u)
f_bar = CompiledExpressionData(ic_expr).volume_interpolant(0, discr)
from grudge.discretization import ExponentialFilterResponseFunction
from grudge.symbolic.operators import FilterOperator
mode_filter = FilterOperator(
ExponentialFilterResponseFunction(min_amplification=0.9, order=4))\
.bind(discr)
final_time = 1000
try:
lbm_dt = op.lbm_delta_t
dg_dt = op.estimate_timestep(discr, stepper=stepper)
print(dg_dt)
dg_steps_per_lbm_step = int(np.ceil(lbm_dt / dg_dt))
dg_dt = lbm_dt / dg_steps_per_lbm_step
lbm_steps = int(final_time // op.lbm_delta_t)
for step in range(lbm_steps):
t = step*lbm_dt
if step % 100 == 0 and write_output:
visf = vis.make_file("fld-%04d" % step)
rho = get_rho(f_bar)
rho_u = get_rho_u(f_bar)
vis.add_data(visf,
[ ("fbar%d" %i,
discr.convert_volume(f_bar_i, "numpy")) for i, f_bar_i in enumerate(f_bar)]+
[
("rho", discr.convert_volume(rho, "numpy")),
("rho_u", discr.convert_volume(rho_u, "numpy")),
],
time=t,
step=step)
visf.close()
print("step=%d, t=%f" % (step, t))
f_bar = collision_update(f_bar)
for substep in range(dg_steps_per_lbm_step):
f_bar = stepper(f_bar, t + substep*dg_dt, dg_dt, stream_rhs)
#f_bar = mode_filter(f_bar)
finally:
if write_output:
vis.close()
discr.close()
if __name__ == "__main__":
main(True)
examples/gas_dynamics/naca.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/gas_dynamics/navierstokes/shearflow.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/gas_dynamics/square.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/gas_dynamics/wing.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/geometry.py
View file @ 49b9b7c7
"""Minimal example of a grudge driver."""
"""Minimal example of viewing geometric quantities."""
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
......@@ -24,37 +25,38 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import numpy as np # noqa
import pyopencl as cl
from grudge import sym, bind, Discretization, shortcuts
import pyopencl.tools as cl_tools
from grudge import geometry, shortcuts
from grudge.array_context import PyOpenCLArrayContext
from grudge.discretization import make_discretization_collection
def main(write_output=True):
def main(write_output: bool = True) -> None:
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
from meshmode.mesh.generation import generate_warped_rect_mesh
mesh = generate_warped_rect_mesh(dim=2, order=4, n=6)
allocator = cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
actx = PyOpenCLArrayContext(queue, allocator=allocator)
discr = Discretization(cl_ctx, mesh, order=4)
from meshmode.mesh import BTAG_ALL
from meshmode.mesh.generation import generate_warped_rect_mesh
sym_op = sym.normal(sym.BTAG_ALL, mesh.dim)
#sym_op = sym.nodes(mesh.dim, where=sym.BTAG_ALL)
print(sym.pretty(sym_op))
op = bind(discr, sym_op)
print()
print(op.code)
mesh = generate_warped_rect_mesh(dim=2, order=4, nelements_side=6)
dcoll = make_discretization_collection(actx, mesh, order=4)
vec = op(queue)
nodes = actx.thaw(dcoll.nodes())
bdry_nodes = actx.thaw(dcoll.nodes(dd=BTAG_ALL))
bdry_normals = geometry.normal(actx, dcoll, dd=BTAG_ALL)
vis = shortcuts.make_visualizer(discr, 4)
vis.write_vtk_file("geo.vtu", [
])
if write_output:
vis = shortcuts.make_visualizer(dcoll)
vis.write_vtk_file("geo.vtu", [("nodes", nodes)])
bvis = shortcuts.make_boundary_visualizer(discr, 4)
bvis.write_vtk_file("bgeo.vtu", [
("normals", vec)
])
bvis = shortcuts.make_boundary_visualizer(dcoll)
bvis.write_vtk_file("bgeo.vtu", [("bdry normals", bdry_normals),
("bdry nodes", bdry_nodes)])
if __name__ == "__main__":
......
examples/heat/heat.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/hello-grudge.py 0 → 100644
View file @ 49b9b7c7
This diff is collapsed.
examples/maxwell/.gitignore deleted 100644 → 0
View file @ d78fe727
bessel_zeros.py
2d_cavity
examples/maxwell/cavities.py
View file @ 49b9b7c7
This diff is collapsed.
examples/maxwell/dipole.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/maxwell/generate-bessel-zeros.py deleted 100644 → 0
View file @ d78fe727
from __future__ import absolute_import
from six.moves import range
def main():
import scipy.special
maxnu = 10
n_zeros = 20
zeros = []
for n in range(0,maxnu+1):
zeros.append(list(scipy.special.jn_zeros(n, n_zeros)))
outf = open("bessel_zeros.py", "w").write("bessel_zeros = %s" % zeros)
if __name__ == "__main__":
main()
examples/maxwell/inhom-cavity.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/maxwell/inhomogeneous_waveguide.mac deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/maxwell/maxwell-2d.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/maxwell/maxwell-pml.py deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.
examples/maxwell/notes.tm deleted 100644 → 0
View file @ d78fe727
This diff is collapsed.