1-dimensional variable-velocity advection¶
One-dimensional advection with variable velocity¶
Solve the non-conservative variable-coefficient advection equation (also known as the color equation):
\[q_t + u(x) q_x = 0.\]
Here q is the density of some quantity and u(x) is the velocity. The velocity field used is
\[u(x) = 2 + sin(2\pi x).\]
The boundary conditions are periodic. The initial data get stretched and compressed as they move through the fast and slow parts of the velocity field.
Output:¶
Source:¶
#!/usr/bin/env python
# encoding: utf-8
r"""
One-dimensional advection with variable velocity
================================================
Solve the non-conservative variable-coefficient advection equation
(also known as the color equation):
.. math:: q_t + u(x) q_x = 0.
Here q is the density of some quantity and u(x) is the velocity.
The velocity field used is
.. math:: u(x) = 2 + sin(2\pi x).
The boundary conditions are periodic.
The initial data get stretched and compressed as they move through the
fast and slow parts of the velocity field.
"""
import numpy as np
def qinit(state):
# Initial Data parameters
ic = 3
beta = 100.
gamma = 0.
x0 = 0.3
x1 = 0.7
x2 = 0.9
x =state.grid.x.centers
# Gaussian
qg = np.exp(-beta * (x-x0)**2) * np.cos(gamma * (x - x0))
# Step Function
qs = (x > x1) * 1.0 - (x > x2) * 1.0
if ic == 1: state.q[0,:] = qg
elif ic == 2: state.q[0,:] = qs
elif ic == 3: state.q[0,:] = qg + qs
def auxinit(state):
# Initialize petsc Structures for aux
xc=state.grid.x.centers
state.aux[0,:] = np.sin(2.*np.pi*xc)+2
def setup(use_petsc=False,solver_type='classic',kernel_language='Python',outdir='./_output'):
from clawpack import riemann
if use_petsc:
import clawpack.petclaw as pyclaw
else:
from clawpack import pyclaw
if solver_type=='classic':
if kernel_language == 'Fortran':
solver = pyclaw.ClawSolver1D(riemann.advection_color_1D)
elif kernel_language=='Python':
solver = pyclaw.ClawSolver1D(riemann.vc_advection_1D_py.vc_advection_1D)
elif solver_type=='sharpclaw':
if kernel_language == 'Fortran':
solver = pyclaw.SharpClawSolver1D(riemann.advection_color_1D)
elif kernel_language=='Python':
solver = pyclaw.SharpClawSolver1D(riemann.vc_advection_1D_py.vc_advection_1D)
solver.weno_order=weno_order
else: raise Exception('Unrecognized value of solver_type.')
solver.kernel_language = kernel_language
solver.limiters = pyclaw.limiters.tvd.MC
solver.bc_lower[0] = 2
solver.bc_upper[0] = 2
solver.aux_bc_lower[0] = 2
solver.aux_bc_upper[0] = 2
xlower=0.0; xupper=1.0; mx=100
x = pyclaw.Dimension(xlower,xupper,mx,name='x')
domain = pyclaw.Domain(x)
num_aux=1
num_eqn = 1
state = pyclaw.State(domain,num_eqn,num_aux)
qinit(state)
auxinit(state)
claw = pyclaw.Controller()
claw.outdir = outdir
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.tfinal = 1.0
claw.setplot = setplot
claw.keep_copy = True
return claw
#--------------------------
def setplot(plotdata):
#--------------------------
"""
Specify what is to be plotted at each frame.
Input: plotdata, an instance of visclaw.data.ClawPlotData.
Output: a modified version of plotdata.
"""
plotdata.clearfigures() # clear any old figures,axes,items data
# Figure for q[0]
plotfigure = plotdata.new_plotfigure(name='q', figno=1)
# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.ylimits = [-.1,1.1]
plotaxes.title = 'q'
# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
plotitem.plot_var = 0
plotitem.plotstyle = '-o'
plotitem.color = 'b'
plotitem.kwargs = {'linewidth':2,'markersize':5}
return plotdata
if __name__=="__main__":
from clawpack.pyclaw.util import run_app_from_main
output = run_app_from_main(setup,setplot)