Skip to content

Commit

Permalink
Merge pull request #449 from mandli/fix-up-quadrants-example
Browse files Browse the repository at this point in the history
Add PyClaw Command Line Interface to Euler quadrants example
  • Loading branch information
ketch committed Jul 24, 2014
2 parents a94b1a6 + 5648d5a commit 6ed9a60
Showing 1 changed file with 122 additions and 35 deletions.
157 changes: 122 additions & 35 deletions examples/euler_2d/quadrants.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
#!/usr/bin/env python
# encoding: utf-8
"""
2D Euler Riemann problem
================================
Euler 2D Quadrants example
==========================
Solve the Euler equations of compressible fluid dynamics:
Simple example solving the Euler equations of compressible fluid dynamics:
.. math::
\rho_t + (\rho u)_x + (\rho v)_y & = 0 \\
Expand All @@ -15,37 +15,124 @@
Here :math:`\rho` is the density, (u,v) is the velocity, and E is the total energy.
The initial condition is one of the 2D Riemann problems from the paper of
Liska and Wendroff.
"""

from clawpack import pyclaw
from clawpack import riemann

solver = pyclaw.ClawSolver2D(riemann.euler_4wave_2D)
solver.all_bcs = pyclaw.BC.extrap

domain = pyclaw.Domain([0.,0.],[1.,1.],[100,100])
solution = pyclaw.Solution(solver.num_eqn,domain)
gamma = 1.4
solution.problem_data['gamma'] = gamma
solver.dimensional_split = False
solver.transverse_waves = 2

# Set initial data
xx,yy = domain.grid.p_centers
l = xx<0.8; r = xx>=0.8; b = yy<0.8; t = yy>=0.8
solution.q[0,...] = 1.5*r*t + 0.532258064516129*l*t + 0.137992831541219*l*b + 0.532258064516129*r*b
u = 0.*r*t + 1.206045378311055*l*t + 1.206045378311055*l*b + 0.*r*b
v = 0.*r*t + 0.*l*t + 1.206045378311055*l*b + 1.206045378311055*r*b
p = 1.5*r*t + 0.3*l*t + 0.029032258064516*l*b + 0.3*r*b
solution.q[1,...] = solution.q[0,...] * u
solution.q[2,...] = solution.q[0,...] * v
solution.q[3,...] = 0.5*solution.q[0,...]*(u**2+v**2) + p/(gamma-1.)

claw = pyclaw.Controller()
claw.tfinal = 0.8
claw.solution = solution
claw.solver = solver

status = claw.run()

#pyclaw.plot.interactive_plot()
def setplot(plotdata):
r"""Plotting settings
Should plot two figures both of density.
"""

from clawpack.visclaw import colormaps

plotdata.clearfigures() # clear any old figures,axes,items data

# Figure for density - pcolor
plotfigure = plotdata.new_plotfigure(name='Density', figno=0)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = 'auto'
plotaxes.ylimits = 'auto'
plotaxes.scaled = True
plotaxes.title = 'Density'

# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
plotitem.plot_var = 0
plotitem.pcolor_cmap = colormaps.yellow_red_blue
plotitem.pcolor_cmin = 0.
plotitem.pcolor_cmax = 2.
plotitem.add_colorbar = True

# Figure for density - Schlieren
plotfigure = plotdata.new_plotfigure(name='Schlieren', figno=1)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = 'auto'
plotaxes.ylimits = 'auto'
plotaxes.title = 'Density'
plotaxes.scaled = True # so aspect ratio is 1

# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='2d_schlieren')
plotitem.schlieren_cmin = 0.0
plotitem.schlieren_cmax = 1.0
plotitem.plot_var = 0
plotitem.add_colorbar = False

return plotdata


def setup(use_petsc=False):
r"""Setup for Euler 2D quadrants example
Simple example solving the Euler equations of compressible fluid dynamics:
.. math::
\rho_t + (\rho u)_x + (\rho v)_y & = 0 \\
(\rho u)_t + (\rho u^2 + p)_x + (\rho uv)_y & = 0 \\
(\rho v)_t + (\rho uv)_x + (\rho v^2 + p)_y & = 0 \\
E_t + (u (E + p) )_x + (v (E + p))_y & = 0.
Here :math:`\rho` is the density, (u,v) is the velocity, and E is the total energy.
The initial condition is one of the 2D Riemann problems from the paper of
Liska and Wendroff.
Currently the only setup option is the *use_petsc* argument.
"""

if use_petsc:
import clawpack.petclaw as pyclaw
else:
from clawpack import pyclaw
from clawpack import riemann

solver = pyclaw.ClawSolver2D(riemann.euler_4wave_2D)
solver.all_bcs = pyclaw.BC.extrap

domain = pyclaw.Domain([0.,0.],[1.,1.],[100,100])
solution = pyclaw.Solution(solver.num_eqn,domain)
gamma = 1.4
solution.problem_data['gamma'] = gamma
solver.dimensional_split = False
solver.transverse_waves = 2

# Set initial data
xx, yy = domain.grid.p_centers
l = xx < 0.8
r = xx >= 0.8
b = yy < 0.8
t = yy >= 0.8
solution.q[0,...] = 1.5 * r * t + 0.532258064516129 * l * t \
+ 0.137992831541219 * l * b \
+ 0.532258064516129 * r * b
u = 0.0 * r * t + 1.206045378311055 * l * t \
+ 1.206045378311055 * l * b \
+ 0.0 * r * b
v = 0.0 * r * t + 0.0 * l * t \
+ 1.206045378311055 * l * b \
+ 1.206045378311055 * r * b
p = 1.5 * r * t + 0.3 * l * t + 0.029032258064516 * l * b + 0.3 * r * b
solution.q[1,...] = solution.q[0, ...] * u
solution.q[2,...] = solution.q[0, ...] * v
solution.q[3,...] = 0.5 * solution.q[0,...]*(u**2 + v**2) + p / (gamma - 1.0)

claw = pyclaw.Controller()
claw.tfinal = 0.8
claw.solution = solution
claw.solver = solver

claw.output_format = 'ascii'
claw.outdir = "./_output"
claw.setplot = setplot

return claw

if __name__ == "__main__":
from clawpack.pyclaw.util import run_app_from_main
output = run_app_from_main(setup, setplot)

0 comments on commit 6ed9a60

Please sign in to comment.