Newer
Older
# 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
class SteadyShearFlow:
def __init__(self):
self.gamma = 1.5
self.mu = 0.01
self.prandtl = 0.72
self.spec_gas_const = 287.1
def __call__(self, t, x_vec):
# JSH/TW Nodal DG Methods, p.326
rho = numpy.ones_like(x_vec[0])
rho_u = x_vec[1] * x_vec[1]
rho_v = numpy.zeros_like(x_vec[0])
e = (2 * self.mu * x_vec[0] + 10) / (self.gamma - 1) + x_vec[1]**4 / 2
return join_fields(rho, e, rho_u, rho_v)
def properties(self):
return(self.gamma, self.mu, self.prandtl, self.spec_gas_const)
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):
result = discr.convert_boundary(
self(t, discr.get_boundary(tag).nodes.T
.astype(discr.default_scalar_type)),
tag=tag, kind=discr.compute_kind)
return result
def main():
rcon = guess_run_context(
#["cuda"]
)
eoc_rec = EOCRecorder()
def boundary_tagger(vertices, el, face_nr, all_v):
return ["inflow"]
if rcon.is_head_rank:
make_centered_regular_rect_mesh
#mesh = make_rect_mesh((0,0), (10,1), max_area=0.01)
refine = 1
mesh = make_centered_regular_rect_mesh((0,0), (10,1), n=(20,4),
#periodicity=(True, False),
post_refine_factor=refine,
boundary_tagger=boundary_tagger)
mesh_data = rcon.distribute_mesh(mesh)
else:
mesh_data = rcon.receive_mesh()
for order in [3]:
discr = rcon.make_discretization(mesh_data, order=order,
default_scalar_type=numpy.float64)
from grudge.visualization import SiloVisualizer, VtkVisualizer
#vis = VtkVisualizer(discr, rcon, "shearflow-%d" % order)
vis = SiloVisualizer(discr, rcon)
shearflow = SteadyShearFlow()
fields = shearflow.volume_interpolant(0, discr)
gamma, mu, prandtl, spec_gas_const = shearflow.properties()
from grudge.models.gas_dynamics import GasDynamicsOperator
op = GasDynamicsOperator(dimensions=2, gamma=gamma, mu=mu,
prandtl=prandtl, spec_gas_const=spec_gas_const,
bc_inflow=shearflow, bc_outflow=shearflow, bc_noslip=shearflow,
inflow_tag="inflow", outflow_tag="outflow", noslip_tag="noslip")
navierstokes_ex = op.bind(discr)
max_eigval = [0]
def rhs(t, q):
ode_rhs, speed = navierstokes_ex(t, q)
max_eigval[0] = speed
return ode_rhs
# needed to get first estimate of maximum eigenvalue
rhs(0, fields)
if rcon.is_head_rank:
print("---------------------------------------------")
print("order %d" % order)
print("---------------------------------------------")
print("#elements=", len(mesh.elements))
stepper = RK4TimeStepper()
# diagnostics setup ---------------------------------------------------
from pytools.log import LogManager, add_general_quantities, \
add_simulation_quantities, add_run_info
logmgr = LogManager("navierstokes-cpu-%d-%d.dat" % (order, refine),
"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:
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
step_it = times_and_steps(
final_time=0.3,
#max_steps=500,
logmgr=logmgr,
max_dt_getter=lambda t: op.estimate_timestep(discr,
stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
for step, t, dt in step_it:
if step % 10 == 0:
#if False:
visf = vis.make_file("shearflow-%d-%04d" % (order, step))
#true_fields = shearflow.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")),
],
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 = shearflow.volume_interpolant(t, discr)
l2_error = 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"))